SIMPLIFIED DESCRIPTION OF TRANSPORT IN TWO-PHASE ENVIRONMENTS

The concept of retardation is usually applied successfully for the description of transport in two-phase systems of different kind with one moving and one stagnant phase. It simplifies understanding the multitude of processes in the two phases. Retardation factors can be derived from measurements and also mathematically from basic transport equations in the single phases, and thus form a link between practice and theory. The author recently presented a generalization of the theoretical derivations for systems with two moving and/or diffusive phases, in which so called R-factors were introduced. Like retardation factors, R-factors reduce the complexity of the description of a two-phase system. Here conditions are examined in detail under which the generalized R-factor approach is applicable. A demonstration example shows steady states depending on phase dependent reaction rates. A final example application demonstrates how the approach works interpreting real world data.


Introduction
There are very different real world situations, in which transport in 2-phase environments is of relevance.Usually one may think of two liquid phases, like water and oil, or one gaseous and one liquid phase, like oil and gas.However, there are much more other examples.A porous medium is by definition a 2-phase system, which consists of a solid and a fluid part.In most applications the fluid is moving, and thus subject to the processes of advection and diffusion or dispersion.In contrast the solid part can often be conceived as fixed, but there are also important instances, in which the solid phase is also moving or subject to alterations that can be described as diffusion.Sediments are a prominent example of such an environment, which is outlined in more detail below.Here a 2-phase system is understood in the outlined general sense, and a mathematical analytical framework is developed.
Mass transport in 2-phase systems can be generally described by two differential equations, each one including storage, diffusion, advection, reaction and phase exchange for a single phase, where subscripts denote phases 1 and 2.
1 and 2 represent volumetric shares of the phases, D1 and D2 are the diffusivities, V1 and V2 the velocity vectors and λ1 and λ2 the linear production or consumption terms of the phases.For positive λ1 or λ2 there is production, for negative values there are losses.The inter-phase exchange terms are denoted by q12 and q21.The equations describe the mass balance for each phase.
In the sequel we discuss conditions to be fulfilled that the system can be described by only one equation (eqn.3) with additional factors, denoted here by R, Rdiff, Radv and Rreac in the following called R-factors.The term R-factor was first introduced by Holzbecher 1 in a study on the transient change of equilibria in a chain of radionuclides.Holzbecher 2 generalized the approach for multi-process 2-phase systems.Explicit formulae for the R-factors are derived below.As it was derived by adding the concentrations in both phases, equation (3) is the conservation equation for the total amount of the concerned chemical species in the system.
The advantage of eqn.( 3) is that it is much easier to understand than eqns.( 1) and (2).The number of parameters is reduced, the exchange parameter cancelled out.Moreover for the single processes the coefficients can be combined to effective parameters.For example, by comparing eqn.(3)  with the common transport equation, 1D1Rdiff can be understood as an effective diffusivity.Parameter discussions and evaluations become easier that way.
Eqn. (3) as a single equation simplifies the original system consisting of two equations.Thus the wide variety of approaches and solutions for the single transport equation can be utilized.There is a multiplicity of numerical tools.Application example 2 demonstrates the use of a numerical tool.Concerning analytical solutions Ogata and Banks 3 provided a formula for the unsteady 1D advection-diffusion case.Van Genuchten and Alves 4 list much more solutions valid for different boundary conditions, and including reactions.Wexler 5 additionally presents solutions for transport in 2D and 3D.
A plethora of solutions of all kind can be found in the book of Bruggeman. 6New numerical techniques can be employed easily, such as presented by Wang et al. 7 for high Pécletand Courant numbers.
Here conditions are discussed in detail, which have to be fulfilled to make the transition from description (1) and (2) to eqn.(3).At first, the underlying idea is illustrated best by re-calling the concept of retardation factors.Further, the final part shows the applicability of the R-factor approach for a real world situation.

Retardation factors
To illustrate the idea of the simplified description the concept of retardation is briefly re-called here.This concept can be derived rigorously for a 2-phase system, if there is one immobile (static) and one mobile (dynamic) phase.The immobile phase is fixed in the chosen coordinate system.The mobile phase is subject to dispersion and advection.Moreover it is assumed that the inter-phase exchange is governed by equilibrium sorption.In analogy to eqns.( 1) and ( 2) for both phases, separate mass balances can be set up as follows, where c = concentrations,  = volumetric share of the mobile phase, D = dispersion tensor, V = velocity, and qmi = inter-phase exchange term.The subscripts are: m for the mobile phase, i for the immobile phase., D and V in equation ( 4) may vary in space.The term qmi denotes the net flux from the mobile to the immobile phase.
As qmi can hardly be quantified in practice, especially if sorption processes are fast, it is convenient to proceed with a formulation in which the exchange term disappears.This is achieved by adding both equations of system (4).
On the right side of the equation, the concentration of the mobile phase cm remains as the only unknown.The same can be achieved on the left side, if the sorption equilibrium can be described by an isotherm ci(cm).Then the differential equation ( 5) can be rewritten as eqn.(6).(6)   The chain rule is applied for the function ci(cm(t)).In comparison with the transport equation for a tracer component only the term in the brackets on the left side makes a difference.
Thus we may interpret R as parameter that describes the change in the time scale.In case of time independent R one can formally express this by re-writing R∂/∂t=∂/∂(t/R).Because R obviously is greater than 1, there is always a prolongation of the time scale.Thus R is called the retardation factor.Retardation is not a physical process itself, but is a combined result from single phases processes, coupled by phase interaction.
In case of linear isotherms formula (7) can be modified.In porous media studies, it is common to measure fluid phase concentrations cf (=cm) in relation to fluid volume, and solid phase concentrations cs in relation to the rock mass (ci=scs with solid phase density s).Using bulk density ρb=(1-)s the linear isotherm cs=Kdcf leads to the well-established formula (for example Roberts et al., 8 Fitts 9 ): The given derivation follows basic texts on transport in groundwater, for example Kinzelbach 10 or Goode and Konikow. 11Bouwer 12 gives a derivation of the retardation factor without using the chain rule.
The concept of retardation has turned out to be useful in many porous media and groundwater applications.Extensions of the presented basic approach have been derived with respect to upscaling, 13,14 time-dependence, 15 heterogeneity 16 or colloid transport. 17,18For the noted cases these extensions can be performed on top of the generalized approach suggested here, eventually considering some further conditions.

Rfactors for two mobile phases
R-factors can be introduced for two phase systems with two mobile phases in the same way as presented above for R in the mobile-immobile system.Before we come to the general description in this section, it should be noted that additional processes like diffusion and reaction are not taken into account, in order to keep focused on the conditions.Examples are environmental systems with two fluids, like a gas and a liquid, or with two liquids, with advection as dominating intra-phase process.Sediments can be considered as systems with two mobile phases, if sedimentation is taken into account as due to sedimentation, the solid phase has to be described by an advective term, in the Lagrange coordinate system. 19The considered processes are storage, advection and inter-phase exchange.The basic expression of the 2-phase system reads as system (9), where subscripts denote phases 1 and 2, 1 and 2 represent volumetric shares of the phases, and V1 and V2 are velocities of the phases.The inter-phase exchange terms are denoted by q12 and q21.Summing up eqns.( 9) and ( 10) leads to a single differential equation.The exchange terms disappear from the formulation

Condition 2: No temporal changes of volumetric shares (i.e. porosity)
The volumetric shares can be taken out of the time derivative on the left hand side.
The term in the brackets on the left side of the equation is equivalent to the well-known retardation factor R, corresponding to the definition in equation (7).The right hand side has been changed using the product rule in order to prepare the next step.

Condition 4: Steady state flow
Mathematically it can be stated as ∇ × ( 1  1 ) and ∇ × ( 2  2 ) and finally the equation simplifies to the formulation (14), with matrix Radv, containing advection R-factors.In 2D, using velocity components v1x and v1y of V1 and v2x and v2y of V2, the detailed formulation is given by formulation (15).
adv = (  advx 0 0  advy ) with The new factors Radvx, Radvy (and Radvz in 3D) are diagonal components of the Radv matrix appearing in the advection term.The factors are bigger or equal than 1.Radv reduces to the unity matrix in case of V2 = 0, i.e. if the second phase is immobile.The ratio R/Radvx can be interpreted as the factor changing the time scale of the advection process in xdirection.Radvx depends on the ratio between velocities in xdirection of the two phases.
In porous media problems the velocities in the phases are often determined by the hydraulic gradient, but with different hydraulic conductivities.Under this condition the ratios v2x/v1x and v2y/v1y are identical and one ends up with only one advection R-factor, denoted by Radv: For linear exchange, characterized by the constant K=c2/c1, follows: R=1+K1/2 and and Radv=1+K22/V11.With the usual definition of Kd (see above) in the fluid-solid system one obtains: Radv=1+KdbVs/Vf.If both phases move with the same velocity, holds: Radv=R.

Formulation of the general case
A general situation is considered in which both phases may involve the storage, diffusion, advection, reaction and interphase exchange.The basic mathematical description of this situation is given by eqns.( 1) and ( 2).Concerning reactions we consider linear production or consumption of the species.
In the case of a chemical equilibrium (condition 1) one can reduce the system of eqns.( 1) and (2), to a single differential equation.On adding both the equations, the exchange terms vanish, resulting in eqn.(17).
If the phase space does not change temporarily (condition 2), one obtains eqn.(18).
If an isotherm exists (condition 3), one can write If ∇ × ( 1  1 ), ∇ × ( 2  2 ), ∇ × ( 1  1 ), and ∇ × ( 2  2 ) are negligible (condition 4) then which can be re-written as eqn.(22), in line with eqn.(3), using the additional R-factors For the linear isotherm in the fluid-solid system with the diffusivity Ds in the solid phase and Df in the fluid phase the diffusion R-factor can be written as Rdiff=1+KdbDs/Df.When both diffusivities are in the same range, one obtains Rdiff in the same range as the retardation factor R. In case of strong sorption the factor may be several orders of magnitude greater than 1.If velocity ratios in the coordinate directions are the same (v2x/v1x=v2y/v1y=v2x/v1x=v2z/v1z) one can re-write eqn.(22) as eqn.(23).
Equation ( 23) provides the clues how the relevance of a process changes in relation to a single-phase system.There are three process-dependent time scaling factors namely R/Rdiff for diffusion, R/Radv for advection, and R/Rreac for linear reactions.If there is no diffusion in the second phase, Rdiff is equal to 1 and retardation by factor R is valid for diffusion.When there is no advection in the second phase, Radv is equal to 1 and the retardation by R is valid for advection.When there is no production or degradation in the second phase, Rreac is equal to 1 and the retardation factor R is valid concerning reactions.The concept of genuine retardation, as described, is valid only if all factors are equal to unity i.e., Rdiff = Radv = Rreac = 1.The derived single differential equation delivers a tool to understand the relevance of the involved processes for the transport dynamics.

Steady-state phase dependent reactions
Here we give an example for the 1D steady state, if all coefficients are constants.The general solution can be obtained easily with the following ansatz.
Putting this in the differential eqn.(3) one obtains The solution ( 24) and ( 25) can also be expressed in terms of the non-dimensional Péclet and Damköhler numbers.In the presented approach both non-dimensional numbers depend on R-factors: and ( 26) The solution is then given by eqn.(27). 20() =  exp(̅ ) +  1 ̅ exp( ̅ ̅ ) where x = x / L denotes the non-dimensional space variable.
The Péclet and Damköhler numbers determine the steady state solutions, as illustrated for example by Holzbecher. 21e integration constants c1 and As an example a detailed look is presented into an environment, in which one of the phases is subject to degradation and inter-phase exchange only.A typical example for such a system is a porous medium with a fluid and a solid phase, the latter fixed in space.We denote the fluid as phase 1 and the solid as phase 2. As there is no advection and no diffusion in the solid phase, Radv=1, and Rdiff=1 (according to eqns.( 15) and ( 20)).Hence eqn.147 (30) The degradation rate in the mobile phase is λ1 > 0 and in the immobile phase it is λ2 > 0. It is usual that the degradation characteristic is different, i.e., phase dependent.In case of radio-nuclides the degradation rates, i.e., rates for radioactive decay are identical, which is a special case of phase independent reaction.Altogether we can distinguish 5 different cases.In the example in case 2 with 0<λ2<λ1 holds Da = 1.5 with second from top profiles shown in the sub-figures.In the special case of phase independent decay rates Da becomes equal to 2 (cf. the corresponding curves in Figure 1).
In the example case 4 with λ2>λ1 holds Da = 3.The graphs for the latter cases for four different Péclet numbers are the ones with lowest concentrations in the sub-figures.A comparison with case 5 in this series is not possible without further parameter specifications.
The illustrations in each sub-plot of Figure 1 demonstrate the decrease of the concentrations with increasing decay constant, provided all other parameters remain at a constant value.The Damköhler-number is different in the cases defined above, and the change in it is due to the different Rreac factor, on which it depends.

Transient phase-dependent diffusion
As was mentioned above aquatic sediments constitute a physical system in which diffusive and advective processes are relevant in both the fluid and the solid phase.Diffusion is not only caused by processes on molecular scale but also on a larger scale due to the activity of biological species, such as worms, which is referred to as bio-turbation. 22,23orm activity results in re-arrangements in the sediment structure and it is assumed that this can be described as a diffusion process.
Work et al. 24 documented laboratory experiments with aquatic sediments.The experiments were performed in a box with a sediment layer of 23 cm thickness, consisting of medium sand, kaolinite, and topsoil.The overlying free water column was 15 cm deep.The experimental test box was initially filled by a 1 ppt NaCl tracer solution, and then circulated by tracer-free water.Experiments were first done with pure sediments, then worms of type Lumbriculus variegatus were introduced as bio-turbators.Several of such experiment series were performed with different flow speed in the free water column and with different bed materials.
Salinity was measured at several locations of a sampling grid in several depths between 2.5 and 21 cm below the sediment-water interface.Sampling times were selected differently, but typically included 1, 3, 6, 12, and 24 h.Evaluating the different profiles from both situations, without and with worms, the effect of bio-turbation could be studied.For two experiments (B9 without worms and B10 with worms, using the notation of Work et al. 24 ) the measured profiles for NaCl in the sediments after 24 h are depicted by markers in Figure 2. Moreover Work et al. 24 compared the observed concentration profiles with analytical solutions for situations in which only diffusion is a driving process.Reducing the RMS error between measurements and the analytical solution they were able to obtain best-fit values for diffusivities.Their results are listed in Table 1.De denotes the effective diffusivity in case without bio-turbation, Db the effective diffusivity in case of bio-turbation.In all cases the effective diffusivity is increased after the introduction of the worms.The Péclet number, as given in eqn.(26), is a measure for the relative importance of advective in comparison to diffusive processes.If there is diffusion in the solid phase we have to consider a factor Rdiff in the denominator of equation ( 26).The ratio between the Péclet number without bio-turbation, Pe, and with bio-turbation, Peb, is given by eqn.(31).Table 1 shows the resulting values for Rdiff for the different experiments.In correspondence with the derivations here the values are all above 1, ranging between a minimum value of 1.31 and a maximum of 6.3.The experimental results confirm the validity of the theoretical approach.
Based on the mathematical analytical approach given above a 1D numerical model was set up to simulate the concentration distribution in the sediment columns.The initial condition is given by c = 1, representing high salinity in non-dimensional units.At t = 0, we assume a fast fresh water inflow in the surface water and thus set the boundary condition c = 0 at the upper boundary, the sediment-water interface.At the lower boundary a no-flow condition for mass transport is required.The numerical model was set-up using the finite element code COMSOL Multiphysics. 25e model was first run at constant diffusivity for the entire column.For a second run the model region was divided into two parts.The upper sediments are assumed to be affected by worm activity, in our approach implicating an increased diffusivity.The thickness of the upper layer is 6 cm, following Work et al. 24 The value stems from the usual length of the worms, which ranges between 2-5 cm.The value of diffusivity is 10 -4 cm 2 s -1 without bio-turbation and is 410 -4 cm 2 s -1 with bio-turbation, corresponding to a value of 4 for Rdiff (eqn.31).In this way two related experiments were modelled, one without and one with diffusion in the solid phase, i.e., without and with bio-turbation.
Figure 2 depicts the results of the numerical modelling for experiments B9 (without bio-turbation) and B10 with bioturbation.The comparison with the measurements shows a rough agreement.For the B9 experiment the uppermost observation is highly overestimated, while others are underestimated.Thus the slope within the profile is not matched, while the mean square error is minimal.Probably one reason for the underestimation is that the c = 0 boundary condition is not matched in reality as it neglects the time to wash out high salinity at the bottom of the surface water.The B10 experiment with bio-turbation seems to match better.Concentrations in the uppermost layer are underestimated, while there is an overestimation by the numerical model in the lower part.The latter may not be taken as serious, because the initialisation of the physical model with equally high NaCl concentrations may not have fully completed, as mentioned by the experimenters, 24 and as seen in the fluctuations of the measured values.Thus it seems to be more probable that the slope of the concentration profile is shifted to lower concentration values by the numerical model.
In order to obtain a better fit in the mentioned sense, a third model run was performed in which the bio-turbation model was varied from a high value at the top to a low background value at the bottom.This takes into account that the worms are of different size and they may not drill vertically to reach deeper levels.The re-arrangements in the sediment structure can thus be expected to be highest near to the sediment-water interface and decrease gradually with distance from it.Thus the diffusivity is assumed to vary linearly from a high value at the upper boundary to the bottom of the upper sediment layer.Correspondingly Rdiff changes from a value of 4 at the sediment-water interface to 1 to the interior sediment interface.
Figure 2 also depicts the simulation results of the third model run.Obviously these match better than the former results in the sense noted above.The deviations between measurements and model at the end of the column can be attributed to difficulties of the experimental set-up.However the slope of the concentration profile is predicted exactly by the model.

Conclusion
The presented approach reduces the complexity of transport processes in general two-phase multi-process environments.Assuming the validity of chemical equilibrium the description can be reduced to a single differential equation, in which coefficients are modified by R-factors.For the involved processes R-factors were introduced in analogy to the derivation of the well-known retardation factor for the storage process.149 The reduction to a single equation allows the utilization of classical approaches and results concerning the advectiondiffusion-reaction equation.As example we present analytical solutions for the 1D steady state.In contrast to the retardation factor the R-factors have an effect not only on the transient development but also on the steady state.The effect can easiest be evaluated by the use of extended definition of Pécletand Damköhler numbers.
In an application concerning aquatic sediments the theoretical results are contrasted with the outcome from laboratory experiments.The comparison shows that the theory gives a correct explanation for the behaviour observed in the sediments.It also shows how a quantification of the effect of bio-turbation is given by the Rdiff factor.The utilization of the R-factor approach in a numerical model with spatially varying bio-turbation delivered a better match with measurements.

Figure 2 .
Figure 2. Concentration profiles of NaCl in aquatic sediments with and without worm activity; normalized concentration is plotted versus depth; markers indicate measured values and curves visualize numerical solutions.

Table 1 .
Changes of effective diffusivities due to bio-turbation in the experiments of Work et al.