Hostname: page-component-78c5997874-fbnjt Total loading time: 0 Render date: 2024-11-03T20:01:52.011Z Has data issue: false hasContentIssue false

NUMERICAL ANALYSIS OF APPARATUS-INDUCED DISPERSION FOR DENSITY-DEPENDENT SOLUTE TRANSPORT IN POROUS MEDIA

Published online by Cambridge University Press:  31 August 2023

H. ZHANG*
Affiliation:
School of Engineering and Built Environment, Griffith University, Gold Coast Campus, QLD 4222, Australia
D. A. BARRY
Affiliation:
Laboratoire de technologie écologique, Institut d’ingénierie de l’environnement, Faculté de l’environnement naturel, architectural et construit (ENAC), Station 2, Ecole polytechnique fédérale de Lausanne (EPFL), 1015 Lausanne, Switzerland; e-mail: [email protected]
B. SEYMOUR
Affiliation:
Department of Mathematics, University of British Columbia, Vancouver V6T 1Z2, Canada; e-mail: [email protected]
G. HOCKING
Affiliation:
Mathematics and Statistics, Murdoch University, Perth, WA 6150, Australia; e-mail: [email protected]
Rights & Permissions [Opens in a new window]

Abstract

The effects of apparatus-induced dispersion on nonuniform, density-dependent flow in a cylindrical soil column were investigated using a finite-element model. To validate the model, the results with an analytical solution and laboratory column test data were analysed. The model simulations confirmed that flow nonuniformities induced by the apparatus are dissipated within the column when the distance to the apparatus outlet exceeds $3R/2$, where R represents the radius of the cylindrical column. Furthermore, the simulations revealed that convergent flow in the vicinity of the outlet introduces additional hydrodynamic dispersion in the soil column apparatus. However, this effect is minimal in the region where the column height exceeds $3R/2$. Additionally, it is found that an increase in the solution density gradient during the solute breakthrough period led to a decrease in flow velocity, which stabilized the flow and ultimately reduced dispersive mixing. Overall, this study provides insights into the behaviour of apparatus-induced dispersion in nonuniform, density-dependent flow within a cylindrical soil column, shedding light on the dynamics and mitigation of flow nonuniformities and dispersive mixing phenomena.

Type
Research Article
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2023. Published by Cambridge University Press on behalf of Australian Mathematical Publishing Association Inc.

1 Introduction

Miscible, density-dependent flow in porous media is important in many engineering fields, such as groundwater contamination remediation, petroleum extraction and sea water intrusion [Reference Al-Dousari, Garrouch and Guedouar2, Reference Holm15, Reference Koch and Zhang17, Reference Nick, Schotting, Gutierrez-Neri and Johannsen18, Reference Quintard, Bertin and Prouvost22, Reference Tran, Bajracharya and Barry26, Reference Watson and Barry29, Reference Zhang and Hocking34Reference Zhang, Hocking and Seymour36]. Laboratory-scale column experiments are commonly carried out to obtain estimates of the properties of the soil and the transport parameters that characterize contaminant migration through it. Probably, the most straightforward way to perform experiments on nonreactive solute transport in laboratory columns is to first establish a steady-state flow of water through the column, then to switch the influent to a solution with a known chemical composition, leaving the flow rate unchanged (see, for example, [Reference Boon, Bijeljic, Niu and Krevor9, Reference Yong, Mohamed and Warkentin33]). The effluent can be sampled to measure the concentration or visualized with an X-ray medical CT-scanner, and the data obtained can used in conjunction with a model to provide estimates of soil and transport parameters [Reference Bajracharya and Barry3]. This experimental design allows the use of one-dimensional or two-dimensional solute transport models [Reference Bajracharya, Tran and Barry4, Reference Barry6, Reference Chen, Liu, Liang, Liu and Lin10, Reference Kanwar, Johnson and Kirkham16, Reference Overman, Chu and Le20, Reference Rao, Rolston, Jessup and Davidson23, Reference Tran, Bajracharya and Barry26, Reference Tran, Barry and Bajracharya27, Reference Wissmeier, Barry and Phillips32], which simplifies the analysis considerably.

A key parameter is the dispersivity of the solute(s) under consideration. If a one-dimensional solute transport model is used, then the water flow is assumed to be spatially and temporally uniform. For flow through an unconsolidated porous medium, the relationship between flow rate and dispersivity is well known for nonreactive solutes [Reference Bear8, Reference Pfannkuch21]. However, any additional dispersion caused by the apparatus is not included in these estimates. Apparatus-induced dispersion can arise if components of the apparatus cause additional mixing and spreading of the solute front, particularly at the ends of the column apparatus. Apparatus-induced dispersion can also arise if the column apparatus causes nonuniform flow to develop within the column, even when the test soil is homogeneous and uniformly packed.

It is common in column experiments for flow to enter and leave the column through a hole much smaller than the column itself. This occurs for a simple practical reason: the influent or effluent flow occurs in tubes with diameters much less than the soil column diameter. This case was analysed by Barry [Reference Barry6], who developed a general analytical solution for steady flow in such a column, subject to arbitrary head or flux boundary conditions. Although the results of [Reference Barry6] are useful for the analysis of nonuniform flow in a column, they do not consider (i) the effect of variable density flow or (ii) the impact of the flow field on solute dispersion within the column. The purpose of this paper is to quantify both of these aspects of solute transport in soil column experiments.

In the present study, the problem is investigated numerically using a finite-element modelling package, COMSOL Multiphysics [11], which has been employed to solve various partial differential equations in engineering fields [Reference Al-Ali, Hocking, Farrow and Zhang1, Reference Barletta, Zanchini, Lazzari and Terenzi5, Reference Cristea, Bagiu and Agachi12, Reference Sahu and Bhattacharyay24, Reference Wissmeier and Barry31]. The impacts of the dimensions of the apparatus and the soil properties on the dispersion of density-dependent solute transport in porous media are analysed numerically.

2 Problem formulation

2.1 Governing equations

We consider a typical experiment carried out in cylindrical soil columns [Reference Watson, Barry, Schotting and Hassanizadeh30] containing uniformly packed isotropic soil, with influent and effluent orifices aligned along the column centreline. Such experiments produce axisymmetric flow in a cylinder containing soil with porosity n and intrinsic permeability $\kappa $ , as shown in Figure 1. The column radius is R, with height is H and outlet orifice radius is $r_n$ . The origin of the cylindrical coordinate system is located at the centre of the bottom of the column $( r,z ) =( 0,0 )$ . The transported fluid has density $\rho $ , dynamic viscosity $\mu $ and salt mass fraction $\omega $ .

The flow and the solute satisfy Darcy’s law and the mass conservation equations, which are summarized as follows. (See [Reference Watson, Barry, Schotting and Hassanizadeh30].)

Darcy’s law and Fick’s law

(2.1) $$ \begin{align} &\mathbf{q}=-\frac{k}{\mu}( \nabla p+\rho \mathbf{g} ), \kern210pt\end{align} $$
(2.2) $$ \begin{align} &\mathbf{J}=-\rho \boldsymbol{D}\cdot \nabla \omega, \nonumber \\ &\text{where }\mathbf{D}=( D_{ij} ) =\bigg( n( D_m+\alpha _T| \boldsymbol{q} | ) \delta _{ij}+n( \alpha _L-\alpha _T ) \frac{q_iq_j}{| \boldsymbol{q} |} \bigg),\quad i=1, 2, 3. \end{align} $$

Here, $\mathbf {q}$ denotes the specific discharge, p is the fluid pressure, $\mathbf {g}=(0,0,g)$ is the gravitational acceleration, n is the porosity, $\mathbf {h}$ is the dispersive mass flux and $\mathbf {D}=(D_{ij})$ the hydrodynamic dispersion tensor [Reference Bear8], where $D_m$ is the mass diffusivity and $\alpha _L$ and $\alpha _T$ are the longitudinal and transversal dispersivity, respectively.

Figure 1 Two different soil column study configuration diagrams: (a) experimental set-up of Watson et al. [Reference Watson, Barry, Schotting and Hassanizadeh30]; and (b) set-up of the flow convergence study by Barry [Reference Barry6].

Conservation of mass of the water and solute

(2.3) $$ \begin{align} \frac{\partial ( n\rho )}{\partial t}+\nabla \cdot (\rho \mathbf{q})=0, \end{align} $$
(2.4) $$ \begin{align} \frac{\partial ( n\rho \omega )}{\partial t}+\nabla \cdot (\rho \omega \mathbf{q})+\nabla \cdot \mathbf{J}=0, \end{align} $$

where $\rho $ is the fluid density, $\mu $ is the dynamic viscosity of fluid and $\omega $ is the fractional solute concentration. The saltwater density and viscosity, respectively, are given by,

(2.5) $$ \begin{align} \rho &=\rho _0\exp ( \gamma \omega ) , \end{align} $$
(2.6) $$ \begin{align} \mu &=\mu _0(1+1.85\omega -4.10\omega ^2+44.5\omega ^3). \end{align} $$

In the equations of state (2.5) and (2.6) for brine, the constant $\gamma = 0.6923$ [Reference Hassanizadeh and Leijnse14], while, for the freshwater reference state, $\rho _0 = 998.23$ kg/m $^3$ and $\mu _0 = 0.001 \ \mathrm {N_S}$ /m $^2$ at 20 $^{\circ }$ C and 1 atm.

The full set of equations is coupled between density $\rho $ , viscosity $\mu $ , specific discharge $\mathbf {q}$ and dispersive flux $\mathbf {J}$ . All variables with dimensions of length are scaled with R, the radius of the column, so that

$$ \begin{align*}z^*=\frac{z}{R},\quad r^*=\frac{r}{R},\quad {r_n}^*=\frac{r_n}{R},\quad H^*=\frac{H}{R}.\end{align*} $$

The other nondimensional variables are

$$ \begin{align*}\rho ^*=\frac{\rho}{\rho _0}, \quad \mu ^*=\frac{\mu}{\mu _0}, \quad p^*=\frac{p}{\rho _0gR}\quad \text{and}\quad\mathbf{g}^*=\frac{\mathbf{g}}{g}.\end{align*} $$

The reference hydraulic conductivity, $K=\kappa \rho _0g/\mu _0$ with $\kappa $ as the permeability, is used to scale time $t^*=tK/nR$ , the discharge flux $\mathbf {q}^*=\mathbf {q}/q_{in}$ , the hydrodynamic dispersion tensor $\boldsymbol {D}^*=\boldsymbol {D}/Rq_{in}$ and the dispersive mass flux $\mathbf {J}^*=\mathbf {J}/\rho _0Rq_{in}$ , where $q_{in}$ is the discharge flux for $0<r<R\text { and }z=0$ . Then equations (2.1)–(2.6) can be rewritten as follows.

Darcy’s law and Fick’s law

(2.7) $$ \begin{align} \mathbf{q}^*&=-\frac{K}{q_{in}}\frac{1}{\mu ^*}( \nabla p^*+\rho ^*\mathbf{g}^* ), \end{align} $$
(2.8) $$ \begin{align}\mathbf{J}^*&=-\rho ^*\mathbf{D}^*\cdot \nabla \omega. \kern35pt\end{align} $$

Conservation of mass

(2.9) $$ \begin{align} \frac{\partial ( \rho ^* )}{\partial t^*}+\nabla \cdot ( \rho ^*\mathbf{q}^* ) =0, \end{align} $$
(2.10) $$ \begin{align} \frac{\partial ( \rho ^*\omega )}{\partial t^*}+\nabla \cdot ( \rho ^*\omega \mathbf{q}^* ) +\nabla \cdot \mathbf{J}^*=0, \end{align} $$

with

(2.11) $$ \begin{align} \rho ^*&=\exp ( \gamma \omega ), \end{align} $$
(2.12) $$ \begin{align} \mu ^*&=1+1.85\omega -4.10\omega ^2+44.5\omega ^3. \end{align} $$

The solute transport depends on the flow, as the velocity $\mathbf {q}^*$ determines both advective and dispersive transport. The density $\rho ^*( \omega )$ and viscosity $\mu ^*( \omega )$ of the solution are coupled back to the flow equation through their dependencies on the mass fraction of the solution, which produces two nonlinear diffusive partial differential equations for the variables $\omega $ and $p^*$ as functions of $z*$ , $r^*$ and $t^*$ .

2.2 Initial and boundary conditions

Following the experiment set-up of [Reference Watson and Barry29], it is assumed that the fluid enters the column from the bottom and exits from the orifice at the top. The initial conditions for the simulations are that the fluid in the column is under hydrostatic condition and there is no solute, so that

(2.13) $$ \begin{align} \begin{cases} p^*=z^*\\ \omega ^*=0\\ \end{cases} \text{for } t^*=0, 0\leq z^*\leq H^*\text{ and } 0\leq r^*\leq 1. \end{align} $$

Measured flow rate data at the influent boundary and measured pressure data at the effluent boundary are used for the boundary conditions applied to the fluid mass balance matrix equation for comparison with the experiment results. At the inlet, a first-type of concentration-type input boundary condition (see, for example, [Reference Barry and Sposito7, Reference Van Genuchten and Parker28]) was used, where the measured flow rates $q^*_{in}$ and solute concentration $\omega _{in}$ are specified: that is,

(2.14) $$ \begin{align} \mathbf{q}^*=( 0,q_{in}^{*} ) ,\quad \omega =\omega _{in}\quad \text{for } 0\le r^*<1, \quad z^*=0\ \text{and all } t^*. \end{align} $$

At the orifice, the most realistic assumption for the finite column apparatus is the Danckwerts boundary condition [Reference Danckwerts13, Reference Schwartz, McInnes, Juo, Wilding and Reddell25, Reference Van Genuchten and Parker28] for the solute concentration and constant pressure for the fluid: that is,

(2.15) $$ \begin{align} p^*=0,\quad \frac{\partial \omega}{\partial z^*}=0\quad\text{for } 0\le r^*<r_{n}^{*},\quad z^*=H^*\text{ at all } t^*. \end{align} $$

In the present study, COMSOL Multiphysics modules of Darcy’s law and solute transport are coupled for a numerical study on the above apparatus-induced dispersion in porous media.

3 Simulations and analysis

3.1 Convergent flow

The problem described above produces a two-dimensional axisymmetric flow regime. When the outlet radius is much less than the column radius, that is, $r_n\leq R$ , the convergent flow forms near the outlet zone. The solute dispersion coefficient depends on the flow rate, that is, the apparatus-induced dispersion in laboratory apparatus is caused by the nonuniform fluid.

Barry [Reference Barry6] developed an analytical solution for the steady-state case with a constant $\rho _*$ and

$$ \begin{align*} \omega = 0, \quad \frac{\partial ^2\phi ^*}{{\partial z^*}^2}+\frac{1}{r^*}\frac{\partial}{\partial r^*}\bigg( r^*\frac{\partial \phi ^*}{\partial r^*} \bigg) =0,\end{align*} $$

where $\phi $ is the hydraulic head, $q_{r}^{*}=-K(\partial \phi ^*/\partial r^*)$ and $q_{z}^{*}=-K(\partial \phi ^*/\partial z^*)$ . The orientation of the soil column is reversed, as shown in Figure 1, where the flow comes in from the top and exits from the orifice at the bottom. The boundary conditions at $z^*=0$ are [Reference Barry6]

(3.1) $$ \begin{align} \frac{\partial \phi ^*}{\partial z^*}=\bigg( \frac{R^*}{r_{n}^{*}} \bigg) ^2\quad\text{for } 0\le r^*<r_{n}^{*}\quad \text{and} \quad \frac{\partial \phi ^*}{\partial z^*}=0\quad\text{for }r_{n}^{*}\le r^*<1. \end{align} $$

The analytical solutions for the case where fluid exits the column through an orifice with a radius less than the column radius, that is, ${r_n}^*$ = 0.01, 0.25, 0.5 and 0.75, and with no variability in the head condition at the column entry, are compared with the present numerical results. For the comparison, the parameters used in the simulations include $H^*=2$ , $R^*=1$ and $\phi ^* (z^*=H^* )=20$ , which are consistent with those in [Reference Barry6]. The nonuniformities introduced by the orifice radius of the typical soil column in Figure 1(b) are further investigated for the boundary conditions of $p^*=0$ at $z^*=H^*$ and $q^*={q_{in}}^*=1$ at $z^*=0$ . To determine the converging flow zone, both the radial flux $q_r^*$ and the longitudinal flux $q_z^*$ are examined. For an ideal uniform flow in a column, the steady discharge rate is $\mathbf {q}^*=( 0, q_{in}^{*} )$ . The disturbance to $\mathbf {q}^*$ is defined as $\Delta \mathbf {q}^*=( \Delta q_{r}^{*},\Delta q_{z}^{*} )$ , where $\Delta q_{r}^{*}=| q_{r}^{*}-0 |\text { and }\Delta q_{z}^{*}=| q_{z}^{*}-q_{in}^{*} |$ . Following [Reference Barry6], we define the disturbed region of fluid as that in which $| \Delta \mathbf {q}^* |>0.02q_{in}^{*}$ , and the dissipation height of the orifice disturbance is defined as the maximum vertical height $d^*$ between the centre of the orifice and the location where $| \Delta \mathbf {q}^*\text {(}r^*=0) |=0.02q_{in}^{*}$ . The disturbance factor 0.02 was selected to be consistent with that used by Barry [Reference Barry6].

To identify the impact of the orifice on the flow pattern in the column, the flow distributions for ${r_n}^*$ = 0.1, 0.25, 0.5 and 0.75 were simulated and presented in Figure 2. The vectors present the discharge rates and directions, the approximately horizontal black contour lines are hydraulic head and the black lines parallel to discharge rate vectors are streamlines. The red-cross lines represent where $| \Delta \mathbf {q}^* |=0.02q_{in}^{*}$ , and the domain under it is the disturbed zone, that is, $| \Delta \mathbf {q}^* |>0.02q_{in}^{*}$ . It can be seen that the dissipation height $d^*$ is larger for a smaller orifice radius, that is, there is a larger converging flow zone. The nonuniform domain size increases when ${r_n}^*$ decreases. When $r_{n}^{*}\rightarrow 0.1$ , $d^*\rightarrow 3/2$ . On the other hand, when $r_{n}^{*}\rightarrow R$ , no converging flow zone exists, that is, the flow in the column is uniform. These results are in a good agreement with the analytical solution in [Reference Barry6], which provides some supporting evidence that COMSOL Multiphysics is a valid numerical tool for the study of Darcy flows.

Figure 2 Steady flow for $H^*=2$ , $R^*=1$ and various orifice radii: (a) ${r_n}^* = 0.1$ ; (b) ${r_n}^* =0.25$ ; (c) ${r_n}^* = 0.5$ ; and (d) ${r_n}^* = 0.75$ . The vectors present the discharge rates and directions, the approximately horizontal black contour lines are hydraulic head and the black lines parallel to discharge rate vectors are streamlines. The red-cross lines represent where $| \Delta \mathbf {q}^* |=0.02q_{in}^{*}$ , the domain under it is the disturbed zone, that is, $| \Delta \mathbf {q}^* |>0.02q_{in}^{*}$ .

Figure 3 Steady flow for ${r_n}^*=0.1$ , ${R}^*=0.1$ : (a) ${H}^*=4$ ; and (b) ${H}^*=2$ .

The impact of the soil column height on the nonuniform flow was also examined. For the orifice radius ${r_n}^*=0.1$ , the same boundary conditions were used, but column heights of $H^*=2$ and 4 are simulated. Figure 3 shows the flow distributions in the converging zones. It can be seen that the disturbance height $d^*\approx 1.50$ for both cases, and when $z^*>3/2$ in the column, the flow is always uniform. The disturbance distance due to the orifice is summarized in Figure 5. This confirms that the maximum dissipation length scale is less than $3/2$ of the column radius. Therefore, once the column dimension and the orifice size are known, the converging zone height can be determined, which is a useful reference for the design of laboratory transport experiments.

The apparatus-induced dispersion occurs mainly in this converging zone. In other words, outside the convergent flow zone, solute breakthrough curves measured in-situ can be used to obtain accurate soil parameter estimates.

3.2 Solute transport model

The magnitude of the hydrodynamic dispersion results from external forces acting on the liquid and variation in liquid properties, such as density and viscosity. Consider a solute transported in a saturated, uniform flow through a porous medium in a column, as in Figure 1(a) (that is, ${r_n}^*=R^*=1$ ). In [Reference Ogata and Banks19] an analytical solution was developed that described the advection-dispersion transport for uniform flow, that is, the solution of equations (2.7)–(2.12) for a constant $\rho ^*$ and $R^*/H^*\leq 1$ , (or infinite $H^*$ ), as

(3.2) $$ \begin{align} \omega (z,t)=\frac{\omega _{in}}{2}\bigg[ \mathrm{erfc}\bigg( \frac{z^*-q_{z}^{*}t^*}{\sqrt{4D_{z}^{*}t^*}} \bigg) +\exp \!\:\bigg( \frac{q_{z}^{*}z^*}{D_{z}^{*}} \bigg) \mathrm{erfc}\bigg( \frac{z^*+q_{z}^{*}t^*}{\sqrt{4D_{z}^{*}t^*}} \bigg) \bigg], \end{align} $$

where ${D_z}^*$ is the longitudinal dispersion coefficient. Substituting equation (2.15) into the solute dispersion equation (2.8), gives the longitudinal dispersive flux

(3.3) $$ \begin{align} {J_z}^*(t)=-\rho ^*{D_z}^*\frac{\omega _{in}}{\sqrt{4{\pi D_z}^*t^*}} \begin{cases} -\exp \bigg( \dfrac{(z^*-{q_z}^*t^*)^2}{4{D_z}^*t^*} \bigg) -\exp \bigg( \dfrac{z^*{q_z}^*}{{D_z}^*}-\dfrac{(z^*+{q_z}^*t^*)^2}{4{D_z}^*t^*} \bigg)\\[4pt] +\sqrt{\dfrac{4\pi t^*}{{D_z}^*}}{q_z}^*\exp \bigg( \dfrac{z^*{q_z}^*}{{D_z}^*} \bigg) \mathrm{erfc}\bigg( \dfrac{z^*+{q_z}^*t^*}{\sqrt{4{D_z}^*t^*}} \bigg). \end{cases}\nonumber\\[-16pt] \end{align} $$

To validate the COMSOL solute transport model, the breakthrough curve (BTC) for uniform flow and its dispersive flux are calculated and compared with those estimated using the analytical solution of equations (3.2) and (3.3) for the parameters in Table 1.

The analytical and numerical results for the BTC for uniform flow and the dispersive flux at $z^* = 2$ are shown in Figure 4, where reasonably good agreement is evident.

Watson et al. [Reference Watson, Barry, Schotting and Hassanizadeh30] designed laboratory column apparatus for an apparatus-induced dispersion study that consisted of an acrylic column with a detachable base and piston top cap. The data in [Reference Watson, Barry, Schotting and Hassanizadeh30] are also used to verify the COMSOL model for the solute transport in a nonuniform flow. The case CS1-5U in [Reference Watson, Barry, Schotting and Hassanizadeh30] is selected, which has the following parameters in Table 2.

Table 1 Parameters for a uniform flow simulation.

Figure 4 The breakthrough curves, that is, the salt concentration at the exit, for uniform flow are calculated using the analytical and numerical methods.

Table 2 Parameters for Case CS1-5U in [Reference Watson, Barry, Schotting and Hassanizadeh30].

Figure 5 Comparison of the experimental cumulative effluent of solute with the analytical and numerical solutions.

The experimental data for cumulative effluent was compared with the analytical solution using equations (3.2) and (3.3) and the numerical solution (Figure 5). The analytical solution, equation (3.2), derived by Ogata and Banks [Reference Ogata and Banks19] is for a semiinfinite one-dimensional half-space with an imposed uniform flow. It obviously is not able to calculate the nonuniform solute transport accurately, as the additional dispersion caused by the apparatus is not included in these estimates. However, the numerical model provided a relatively accurate estimate of the solute transport, which provides support for the use of COMSOL for studies on apparatus-induced dispersion.

3.3 Apparatus-induced dispersion

When a solute is transported in a nonuniform, density-dependent flow, the displacing solute concentration influences the velocity field. The flow nonuniformity and the density variations will, consequently, affect hydrodynamic dispersion. The transport process is complex. Therefore, the COMSOL model is used to investigate the effects of the apparatus dimensions, the soil parameters and the solution concentration on the dispersion.

3.3.1 Ratio between orifice radius and column radius, ${r_n}^*$

The flow of nonuniformity near the orifice is determined by the ratio of the orifice radius and column radius ${r_n}^*$ . Owing to the convergence, the flow rate in the orifice centre is much greater than that near the sidewall, as discussed in Section 3.1. The variation in local velocity, both in magnitude and direction, causes the solute mass to spread. When the solute front exits, the sidewall lag reduces the average effluent concentration and increases the time over which breakthrough occurs. The dispersive flux of the effluent for several values of ${r_n}^*$ is shown in Figure 6, where $H^* = 2$ for all cases. When ${r_n}^*$ is small, that is, with large flow nonuniformity, the breakthrough occurs earlier and the dispersive flux increases earlier, but the peak value is small. It also takes longer for the dispersive flux to return to zero. This means that the average breakthrough of the effluent takes longer for smaller ${r_n}^*$ . The total dispersive flux is also compared (Figure 7). When ${r_n}^*$ is small, the dispersive flux is smallest, whereas when the flow is uniform the dispersive flux is the largest.

Figure 6 Dispersive flux time series at the orifice for various orifice radii: top figure ${r_n}^*$ ranging from 0.6 to 1.0 and bottom figure ${r_n}^*$ ranging from 0.1 to 0.5.

Figure 7 The relationship between the total dispersive flux and the orifice radii.

3.3.2 Ratio between orifice radius and column height, ${r_n}^*/H*$

As discussed in Section 3.1, the influence of the converging flow can spread upwards to $3/2$ of the column radius from the orifice: that is, once $H^*> 3/2$ , the column height has no influence on the converging flow pattern. However, equation (3.2) has shown that the BTC depends on the solute travelling time for uniform flow. Therefore, at the orifice, the solute concentration will depend on the column height and, consequently, the total dispersive flux is affected. Figure 8(a) shows the dispersive flux at the orifice in relation to the column height where ${r_n}^*$ = 0.14, which is the same as the experimental setting in [Reference Watson, Barry, Schotting and Hassanizadeh30]. It is found that the taller the column, the longer the breakthrough process and the smaller the maximum dispersive flux. However, the total dispersive flux is very close, as shown in Figure 8(b): that is, the column height has little impact on the total dispersive flux when the orifice radius is fixed.

Figure 8 (a) Dispersive flux time series at the orifice for various column heights. (b) The relationship between the dispersive flux and the column height.

Figure 9 Dispersive flux time series at the orifice for various dispersivities.

3.3.3 Longitudinal dispersivity, ${\alpha _L}^*$

The dispersivity is a material property that determines the spreading of solute in a fluid parcel. Hydrodynamic dispersion for the flow in a column, as shown in Figure 1, will be dominated by the longitudinal dispersivity, $\alpha _L$ , in the uniform flow zone according to equation (2.2). Using the parameters in Table 2, simulations for various $\alpha _L$ were carried out to examine its effect on the dispersion. Figure 9 shows that when ${\alpha _L}^*$ increases (0.01275–0.02775), the dispersive flux decreases. However, when ${\alpha _L}^*$ keeps increasing from 0.03225 to 0.05525, the dispersive flux increases and flattens out at the end.

Figure 10 (a) Dispersive flux time series and (b) Vertical flow velocity time series at the orifice for various influent concentrations.

3.3.4 Solute concentration of the influent, $\omega _{in}$

Variations in solute concentration, such as for brine, introduce density and viscosity changes in the fluid. Consequently, this will produce convective currents and affect the flow pattern in the column. The impact of the influent solute concentration, $\omega _{in}$ , is examined using the parameters in Table 1 with $H* = 2$ and ${r_n}^* = 0.14$ : $\omega _{in}$ varies from 0.5% to 20%.

The nondimensional dispersive flux is illustrated in Figure 10(a). It can be seen that the breakthrough occurs earliest and lasts longest for the smallest $\omega_{in}$ . Conversely, for the largest $\omega _{in}$ , the breakthrough curve will be the steepest. The nondimensional dispersive fluxes for various $\omega _{in}$ are listed as in Table 3, where the influent concentration has little impact on the dispersive flux.

Table 3 Normalized cumulative dispersive flux for various $\omega _{in}$ .

When $\omega _{in}$ is high, the density variation in the miscible zone is large. The denser fluid tends to flow against the inflow direction. This is illustrated in Figure 10(b), which shows the longitudinal velocity at the centre of the orifice. From equations (2.5) and (2.6), the density and viscosity depend on the solute concentration. The ratio $\rho /\mu $ can be calculated using equations (2.5) and (2.6), which apparently decreases with the increase in $\omega $ . Consequently, the hydraulic conductivity decreases with the increase in $\omega $ . Hence, the nondimensional velocity for a higher concentration of brine is greater if all other conditions remain the same, as shown in Figure 10(b).

4 Conclusions

In this study, a multiphysics finite-element model was employed to investigate apparatus-induced dispersion in nonuniform, density-dependent flow within a cylindrical column. The accuracy of the model was validated through comparisons with an analytical solution and laboratory column test data. Our findings can be summarized as follows.

  1. (i) The simulations conducted using the model confirmed that flow nonuniformities induced by the apparatus are effectively dissipated within the column when the distance to the outlet exceeds $3R/2$ , where R represents the radius of the column. Additionally, we observed that convergent flow in the column apparatus introduced additional hydrodynamic dispersion, primarily when the height of the soil column was relatively small. However, the influence of convergent flow on dispersion diminished when the height of the column exceeded $3R/2$ .

  2. (ii) The apparatus-induced dispersion characteristics were found to be influenced by soil properties, such as dispersivity. Moreover, we discovered that an increase in the density gradient of the solution resulted in a decrement in flow velocity at the breakthrough point. This phenomenon played a stabilizing role in the flow dynamics and ultimately led to a reduction in dispersive mixing.

  3. (iii) As the dispersion coefficient is known, the dispersive flux increase (or decrease) due to the apparatus, as well as density, can be quantified.

In conclusion, this study provides valuable insights into the phenomenon of apparatus-induced dispersion in nonuniform, density-dependent flow within a cylindrical column. By understanding the dissipation of flow nonuniformities, the effects of convergent flow and the impact of solution density gradient, we can improve our understanding of dispersive mixing phenomena and contribute to the development of more effective strategies for managing and mitigating dispersion in practical applications.

Acknowledgements

This paper is part of a special issue dedicated to honouring the remarkable contributions of Professor Graeme Hocking to the ANZIAM Society, with a specific focus on applied and industrial mathematics.

All the authors of this article have had the distinct privilege of collaborating with Graeme for over 20 years. Hong Zhang embarked on her PhD journey under Graeme’s guidance, studying fluid flow in porous media nearly 30 years ago. Throughout her doctoral studies, Graeme’s invaluable insights and profound knowledge in applied mathematics and fluid dynamics guided and shaped Hong’s research. Following the completion of her PhD, Hong has been incredibly fortunate to continue her association with Graeme, benefiting from his mentorship and enjoying a fruitful professional relationship.

References

Al-Ali, S., Hocking, G. C., Farrow, D. E. and Zhang, H., “A spectral modelling approach for fluid flow into a line sink in a confined aquifer”, European J. Appl. Math. 33 (2022) 960981; doi:10.1017/S0956792521000310.CrossRefGoogle Scholar
Al-Dousari, M. M., Garrouch, A. A and Guedouar, L. H., “Analysis of the convective dispersive transport in porous media”, J. Petroleum Sci. Eng. 66 (2009) 1523; doi:10.1016/j.petrol.2009.01.004.CrossRefGoogle Scholar
Bajracharya, K. and Barry, D. A., “MCMFIT: efficient optimal fitting of a generalized nonlinear advection-dispersion model to experimental data”, Comput. Geosci. 21 (1995) 6176; doi:10.1016/0098-3004(94)00060-8.CrossRefGoogle Scholar
Bajracharya, K., Tran, Y. T. and Barry, D. A., “Cadmium adsorption at different pore water velocities”, Geoderma 73 (1996) 197216; doi:10.1016/0016-7061(96)00048-1.CrossRefGoogle Scholar
Barletta, A., Zanchini, E., Lazzari, S. and Terenzi, A., “Numerical study of heat transfer from an offshore buried pipeline under steady-periodic thermal boundary conditions”, Appl. Therm. Eng. 28 (2008) 11681176; doi:10.1016/j.applthermaleng.2007.08.004.CrossRefGoogle Scholar
Barry, D. A., “Effect of nonuniform boundary conditions on steady flow in saturated homogeneous cylindrical soil columns”, Adv. Water Resources 32 (2009) 522531; doi:10.1016/j.advwatres.2009.01.003.CrossRefGoogle Scholar
Barry, D. A. and Sposito, G., “Analytical solution of a convection-dispersion model with time-dependent transport coefficients”, Water Resources Res. 25 (1989) 24072416; doi:10.1029/WR025i012p02407.CrossRefGoogle Scholar
Bear, J., Dynamics of fluids in porous media (Dover Publication, New York, 1972).Google Scholar
Boon, M., Bijeljic, B., Niu, B. and Krevor, S., “Observations of 3-D transverse dispersion and dilution in natural consolidated rock by X-ray tomography”, Adv. Water Resources 96 (2016) 266281; doi:10.1016/j.advwatres.2016.07.020.CrossRefGoogle Scholar
Chen, J. S., Liu, Y. H., Liang, C. P., Liu, C. W. and Lin, C. W., “Exact analytical solutions for two-dimensional advection-dispersion equation in cylindrical coordinates subject to third-type inlet boundary condition”, Adv. Water Resources 34 (2011) 365374; doi:10.1016/j.advwatres.2010.12.008.CrossRefGoogle Scholar
COMSOL Inc., COMSOL Multiphysics reference manual, version 6.0. www.comsol.com/release/6.0.Google Scholar
Cristea, V. M., Bagiu, E. D. and Agachi, P. S., “Simulation and control of pollutant propagation in Someş river using COMSOL Multiphysics”, Comput. Aided Chem. Eng. 28 (2010) 985990; doi:10.1016/S1570-7946(10)28165-8.CrossRefGoogle Scholar
Danckwerts, P. V., “Continuous flow systems: distribution of residence times”, Chem. Eng. Sci. 2 (1953) 113; doi:10.1016/0009-2509(53)80001-1.CrossRefGoogle Scholar
Hassanizadeh, S. M. and Leijnse, T., “On the modeling of brine transport in porous media”, Water Resources Res. 24 (1988) 321330; doi:10.1029/WR024i003p00321.CrossRefGoogle Scholar
Holm, L. W., “Miscibility and miscible displacement”, J. Petroleum Technol. 38 (1986) 817818; doi:10.2118/15794-PA.CrossRefGoogle Scholar
Kanwar, R. S., Johnson, H. P. and Kirkham, D., “Transport of nitrate and gaseous denitrification in soil columns during leaching”, J. Hydrology 55 (1982) 171184; doi:10.1016/0022-1694(82)90128-7.CrossRefGoogle Scholar
Koch, M. and Zhang, G., “Numerical simulation of the effects of variable density in a contaminant plume”, Groundwater 30 (1992) 731742; doi:10.1111/j.1745-6584.1992.tb01559.x.CrossRefGoogle Scholar
Nick, H. M., Schotting, R., Gutierrez-Neri, M. and Johannsen, K., “Modeling transverse dispersion and variable density flow in porous media”, Transp. Porous Media 78 (2009) 1135; doi:10.1007/s11242-008-9277-x.CrossRefGoogle Scholar
Ogata, A. and Banks, R. B., A solution of the differential equation of longitudinal dispersion in porous media (US Government Printing Office, Washington, 1961); https://pubs.usgs.gov/pp/0411a/report.pdf.CrossRefGoogle Scholar
Overman, A. R., Chu, R. L. and Le, Y., “Kinetic coefficients for phosphorus transport in packed-bed reactor”, J. Water Pollut. Control Federat. 50 (1978) 19051910; http://www.jstor.org/stable/25040371.Google Scholar
Pfannkuch, H. O., “Contribution à l’étude des déplacements de fluides miscible dans un milieu poreux”, Rev. Inst. Français Pétrole 18 (1963) 215270.Google Scholar
Quintard, M., Bertin, H. and Prouvost, L., “Criteria for the stability of miscible displacement through porous columns”, Chem. Ingenieur Technik 59 (1987) 354355; doi:10.1002/cite.330590424.CrossRefGoogle Scholar
Rao, P. S. C., Rolston, D. E., Jessup, R. E. and Davidson, J. M., “Solute transport in aggregated porous media: theoretical and experimental evaluation”, Soil Sci. Soc. Amer. J. 44 (1980) 11391146; doi:10.2136/sssaj1980.03615995004400060003x.CrossRefGoogle Scholar
Sahu, S. and Bhattacharyay, R., “Validation of COMSOL code for analyzing liquid metal magnetohydrodynamic flow”, Fusion Eng. Design 127 (2018) 151159; doi:10.1016/j.fusengdes.2018.01.009.CrossRefGoogle Scholar
Schwartz, R. C., McInnes, K. J., Juo, A. S. R., Wilding, L. P. and Reddell, D. L., “Boundary effects on solute transport in finite soil columns”, Water Resources Res. 35 (1999) 671681; doi:10.1029/1998WR900080.CrossRefGoogle Scholar
Tran, Y. T., Bajracharya, K. and Barry, D. A., “Anomalous cadmium adsorption in flow interruption experiments”, Geoderma 84 (1998) 169184; doi:10.1016/S0016-7061(97)00127-4.CrossRefGoogle Scholar
Tran, Y. T., Barry, D. A. and Bajracharya, K., “Cadmium desorption in sand”, Environ. Internat. 28 (2002) 493502; doi:10.1016/S0160-4120(02)00077-6.CrossRefGoogle ScholarPubMed
Van Genuchten, M. T. and Parker, J. C., “Boundary conditions for displacement experiments through short laboratory soil columns”, Soil Sci. Soc. Amer. J. 48 (1984) 703708; doi:10.2136/sssaj1984.03615995004800040002x.CrossRefGoogle Scholar
Watson, S. J. and Barry, D. A., “Numerical analysis of stable brine displacements for evaluation of density-dependent flow theory”, Phys. Chem. Earth, Part B 26 (2001) 325331; doi:10.1016/S1464-1909(01)00014-4.CrossRefGoogle Scholar
Watson, S. J., Barry, D. A., Schotting, R. J. and Hassanizadeh, S. M., “Validation of classical density-dependent solute transport theory for stable, high-concentration-gradient brine displacements in coarse and medium sands”, Adv. Water Resources 25 (2002) 611635; doi:10.1016/S0309-1708(02)00022-2.CrossRefGoogle Scholar
Wissmeier, L. and Barry, D. A., “Simulation tool for variably saturated flow with comprehensive geochemical reactions in two- and three-dimensional domains”, Environ. Modell. Software 26 (2011) 210218; doi:10.1016/j.envsoft.2010.07.005.CrossRefGoogle Scholar
Wissmeier, L., Barry, D. A. and Phillips, I. R., “Predictive hydrogeochemical modelling of bauxite residue sand in field conditions”, J. Hazardous Mater. 191 (2011) 306324; doi:10.1016/j.jhazmat.2011.04.078.CrossRefGoogle ScholarPubMed
Yong, R. N., Mohamed, A.-M. O. and Warkentin, B. P., Principles of contaminant transport in soils, Volume 73 of Developments in Geotechnical Engineering (Elsevier, Amsterdam, The Netherlands, 1992).Google Scholar
Zhang, H. and Hocking, G. C., “Withdrawal of layered fluid through a line sink in a porous medium”, ANZIAM J. 38 (1996) 240254; doi:10.1017/S0334270000000631.Google Scholar
Zhang, H., Hocking, G. C. and Barry, D. A., “An analytical solution for critical withdrawal of layered fluid through a line sink in a porous medium”, ANZIAM J. 39 (1997) 271279; doi:10.1017/S0334270000008845.Google Scholar
Zhang, H., Hocking, G. C. and Seymour, B., “Critical and supercritical withdrawal from a two-layer fluid through a line sink in a partially bounded aquifer”, Adv. Water Resources 32 (2009) 17031710; doi:10.1016/j.advwatres.2009.09.002.CrossRefGoogle Scholar
Figure 0

Figure 1 Two different soil column study configuration diagrams: (a) experimental set-up of Watson et al. [30]; and (b) set-up of the flow convergence study by Barry [6].

Figure 1

Figure 2 Steady flow for $H^*=2$, $R^*=1$ and various orifice radii: (a) ${r_n}^* = 0.1$; (b) ${r_n}^* =0.25$; (c) ${r_n}^* = 0.5$; and (d) ${r_n}^* = 0.75$. The vectors present the discharge rates and directions, the approximately horizontal black contour lines are hydraulic head and the black lines parallel to discharge rate vectors are streamlines. The red-cross lines represent where $| \Delta \mathbf {q}^* |=0.02q_{in}^{*}$, the domain under it is the disturbed zone, that is, $| \Delta \mathbf {q}^* |>0.02q_{in}^{*}$.

Figure 2

Figure 3 Steady flow for ${r_n}^*=0.1$, ${R}^*=0.1$: (a) ${H}^*=4$; and (b) ${H}^*=2$.

Figure 3

Table 1 Parameters for a uniform flow simulation.

Figure 4

Figure 4 The breakthrough curves, that is, the salt concentration at the exit, for uniform flow are calculated using the analytical and numerical methods.

Figure 5

Table 2 Parameters for Case CS1-5U in [30].

Figure 6

Figure 5 Comparison of the experimental cumulative effluent of solute with the analytical and numerical solutions.

Figure 7

Figure 6 Dispersive flux time series at the orifice for various orifice radii: top figure ${r_n}^*$ ranging from 0.6 to 1.0 and bottom figure ${r_n}^*$ ranging from 0.1 to 0.5.

Figure 8

Figure 7 The relationship between the total dispersive flux and the orifice radii.

Figure 9

Figure 8 (a) Dispersive flux time series at the orifice for various column heights. (b) The relationship between the dispersive flux and the column height.

Figure 10

Figure 9 Dispersive flux time series at the orifice for various dispersivities.

Figure 11

Figure 10 (a) Dispersive flux time series and (b) Vertical flow velocity time series at the orifice for various influent concentrations.

Figure 12

Table 3 Normalized cumulative dispersive flux for various $\omega _{in}$.