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.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).