**Exploring Compatibility of Sherwood-Gilland NAPL Dissolution Models with Micro-Scale Physics Using an Alternative Volume Averaging Approach**

#### **Scott K. Hansen**

Zuckerberg Institute for Water Research, Ben-Gurion University of the Negev, Beer-Sheva 84105, Israel; skh@bgu.ac.il

Received: 21 May 2019; Accepted: 16 July 2019; Published: 23 July 2019

**Abstract:** The dynamics of NAPL dissolution into saturated porous media are typically modeled by the inclusion of a reaction term in the advection-dispersion-reaction equation (ADRE) with the reaction rate defined by a Sherwood-Gilland empirical model. This stipulates, among other things, that the dissolution rate is proportional to a power of the NAPL volume fraction, and also to the difference between the local average aqueous concentration of the NAPL species and its thermodynamic saturation concentration. Solute source models of these sorts are ad hoc and empirically calibrated but have come to see widespread use in contaminant hydrogeology. In parallel, a number of authors have employed the method of volume averaging to derive upscaled transport equations describing the same dissolution and transport phenomena. However, these solutions typically yield forms of equations that are seemingly incompatible with Sherwood-Gilland source models. In this paper, we revisit the compatibility of the two approaches using a radically simplified alternative volume averaging analysis. We begin from a classic micro-scale formulation of the NAPL dissolution problem but develop some new simplification approaches (including a physics-preserving transformation of the domain and a new geometric lemma) which allow us to avoid solving traditional closure boundary value problems. We arrive at a general, volume-averaged governing equation that does not reduce to the ADRE with a Sherwood-Gilland source but find that the two approaches do align under straightforward advection-dominated conditions.

**Keywords:** NAPL; volume averaging; upscaling; mass transfer

#### **1. Introduction**

We consider the modeling of dissolution of residual non-aqueous phase liquids (NAPL)—that is, liquids that are essentially immiscible in water and which are trapped by capillary forces in porous media—into pore water and the subsequent transport of solute originating in the NAPL. Transport can be modeled without difficulty by means of the advection-dispersion-reaction equation (ADRE):

$$\frac{\partial \mathcal{C}}{\partial t} + \boldsymbol{\sigma} \cdot \nabla \mathcal{C} - \nabla \cdot \mathbf{D}\_{\mathrm{eff}} \nabla \mathcal{c} = \mathbf{Q}\_{\prime} \tag{1}$$

where *v* is the pore water velocity, *D*eff is an effective Fickian dispersion coefficient capturing the effects of local-scale hydrodynamic dispersion and diffusion, and *Q* is a source term. (See the nomenclature in Appendix A for units and definitions of all symbols.) However, it is necessary to have a model of *Q* that captures the physics of NAPL dissolution. To this end, a number of similar empirical models of *Q* have been proposed [1–6] which are of the form *Q* = *K* (*c*sat − *c*), where *c*sat is the thermodynamic saturation aqueous concentration of the NAPL chemical species and *K* is a mass transfer coefficient. These models, known as Sherwood-Gilland models [7], most generally express *K* as

$$K = \kappa \frac{\partial}{\partial\_m^2} \eta\_n^{\beta} \text{Re}^{\gamma},\tag{2}$$

where *α*, *β* and *γ* are fitting parameters, D is the Fickian diffusivity, *dm* is a characteristic length of the porous media, *η<sup>n</sup>* is the bulk volume fraction of NAPL, and Re is the Reynolds number. Certain empirical studies combine terms and do not calibrate all of *α*, *β* and *γ*, but with the exception of [6], all contain terms *η<sup>β</sup> <sup>n</sup>* , or equivalently *<sup>S</sup><sup>β</sup> <sup>n</sup>* (with NAPL saturation expressed as a fraction of pore volume rather than bulk volume), with calibrated values of 0.6 to 1.24 for *β* ([8], see Table 2) The theory underlying the form of such models is out of scope for this paper but is discussed in more detail in References [1,8]. The importance of the Sherwood-Gilland approach lies in the relation of *Q* to local concentration in the water phase and a small number of measurable system parameters.

By contrast with the Sherwood-Gilland empirical approach, the upscaled dynamics of residual NAPL dissolution have also been considered using volume averaging theory [9]. While there is not room for a full account here, volume averaging operates by defining a superficial volume average operator

$$
\langle \cdot \rangle = \frac{1}{\mathcal{V}} \int\_{V(\mathbf{x})} \cdot dV \,\tag{3}
$$

where *V* is a simply-connected region of volume V centered at *x*. From the superficial average, phase-intrinsic volume average operators are defined for each phase *i* according to · = *η<sup>i</sup>* · *i* , where *i* stands in for any of the three phases in the system: water (*w*), NAPL (*n*), and solid (*s*) and *ηi*(*x*) is the volume fraction of *V*(*x*) occupied by the *i* phase. Inside *V*(*x*), the surface between phases *i* and *j* is denoted by Γ*ij* and an infinitesimal element of that surface is denoted by *dAwn*. The unit vector normal to Γ*ij*, pointing into phase *j*, is denoted by *nij*. Although all averaged quantities are functions of location, *x*, and time, *t*, these parameters are generally suppressed. Employing these concepts, we may state a volume averaging theorem [9–11] for *c* (and analogously for other scalar quantities),

$$
\langle \nabla c \rangle = \nabla \left\langle c \right\rangle + \sum\_{j \neq w} \frac{1}{j'} \int\_{\Gamma\_{wj}} c \mathfrak{n}\_{wj} dA\_{wj'} \tag{4}
$$

and for *v* (and analogously for other vector-valued quantities),

$$
\langle \nabla \cdot \mathbf{v} \rangle = \nabla \cdot \langle \mathbf{v} \rangle + \sum\_{j \neq w} \frac{1}{\mathcal{V}} \int\_{\Gamma\_{wj}} \mathbf{v} \cdot \mathbf{n}\_{wj} dA\_{wj}.\tag{5}
$$

The essence of the volume averaging technique is that spatially-nonuniform dependent variables are partitioned into volume averages and small-scale fluctuations and the governing PDE for the un-averaged dynamics is rewritten in terms of both the volume average quantities and spatial integrals of the fluctuations using (4) and (5). Finally, an equation in terms of the volume averaged terms alone is developed by application of physical and mathematical knowledge to rewrite the fluctuation integrals in terms of the volume averaged quantities [12]. Commonly, the closure is accomplished by specifying boundary value problems (BVPs) for the fluctuations in which the volume average quantities participate as forcing functions, as in Reference [9]. If one is only interested in the form of the volume averaged PDEs, these IBVPs do not need to be solved explicitly. If quantitative results are required, then analytical [13] or numerical [14,15] solution of the IBVP in a characteristic cell may be indicated, using for example, Lattice-Boltzmann methods [16].

The volume averaging approach has been applied to the dissolution of residual NAPL, most notably in a pair of papers by Quintard and Whitaker [13,17], each employing a different micro-scale formulation but both arriving at the same upscaled (i.e., volume averaged) equation, which is apparently incompatible with (1):

$$\begin{split} \eta\_{w} \frac{\partial \left< \boldsymbol{c} \right>^{w}}{\partial t} + \left[ \eta\_{w} \left< \boldsymbol{v} \right>^{w} - \boldsymbol{u} \right] \cdot \nabla \left< \boldsymbol{c} \right>^{w} - \nabla \cdot \left( \mathbf{D}\_{w}^{\bullet} \cdot \nabla \left< \boldsymbol{c} \right>^{w} \right) + \nabla \eta\_{w} \cdot \left< \boldsymbol{\mathcal{O}} \nabla \left< \boldsymbol{c} \right>^{w} \\ &= \nabla \cdot \left[ \mathbf{d}\_{\text{ib}} \left( \left< \boldsymbol{c} \right>^{w} - \mathbf{c}\_{\text{sat}} \right) \right] - \mathbf{x} \left( \left< \boldsymbol{c} \right>^{w} - \mathbf{c}\_{\text{sat}} \right) . \end{split} \tag{6}$$

Equation (6) contains two quantities, *D∗ <sup>w</sup>* and *dw*, that are functions of the formal solutions to two distinct closure problems, along with *u* and *κ* which are forcing quantities in those respective closure problems and which are themselves formally defined in terms of micro-scale surface integrals that must be expressed in terms of larger-scale quantities to fully close the problem. See Reference [13] for full details on the underlying assumptions and the volume averaging procedure. An identical treatment of the water-NAPL mass transfer equations appears in the four-phase (also including biofilm) analysis of Reference [18], with a minor simplification appearing in Reference [19], which embeds an apparent assumption that *dw* is spatially-uniform. An alternative analysis which makes different closure assumptions [20] also arrives at an equation with two sources that seemingly cannot be re-expressed in Sherwood-Gilland form.

In Reference [17], the authors show that an ADRE with Sherwood-Gilland source term can be derived by imposing four scale restrictions that they dub "dominant convective transport." However, some of these scale restrictions involve the quantities *D∗ <sup>w</sup>* and *dw*, which do not have a straightforward interpretation, and so it is difficult to know when or if they are satisfied. In a similar study on biofilm growth (where the "NAPL" grows according to Monod kinetics, rather than shrinks by dissolution) [21], the authors obtained a closed, classic ADRE and in a follow-up study [22], the authors showed that either ADRE or non-classical volume average behavior occurred, depending on the coupling assumptions employed.

Informed by these works, in this note we reconsider the compatibility of volume averaged micro-scale physics and the ADRE with a Sherwood-Gilland source. We begin from the microscopic formulation employed by Reference [13]. However, we perform different simplifying analysis, obviating the need for formal solution of closure boundary value problems relating fluctuations to volume averaged quantities, yielding a slightly different volume averaged governing equation in terms of independent variables *η<sup>n</sup>* and *c <sup>w</sup>*, still ostensibly with multiple sources, which is incompatible with (1). We then show how using any Sherwood-Gilland model as a second governing equation relating *η<sup>n</sup>* and *c <sup>w</sup>* causes one of the source terms to be negligible under advection dominated conditions much more straightforward than those in [17], allowing emergence of the expected form of upscaled ADRE:

$$\frac{\partial \left< \boldsymbol{c} \right>^{w}}{\partial t} + \left< \boldsymbol{\nu} \right>^{w} \cdot \nabla \left< \boldsymbol{c} \right>^{w} - \nabla \cdot \mathbf{D}\_{\text{eff}} \nabla \left< \boldsymbol{c} \right>^{w} \approx \eta\_{u}^{\beta} \left( \mathbf{c}\_{\text{sat}} - \left< \boldsymbol{c} \right>^{w} \right). \tag{7}$$

#### **2. Formulation**

#### *2.1. Governing Equations*

We are attempting to model a three-phase system, a characteristic portion of which is illustrated in Figure 1, each phase of which occupies a respective volume fraction of averaging domain *V*(*x*) *ηw*(**x**, *t*), *ηn*(**x**, *t*), and *ηs*, which satisfy

$$
\eta\_w + \eta\_n + \eta\_5 = 1.\tag{8}
$$

We assume that the Darcy flux, *q*, for the system as a whole is fixed and satisfies ∇ · *q* = 0. Then, it follows from the chain rule and the fact that *v <sup>w</sup>* = *q*/*ηw*, that

$$
\nabla \cdot \langle \boldsymbol{v} \rangle^{\boldsymbol{w}} = -\frac{\boldsymbol{q}}{\eta\_{\boldsymbol{w}}^{2}} \cdot \nabla \eta\_{\boldsymbol{w}}.\tag{9}
$$

The ADE applies in the water phase:

$$\frac{\partial \mathcal{C}}{\partial t} + \nabla \cdot (\mathcal{c}\boldsymbol{\sigma}) = \nabla \cdot \left(\mathcal{O}\nabla \boldsymbol{c}\right). \tag{10}$$

The solid and NAPL phases are homogeneous and do not need internal mass balance equations. However, the NAPL phase shrinks over time, yielding space to the water phase. This mass balance can be expressed by:

$$
\rho\_n \frac{\partial \eta\_n}{\partial t} = \frac{1}{\mathcal{Y}} \int\_{\Gamma\_{\text{un}}} \mathcal{O} \nabla \mathbf{c} \cdot \mathbf{n}\_{\text{un}} dA\_{\text{un}}.\tag{11}
$$

We work extensively with this term and so make the definition

$$Q \equiv -\rho\_n \frac{\partial \eta\_n}{\partial t},\tag{12}$$

where *Q* represents the upscaled source per unit time from the NAPL phase into the water phase.

**Figure 1.** Schematic diagram of blobs of residual NAPL in pores defined by solid media otherwise saturated with water, with averaging volume superimposed.

#### *2.2. Coupling Conditions*

We assume no-slip boundary conditions, local equilibrium at the surface of the NAPL and no flux into the solid. Mathematically, these boundary conditions are expressed:

$$
\sigma = 0 \text{ к} \epsilon \Gamma\_{\text{uv}\prime} \tag{13}
$$

$$
\sigma = 0 \text{ ж}\epsilon \Gamma\_{w \bowtie \prime} \tag{14}
$$

$$c = c\_{\text{sat}} \ge \varepsilon \Gamma\_{\text{var}\_f} \tag{15}$$

$$
\nabla \boldsymbol{\mathcal{L}} \cdot \boldsymbol{\mathfrak{n}}\_{\text{uv}} = \boldsymbol{0} \,\, \mathbf{x} \boldsymbol{\varepsilon} \Gamma\_{\text{uv}} \,. \tag{16}
$$

The solid volume fraction does not change with time, so the volume fractions can be coupled according to the following relationship:

$$\frac{\partial \eta\_w}{\partial t} = -\frac{\partial \eta\_n}{\partial t}.\tag{17}$$

#### **3. Simplified Volume Averaging Analysis**

We begin by applying the superficial averaging operator · to (10), yielding

$$\frac{\partial \langle \mathbf{c} \rangle}{\partial t} + \langle \nabla \cdot (\mathbf{c} \mathbf{v}) \rangle = \langle \nabla \cdot (\mathcal{O} \nabla \mathbf{c}) \rangle \,. \tag{18}$$

We manipulate each of the three terms separately to write them all in terms of intrinsic volume averages, simplifying each. We then assemble them and simplify further.

#### *3.1. Accumulation Term*

There is next to no analysis required here. The accumulation term can be expressed in terms of intrinsic volume averages directly

$$\frac{\partial \langle \mathbf{c} \rangle}{\partial t} = \eta\_w \frac{\partial \langle \mathbf{c} \rangle^w}{\partial t} + \frac{\partial \eta\_w}{\partial t} \left\langle \mathbf{c} \right\rangle^w. \tag{19}$$

*3.2. Advection Term*

From volume averaging theorem (5):

$$
\langle \nabla \cdot (c\mathbf{v}) \rangle = \nabla \cdot \langle (c\mathbf{v}) \rangle + \frac{1}{\mathcal{V}} \int\_{\Gamma\_{\text{uv}}} c\mathbf{v} \cdot \mathfrak{n}\_{\text{uv}} dA\_{\text{uv}} + \frac{1}{\mathcal{V}} \int\_{\Gamma\_{\text{uv}}} c\mathbf{v} \cdot \mathfrak{n}\_{\text{uv}} dA\_{\text{uv}}.\tag{20}
$$

from which it follows from no-slip conditions (13) and (14), that ∇ · (*cv*) = ∇ · *cv* = ∇ · *η<sup>w</sup> cv w*. Further simplification follows from writing *c* and *v* in terms of their volume average and a perturbation: *c* = *c <sup>w</sup>* <sup>+</sup> *<sup>c</sup>*˜, and *<sup>v</sup>* <sup>=</sup> *<sup>v</sup> <sup>w</sup>* <sup>+</sup> *<sup>v</sup>*˜. Then, we may write *<sup>c</sup><sup>v</sup> <sup>w</sup>* <sup>=</sup> *<sup>c</sup> <sup>w</sup> <sup>v</sup> <sup>w</sup>* <sup>+</sup> *<sup>c</sup>*˜*v*˜ *<sup>w</sup>*. We have

$$\begin{split} \left< \nabla \cdot \left( \varepsilon \boldsymbol{\nu} \right) \right> &= \nabla \cdot \left( \eta\_{\boldsymbol{\nu}} \left< \boldsymbol{c} \right>^{\boldsymbol{\mu}} \left< \boldsymbol{\nu} \right>^{\boldsymbol{\mu}} \right) + \nabla \cdot \left( \eta\_{\boldsymbol{\nu}} \left< \boldsymbol{c} \boldsymbol{\mathcal{O}} \right>^{\boldsymbol{\mu}} \right) \\ &= \nabla \eta\_{\boldsymbol{\nu}} \cdot \left< \boldsymbol{c} \right>^{\boldsymbol{\mu}} \left< \boldsymbol{\nu} \right>^{\boldsymbol{\mu}} + \eta\_{\boldsymbol{\nu}} \nabla \left< \boldsymbol{c} \right>^{\boldsymbol{\mu}} \cdot \left< \boldsymbol{\nu} \right>^{\boldsymbol{\mu}} + \eta\_{\boldsymbol{\nu}} \left< \boldsymbol{c} \right>^{\boldsymbol{\mu}} \nabla \cdot \left< \boldsymbol{\nu} \right>^{\boldsymbol{\mu}} + \eta\_{\boldsymbol{\nu}} \nabla \cdot \left< \boldsymbol{c} \boldsymbol{\mathcal{O}} \right>^{\boldsymbol{\mu}} + \nabla \eta\_{\boldsymbol{\nu}} \cdot \left< \boldsymbol{\mathcal{O}} \boldsymbol{\mathcal{O}} \right>^{\boldsymbol{\mu}} \end{split} \tag{21}$$

Note that we can manipulate this by applying (9), so that

$$\eta\_{\rm w} \left< \mathfrak{c} \right>^{\rm w} \nabla \cdot \left< \mathfrak{v} \right>^{\rm w} = \eta\_{\rm w} \left< \mathfrak{c} \right>^{\rm w} \left( -\frac{\mathfrak{q}}{\eta\_{\rm w}^2} \cdot \nabla \eta\_{\rm w} \right) = -\nabla \eta\_{\rm w} \left< \mathfrak{c} \right>^{\rm w} \cdot \left< \mathfrak{v} \right>^{\rm w},\tag{22}$$

which can be back-substituted into (21), canceling and leading to the conclusion

$$
\langle \nabla \cdot (\varepsilon \mathbf{v}) \rangle = \eta\_w \nabla \left< \varepsilon \right>^w \cdot \left< \mathbf{v} \right>^w + \eta\_w \nabla \cdot \left< \varepsilon \overline{\mathbf{v}} \right>^w + \nabla \eta\_w \cdot \left< \varepsilon \overline{\mathbf{v}} \right>^w. \tag{23}
$$

For our purposes, we are primarily interested in the form of the source terms. It may thus be reasonable to employ a mean field approximation and discard the small double fluctuation term *c*˜*v*˜, however we retain it presently, later to be incorporated into an effective dispersion.

#### *3.3. Diffusion Term*

Again applying volume averaging theorem (5), we see that

$$
\langle \nabla \cdot (\mathcal{O} \nabla \mathfrak{c}) \rangle = \nabla \cdot \langle \mathcal{O} \nabla \mathfrak{c} \rangle + \frac{\mathcal{\mathcal{O}}}{\mathcal{Y}} \int\_{\Gamma\_{\text{uv}}} \nabla \mathfrak{c} \cdot \mathfrak{n}\_{\text{uv}} dA\_{\text{uv}} + \frac{\mathcal{\mathcal{O}}}{\mathcal{Y}} \int\_{\Gamma\_{\text{uv}}} \nabla \mathfrak{c} \cdot \mathfrak{n}\_{\text{uv}} dA\_{\text{uv}} \tag{24}$$

or, by application of the no-flux boundary conditions, that

$$
\langle \nabla \cdot (\mathcal{O} \nabla \boldsymbol{c}) \rangle = \nabla \cdot \mathcal{O} \left\langle \nabla \boldsymbol{c} \right\rangle + \frac{\mathcal{\boldsymbol{\mathcal{O}}}}{\mathcal{Y}} \int\_{\Gamma\_{\text{uv}}} \nabla \boldsymbol{c} \cdot \mathfrak{n}\_{\text{uv}} dA\_{\text{uv}\prime} \tag{25}
$$

where we observe by direct inspection that the second term is a volumetric NAPL sink, representing diffusion of solute from the water phase *into* the NAPL. We thus replace it with the symbol −*Q*, as defined in (12), emphasizing this role. By another application of an averaging theorem (4), it follows that

$$
\nabla \cdot \mathcal{O} \left< \nabla \boldsymbol{\varepsilon} \right> = \mathcal{O} \nabla \cdot \nabla \left< \boldsymbol{\varepsilon} \right> + \mathcal{O} \nabla \cdot \left[ \frac{1}{\mathcal{V}} \int\_{\Gamma\_{\mathrm{ws}}} \mathrm{c} \mathbf{u}\_{\mathrm{ws}} dA\_{\mathrm{ws}} \right] + \mathcal{O} \nabla \cdot \left[ \frac{1}{\mathcal{V}} \int\_{\Gamma\_{\mathrm{wn}}} \mathrm{c} \mathbf{u}\_{\mathrm{wn}} dA\_{\mathrm{wn}} \right], \tag{26}
$$

We now strongly diverge from the approach in Reference [13]. We eliminate the middle term on the RHS by noting that the average of *nws* at the elements of Γ*ws* intersecting any plane is **0**, and that *V* can be partitioned into planar slices, each orthogonal to ∇ *c*, on which the expected value of *c* is constant. Assuming ergodicity, we conclude that the portion in square brackets is **0**, independent of location, meaning that its divergence is also zero. To simplify the third term, we apply a result developed in Appendix B ( <sup>1</sup> V <sup>Γ</sup>*wn <sup>n</sup>wndAwn* = −∇*ηw*). Observing that *<sup>c</sup>* = *<sup>c</sup>*sat everywhere on <sup>Γ</sup>*wn*, we see that:

$$
\nabla \cdot \mathcal{O} \left< \nabla \mathbf{c} \right> = \eta\_w \mathcal{O} \nabla^2 \left< \mathbf{c} \right>^w + 2 \mathcal{O} \nabla \eta\_w \cdot \nabla \left< \mathbf{c} \right>^w + \mathcal{O} \nabla^2 \eta\_w \left[ \left< \mathbf{c} \right>^w - \mathbf{c}\_{\text{sat}} \right].\tag{27}
$$

Thus,

$$
\langle \nabla \cdot (\mathcal{O} \nabla \mathfrak{c}) \rangle = \eta\_w \mathcal{O} \nabla^2 \langle \mathfrak{c} \rangle^w + 2 \mathcal{O} \nabla \eta\_w \cdot \nabla \left\langle \mathfrak{c} \right\rangle^w + \mathcal{O} \nabla^2 \eta\_w \left[ \langle \mathfrak{c} \rangle^w - \mathfrak{c}\_{\mathsf{sat}} \right] - Q. \tag{28}
$$

#### *3.4. Assembly and Simplification of Governing Equation*

Inserting (19), (23), and (28) into (18), and placing the variables in each of the terms in a consistent order yields:

$$\eta\_{\overline{w}} \frac{\partial \langle c \rangle^{w}}{\partial t} + \frac{\partial \eta\_{\overline{w}}}{\partial t} \left\langle c \right\rangle^{w} + \eta\_{\overline{w}} \left\langle \mathbf{v} \right\rangle^{w} \cdot \nabla \left\langle c \right\rangle^{w} + \eta\_{\overline{w}} \nabla \cdot \left\langle \vec{c} \vec{\sigma} \right\rangle^{w} + \nabla \eta\_{\overline{w}} \cdot \left\langle \vec{c} \vec{\sigma} \right\rangle^{w} = -\mathcal{J} \nabla^{2} \eta\_{w} c\_{\text{sat}} - Q. \tag{29}$$

We here diverge again from the customary volume averaging approach by immediately making the closure approximation *c*˜*v*˜ *<sup>w</sup>* <sup>=</sup> <sup>−</sup>*k*<sup>∇</sup> *<sup>c</sup> <sup>w</sup>*. Customarily, this relationship would emerge from solution of a closure problem, as in Reference [9]. However, this particular dispersive flux relationship is "common practice" [23] (p. 174) textbook material (see also Reference [24], p. 244) and does not affect the form of the NAPL source term, which is the focus of the analysis in this note. Thus, we take the liberty of skipping its derivation.

$$\begin{split} \eta\_{w}\frac{\partial \langle \boldsymbol{c} \rangle^{w}}{\partial t} + \frac{\partial \eta\_{w}}{\partial t} \langle \boldsymbol{c} \rangle^{w} + \eta\_{w} \langle \boldsymbol{\nu} \rangle^{w} \cdot \nabla \left\langle \boldsymbol{c} \right\rangle^{w} - k \eta\_{w} \nabla^{2} \langle \boldsymbol{c} \rangle^{w} - k \nabla \eta\_{w} \cdot \nabla \left\langle \boldsymbol{c} \right\rangle^{w} \\ = \eta\_{w} \mathcal{O} \nabla^{2} \left\langle \boldsymbol{c} \right\rangle^{w} + 2 \mathcal{O} \nabla \eta\_{w} \cdot \nabla \left\langle \boldsymbol{c} \right\rangle^{w} + \mathcal{O} \nabla^{2} \eta\_{w} \left[ \left\langle \boldsymbol{c} \right\rangle^{w} - \mathfrak{c}\_{\text{sat}} \right] - Q. \tag{30} \end{split}$$

*Water* **2019**, *11*, 1525

This can be rearranged with its LHS containing only standard transport terms, and its RHS "source" terms relating to the dissolution of the NAPL:

$$\begin{split} \eta\_{\boldsymbol{w}} \frac{\partial \left< \boldsymbol{c} \right>^{\boldsymbol{w}}}{\partial t} + \eta\_{\boldsymbol{w}} \left< \boldsymbol{\eta} \right>^{\boldsymbol{w}} \cdot \nabla \left< \boldsymbol{c} \right>^{\boldsymbol{w}} - \eta\_{\boldsymbol{w}} \left< \boldsymbol{\mathcal{O}} + \boldsymbol{k} \right> \nabla^{2} \left< \boldsymbol{c} \right>^{\boldsymbol{w}} \\ &= \left( 2 \boldsymbol{\mathcal{O}} + \boldsymbol{k} \right) \nabla \eta\_{\boldsymbol{w}} \cdot \nabla \left< \boldsymbol{c} \right>^{\boldsymbol{w}} + \boldsymbol{\mathcal{O}} \nabla^{2} \eta\_{\boldsymbol{w}} \left[ \left< \boldsymbol{c} \right>^{\boldsymbol{w}} - \boldsymbol{c}\_{\operatorname{sat}} \right] - \boldsymbol{Q} - \frac{\partial \eta\_{\boldsymbol{w}}}{\partial \boldsymbol{t}} \left< \boldsymbol{c} \right>^{\boldsymbol{w}}. \end{split} \tag{31}$$

We simplify this equation by noting that *Q* = *ρ<sup>n</sup> ∂ηw <sup>∂</sup><sup>t</sup>* , and since the water phase solute concentration is dilute, *c <sup>w</sup> <sup>ρ</sup>n*. It follows immediately that *∂ηw <sup>∂</sup><sup>t</sup> c w* <sup>|</sup>*Q*|:

$$\begin{split} \eta\_{\text{\textquotedblleft}} \frac{\partial \left< \boldsymbol{c} \right>^{\text{w}}}{\partial t} + \eta\_{\text{\textquotedblleft}} \left< \boldsymbol{v} \right>^{\text{w}} \cdot \nabla \left< \boldsymbol{c} \right>^{\text{w}} - \eta\_{\text{\textquotedblleft}} \left< \boldsymbol{\mathcal{O}} + \boldsymbol{k} \right> \nabla^{2} \left< \boldsymbol{c} \right>^{\text{w}} \\ &= \left( 2\boldsymbol{\mathcal{O}} + \boldsymbol{k} \right) \nabla \eta\_{\text{\textquotedblleft}} \cdot \nabla \left< \boldsymbol{c} \right>^{\text{w}} + \boldsymbol{\mathcal{O}} \nabla^{2} \eta\_{\text{\textquotedblleft}} \left[ \left< \boldsymbol{c} \right>^{\text{w}} - \boldsymbol{c}\_{\text{\textquotedblleft}} \right] - Q. \end{split} \tag{32}$$

The LHS of this equation resembles the desired advection-dispersion equation (multiplied through by *ηw*). A necessary condition for this to reduce to an ADRE with Sherwood-Gilland source (which depends on only the concentration difference *c <sup>w</sup>* <sup>−</sup> *<sup>c</sup>*sat ) to apply, the equation must simplify so that the mixed gradient term (proportional to ∇*η<sup>w</sup>* · ∇ *c <sup>w</sup>*) on the RHS is eliminated. Doing so requires special assumptions. We consider one such case below.

#### *3.5. Special Case: Advection Domination*

In light of the failure to find an ADRE with Sherwood-Gilland source under general conditions, we follow Reference [17] in considering strongly advective conditions, which we here define as conditions in which *<sup>v</sup> <sup>w</sup>* (2<sup>D</sup> <sup>+</sup> *<sup>k</sup>*) /*L*. As *<sup>η</sup><sup>w</sup>* is less than unity, has a system minimum value above zero, and is a volume averaged quantity, ∇*ηw* <sup>&</sup>lt; *<sup>L</sup>*<sup>−</sup>1. Thus, under sufficiently strongly advective conditions, *η<sup>w</sup> <sup>v</sup> <sup>w</sup>* (2<sup>D</sup> <sup>+</sup> *<sup>k</sup>*) <sup>∇</sup>*ηw*, allowing neglect of the <sup>∇</sup> *<sup>c</sup> <sup>w</sup>* term on the RHS of (32), as its coefficient is dominated by the one on the LHS. In such a case, one may simplify:

$$\frac{\partial \left< \mathbf{c} \right>^{w}}{\partial t} + \left< \mathbf{v} \right>^{w} \cdot \nabla \left< \mathbf{c} \right>^{w} - \left( \mathcal{O} + k \right) \nabla^{2} \left< \mathbf{c} \right>^{w} = \frac{1}{\eta\_{w}} \left[ \mathcal{O} \nabla^{2} \eta\_{n} \left[ \mathbf{c}\_{\text{sat}} - \left< \mathbf{c} \right>^{w} \right] - Q \right]. \tag{33}$$

In order to close the system, we need a second equation to relate *c <sup>w</sup>* and *ηn*. Note also that we are aiming to upscale the dissolution dynamics, and we only have a microscopic equation for *Q*, which is expressed in terms of the local concentration *c*, rather than the needed, upscaled *c <sup>w</sup>*. To solve both problems, we require an ansatz relating *Q* to *c <sup>w</sup>* and *ηw*, for which we consider the Sherwood-Gilland relation. We stress, lest there be any confusion that we "get out what we put in," we are not *deriving* the Sherwood-Gilland formulation from first principles in this section, only considering whether it is consistent with the microscopic formulation presented in Section 2. In the context of our analysis, this means that substitution into (33) ought to transform it to be similar to (7). For our analysis, we combine certain quantities defining the mass-transfer coefficient *K* in (2) into a single parameter, *ω*, and write

$$Q = -\omega \eta\_n^{\beta} \partial \left[ \mathfrak{c}\_{\text{sat}} - \langle \mathfrak{c} \rangle^w \right]. \tag{34}$$

From dimensional considerations, it is clear that *ω* has units of L−<sup>2</sup> and we observe that it plays the role in the above equation of the squared reciprocal of a characteristic length ("boundary layer") for Fickian diffusion. Substituting (12), the definition for *Q*, into (34) yields:

$$\frac{\partial \eta\_n}{\partial t} = -\eta\_n^{\beta} \frac{\omega}{\rho\_n} \left[ \mathbf{c}\_{\text{sat}} - \left< \mathbf{c} \right>^w \right]. \tag{35}$$

This is a second, independent equation in terms of our dependent variables *c <sup>w</sup>* and *ηw*. To arrive at (7), it then suffices to establish that *η<sup>β</sup> <sup>n</sup><sup>ω</sup>* ∇2*ηn*. For *<sup>β</sup>* = 1, a characteristic value in the middle of the range of proposed values of *β* (recall that these ranged from 0.6 to 1.24), (35) can be solved explicitly using the usual strategy for first-order non-homogeneous equations to generate

$$\eta\_{\rm nl}(\mathbf{x},t) = \eta\_{\rm nl}(\mathbf{x},0) \exp\left\{-\int\_{0}^{t} \frac{\omega \mathcal{Q}}{\rho\_{\rm n}} \left[ \langle c \rangle^{w} \left( \mathbf{x}, \tau \right) - c\_{\rm sat} \right] d\tau \right\}.\tag{36}$$

By differentiating (36), we observe that

$$\nabla^2 \eta\_n = \eta\_n \omega \left[ \omega \left( \frac{\varrho}{\rho\_u} \right)^2 \left( \int\_0^l \nabla \left< \boldsymbol{\varepsilon} \right>^w (\mathbf{x}, \boldsymbol{\tau}) d\tau \cdot \int\_0^l \nabla \left< \boldsymbol{\varepsilon} \right>^w (\mathbf{x}, \boldsymbol{\tau}) d\tau \right) + \left( \frac{\varrho}{\rho\_u} \right) \int\_0^l \nabla^2 \left< \boldsymbol{\varepsilon} \right>^w (\mathbf{x}, \boldsymbol{\tau}) d\tau \right]. \tag{37}$$

Recalling that *L* is the length scale of the averaging volume, we see that

$$\nabla^2 \eta\_n \approx -\eta\_n \omega \left(\frac{1}{L^2}\right) \left[\omega \left(\left(\partial t\right)^2 \left(\frac{c\_{\rm sat}}{\rho\_n}\right)^2 + \left(\partial t\right) \left(\frac{c\_{\rm sat}}{\rho\_n}\right)\right] \tag{38}$$

$$<\langle \eta\_n \omega \left[ (\omega \otimes t)^2 \left( \frac{c\_{\rm sat}}{\rho\_n} \right)^2 + (\omega \otimes t) \left( \frac{c\_{\rm sat}}{\rho\_n} \right) \right] \rangle \tag{39}$$

where the last inequality follows because the length scale of pore scale diffusion is small relative to the length scale of volume averaging. Because NAPL are sparingly soluble, it follows that *<sup>c</sup>*sat *<sup>ρ</sup><sup>n</sup>* 1. We may understand *ω*D*t* as an expected number of boundary layer length scales traversed by a solute particle due to diffusion since source emplacement. As long as (*ω*D*t*) *c*sat *ρn* 1, it follows that *<sup>η</sup>n<sup>ω</sup>* ∇2*η<sup>n</sup>* and thus that

$$\frac{\partial \left< \mathbf{c} \right>^{w}}{\partial t} + \left< \mathbf{v} \right>^{w} \cdot \nabla \left< \mathbf{c} \right>^{w} - \left( \mathcal{O} + k \right) \nabla^{2} \left< \mathbf{c} \right>^{w} = \frac{\omega \eta\_{n} \mathcal{O}}{\eta\_{w}} \left( \mathbf{c}\_{\text{sat}} - \left< \mathbf{c} \right>^{w} \right). \tag{40}$$

Under these specific conditions, to the extent that *η<sup>w</sup>* can be treated as approximately constant (for instance, if everywhere *η<sup>w</sup> ηn*), we do recover the ADRE with Sherwood-Gilland source term.

#### **4. Conclusions**

In this note, we consider the dissolution of residual NAPL in saturated porous media by means of a nonstandard volume averaging analysis. Our approach applies a variety of strategies (including direct application of boundary conditions, a physics-preserving domain simplification, elimination of coefficient-dominated terms, a geometric lemma, and some other physically-grounded approximations) to simplify the phase boundary jump discontinuity terms that arise under superficial volume averaging and perform analysis using only the volume averaged equation. This is in contrast with the traditional volume averaging approach which entails developing a separate equation for concentration fluctuations which contains volume averaged terms such as "sources," and formally closing the system by solving one or more boundary value problems to determine the functional form of the dependence of the fluctuations on the volume averaged quantities and then substituting these back into the volume averaged equation. The analysis here is radically simplified.

The particular problem we consider is whether Sherwood-Gilland empirical models are consistent with the micro-scale physics. Existing volume averaging analyses have led to equations that look much unlike the ADRE with a Sherwood-Gilland-type source, so this was a topic worth considering explicitly. The only attempt to recover a Sherwood-Gilland formula that we are aware of in the previous literature is due to Reference [17]: the form is recovered by applying four separate scale constraints that are defined in terms of the formal solutions to fluctuation boundary value problems formulated as part of the volume averaging process. Unfortunately, it is difficult to predict a priori whether these constraints

are satisfied. Above, we derive the volume-averaged transport equation using the simplified volume averaging approach and find it is not an ADRE with a Sherwood-Gilland source term. We subsequently consider a physical restriction of advection-dominated behavior, but unlike Reference [17], it only requires a single restriction on the strength of advective transport relative to dispersive transport. We explicitly corroborate usage of the Sherwood-Gilland source model with the ADRE under such conditions when *β* = 1.

**Funding:** This research received no external funding; SKH holds the Helen Unger Career Development Chair in Desert Hydrogeology.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Appendix A. Nomenclature**

Here, we summarize the algebraic symbols and special operators introduced in this note. Phase subscripts or superscripts *i* and *j* are placeholders for any of the three phases in the system: *n* for NAPL, *s* for solid, and *w* for water. Bold is used to denote vector and matrix quantities (with the indicated dimensions applying to each component), with ordinary-weight Roman text denoting scalars. Dimensions are encoded in accordance with the SI standard: M for mass, L for length, and T for time, with a one indicating a dimensionless quantity.



#### **Appendix B. Derivation of a Geometric Lemma**

This is a simplification of a "geometric lemma" presented by Whitaker on p. 17 of [9]. The derivation follows from use of a conceptual simplification of the domain—eliminating all NAPL-solid interfaces by imagining an infinitesimally thin water layer between the NAPL and the solid—that we show does not alter the flow or transport in the system or the volume averaging mathematics. Physically, because of the no-slip boundary conditions, there is no flow in the new notional, infinitesimally-thin water layers. Because any solute in these layers is at *c*sat, no net NAPL-water mass transfer takes place there. And because the "new" water has infinitesimal volume, it represents a negligible source of mass into the existing pore water phase at the former NAPL-water-solid three-way interfaces. Mathematically, this simplification eliminates the interface Γ*ns*, therefore eliminating terms of the form <sup>1</sup> V <sup>Γ</sup>*ws f nwsdAws*, for arbitrary *f* , when the appropriate volume averaging theorem is applied. Instead, new (collocated) interfaces Γ∗ *nw* and Γ<sup>∗</sup> *ws* come into being. However, <sup>1</sup> V Γ∗ *ns <sup>f</sup> <sup>n</sup>nwdAnw* <sup>+</sup> <sup>1</sup> V Γ∗ *ws <sup>f</sup> <sup>n</sup>wsdAws* = 0 for any *<sup>f</sup>* , because <sup>Γ</sup><sup>∗</sup> *nw* and Γ<sup>∗</sup> *ws* are the same surface, and the two integrals have opposite-directed normal vectors. Thus, the mathematics of volume averaging are unchanged.

Having established that we are always justified in assuming no NAPL-solid interface, the geometric lemma follows from application of volume averaging theorem (4) to an indicator function, *Iw*, representing the presence of the water phase:

$$
\langle \nabla I\_{\rm w} \rangle = \nabla \left\langle I\_{\rm w} \right\rangle + \frac{1}{\mathcal{V}} \int\_{\rm uw} \mathfrak{n}\_{\rm wn} dA\_{\rm wn} + \frac{1}{\mathcal{V}} \int\_{\rm uw} \mathfrak{n}\_{\rm uw} dA\_{\rm ws} \tag{A1}
$$

noticing that everywhere (except for jump discontinuities) *Iw* is a constant, so ∇*Iw* = 0 and *Iw* = *ηw*. Thus,

$$\nabla \eta\_{\overline{\nu}} = -\frac{1}{\mathcal{Y}} \int\_{\Gamma\_{\text{uv}}} \mathfrak{n}\_{\text{uv}} dA\_{\text{uv}} - \frac{1}{\mathcal{Y}} \int\_{\Gamma\_{\text{uv}}} \mathfrak{n}\_{\text{uv}} dA\_{\text{uv}}.\tag{A2}$$

We may apply identical analysis to the solid phase:

$$\nabla \eta\_{\mathbb{S}} = -\frac{1}{\mathcal{Y}} \int\_{\Gamma\_{\mathrm{sr}}} \mathfrak{n}\_{\mathrm{sn}} dA\_{\mathrm{sn}} - \frac{1}{\mathcal{Y}} \int\_{\Gamma\_{\mathrm{sw}}} \mathfrak{n}\_{\mathrm{sw}} dA\_{\mathrm{sw}} \tag{A3}$$

but notice that because porosity is constant, ∇*η<sup>s</sup>* = 0. Furthermore, because we have notionally eliminated the NAPL-solid interface by imagining an infinitesimal layer of saturated, immobile water phase between them,

$$\frac{1}{\mathcal{V}} \int\_{\Gamma\_{sw}} \mathfrak{n}\_{sw} dA\_{sw} = 0.\tag{A4}$$

Noting that <sup>1</sup> V <sup>Γ</sup>*ws <sup>n</sup>wsdAws* <sup>=</sup> <sup>−</sup> <sup>1</sup> V <sup>Γ</sup>*sw <sup>n</sup>swdAsw* = 0, it follows immediately from (A2) that

$$\frac{1}{\mathcal{V}} \int\_{\Gamma\_{\rm uv}} \mathfrak{n}\_{\rm uv} dA\_{\rm uv} = -\nabla \eta\_{\rm uv}. \tag{A5}$$

#### **References**


c 2019 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Water* Editorial Office E-mail: water@mdpi.com www.mdpi.com/journal/water

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18