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

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.


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 [2,15,17,18,22,26,29,[34][35][36].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, [9,33]).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 [3].This experimental design allows the use of one-dimensional or two-dimensional solute transport models [4,6,10,16,20,23,26,27,32], 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 [8,21].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 [6], 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 [6] 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 [1,5,12,24,31].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.

Governing equations
We consider a typical experiment carried out in cylindrical soil columns [30] 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 κ, 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 ρ, dynamic viscosity μ and salt mass fraction ω.
The flow and the solute satisfy Darcy's law and the mass conservation equations, which are summarized as follows.(See [30].) Darcy's law and Fick's law Here, q denotes the specific discharge, p is the fluid pressure, g = (0, 0, g) is the gravitational acceleration, n is the porosity, h is the dispersive mass flux and D = (D ij ) the hydrodynamic dispersion tensor [8], where D m is the mass diffusivity and α L and α T are the longitudinal and transversal dispersivity, respectively.
Conservation of mass of the water and solute where ρ is the fluid density, μ is the dynamic viscosity of fluid and ω is the fractional solute concentration.The saltwater density and viscosity, respectively, are given by, In the equations of state (2.5) and (2.6) for brine, the constant γ = 0.6923 [14], while, for the freshwater reference state, ρ 0 = 998.23 kg/m 3 and μ 0 = 0.001 N S /m 2 at 20 • C and 1 atm.
The full set of equations is coupled between density ρ, viscosity μ, specific discharge q and dispersive flux J.All variables with dimensions of length are scaled with R, the radius of the column, so that The other nondimensional variables are and g * = g g .
The reference hydraulic conductivity, K = κρ 0 g/μ 0 with κ as the permeability, is used to scale time t * = tK/nR, the discharge flux q * = q/q in , the hydrodynamic dispersion tensor D * = D/Rq in and the dispersive mass flux J * = J/ρ 0 Rq in , where q in is the discharge flux for 0 < r < R and z = 0. Then equations (2.1)-(2.6)can be rewritten as follows.
Darcy's law and Fick's law ) Conservation of mass with ρ * = exp(γω), ( 2 . 1 1 ) The solute transport depends on the flow, as the velocity q * determines both advective and dispersive transport.The density ρ * (ω) and viscosity μ * (ω) 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 ω and p * as functions of z * , r * and t * .

Initial and boundary conditions
Following the experiment set-up of [29], 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 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, [7,28]) was used, where the measured flow rates q * in and solute concentration ω in are specified: that is, At the orifice, the most realistic assumption for the finite column apparatus is the Danckwerts boundary condition [13,25,28] for the solute concentration and constant pressure for the fluid: that is, 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.

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 ≤ 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 [6] developed an analytical solution for the steady-state case with a constant ρ * and where φ is the hydraulic head, q * r = −K(∂φ * /∂r * ) and q * z = −K(∂φ * /∂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 [6] 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 φ * (z * = H * ) = 20, which are consistent with those in [6].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 q * = (0, q * in ).The disturbance to q * is defined as Following [6], we define the disturbed region of fluid as that in which |Δ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 |Δq * (r * = 0)| = 0.02q * in .The disturbance factor 0.02 was selected to be consistent with that used by Barry [6].
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 |Δq * | = 0.02q * in , and the domain under it is the disturbed zone, that is, |Δ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 → 0.1, d * → 3/2.On the other hand, when r * n → 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 [6], which provides some supporting evidence that COMSOL Multiphysics is a valid numerical tool for the study of Darcy flows.
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 * ≈ 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.

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 [19] 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 ρ * and R * /H * ≤ 1, (or infinite H * ), as where D z * is the longitudinal dispersion coefficient.Substituting equation (2.15) into the solute dispersion equation (2.8), gives the longitudinal dispersive flux 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. [30] 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 [30] are also used to verify the COMSOL model for the solute transport in a nonuniform flow.The case CS1-5U in [30] is selected, which has the following parameters in Table 2.
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 [19] 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 [9]

Solute concentration
Time 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 [30].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.

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.

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.

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 [30].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.

Longitudinal dispersivity, α 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, α L , in the uniform flow zone according to equation (2.2).Using the parameters in Table 2, simulations for various α L were carried out to examine its effect on the dispersion.Figure 9 shows that when α L * increases (0.01275-0.02775),    the dispersive flux decreases.However, when α L * keeps increasing from 0.03225 to 0.05525, the dispersive flux increases and flattens out at the end.

Solute concentration of the influent, ω 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, ω in , is examined using the parameters in Table 1 with H * = 2 and r n * = 0.14: ω 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 ω in .Conversely, for the largest ω in , the breakthrough curve will be the steepest.The nondimensional Cumulative dispersive flux 0.0088 0.0086 0.0082 0.0080 0.0078 dispersive fluxes for various ω in are listed as in Table 3, where the influent concentration has little impact on the dispersive flux.
When ω 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 ρ/μ can be calculated using equations (2.5) and (2.6), which apparently decreases with the increase in ω.Consequently, the hydraulic conductivity decreases with the increase in ω.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).

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

FIGURE 1 .
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 2 .FIGURE 3 .
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 |Δq * | = 0.02q * in , the domain under it is the disturbed zone, that is, |Δq * | > 0.02q * in .

4 FIGURE 6 . 3 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.

4 FIGURE 8 .
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 .
FIGURE 9. Dispersive flux time series at the orifice for various dispersivities.

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

TABLE 1 .
Parameters for a uniform flow simulation.

TABLE 3 .
Normalized cumulative dispersive flux for various ω in .