Dispersion in doublet-type flows through highly anisotropic porous formations

Abstract Steady doublet-type flow takes place in a porous formation, where the log-transform $Y = \ln K$ of the spatially variable hydraulic conductivity $K$ is regarded as a stationary random field of two-point autocorrelation $\rho _Y$. A passive solute is injected at the source in the porous formation and we aim to quantify the resulting dispersion process between the two lines by means of spatial moments. The latter depend on the distance $\ell$ between the lines, the variance $\sigma ^2_Y$ of $Y$ and the (anisotropy) ratio $\lambda$ between the vertical and the horizontal integral scales of $Y$. A simple (analytical) solution to this difficult problem is obtained by adopting a few simplifying assumptions: (i) a perturbative solution, which regards $\sigma ^2_Y$ as a small parameter, of the velocity field is sought; (ii) pore-scale dispersion is neglected; and (iii) we deal with a highly anisotropic formation ($\lambda \lesssim 0.1$). We focus on the longitudinal spatial moment, as it is of most importance for the dispersion mechanism. A general expression is derived in terms of a single quadrature, which can be straightforwardly carried out once the shape of $\rho _Y$ is specified. Results permit one to grasp the main features of the dispersion processes as well as to assess the difference with similar mechanisms observed in other non-uniform flows. In particular, the dispersion in a doublet-type flow is observed to be larger than that generated by a single line. This effect is explained by noting that the advective velocity in a doublet, unlike that in source/line flows, is rapidly increasing in the far field owing to the presence there of the singularity. From the standpoint of the applications, it is shown that the solution pertaining to $\lambda \to 0$ (stratified formation) provides an upper bound for the dispersion mechanism. Such a bound can be used as a conservative limit when, in a remediation procedure, one has to select the strength as well as the distance $\ell$ of the doublet. Finally, the present study lends itself as a valuable tool for aquifer tests and to validate more involved numerical codes accounting for complex boundary conditions.


Introduction
We consider a steady flow taking place in a horizontally unbounded, three-dimensional domain Ω of thickness D. The velocity field is generated by a doublet (a pair of injecting/pumping lines of given strength Q [L 2 /T]). The source and the sink are at (− /2, 0, x 3 ) and at ( /2, 0, x 3 ), respectively, where |x 3 | ≤ D/2 (figure 1a). This configuration is typical of field-scale procedures to determine the properties of aquifers or in in situ remediation (typically pump and treat) strategies. Common to these methodologies is the injection at the well of a solute (either passive or reactive), and the recovery of the flux-averaged concentration (the so-called breakthrough curve, BTC) at the pumping well (figure 1b). The pumped water may be eventually used (after treating) to recharge. The question of relevance here is whether a doublet represents a realistic model for a system of injecting/pumping wells. It is well known (Dagan 1978) that this is authorized when the well radius r w is much smaller than the well length. Because r w ∼ O(1-10 cm), whereas the aquifer thickness is O(1-10 m), replacing a well with a line of singularity is a reasonable approximation. As a consequence, the theoretical study of transport in a doublet-type flow becomes of definite interest for applications. The drawback attached to the doublet is the non-uniformity of the flow field, which makes the simulation (and the successive interpretation) of the BTC very complicated. As it will be clarified later on, such a difficulty is tremendously enhanced by the heterogeneity of the aquifer, and this is (partly) the reason for the very limited analytical studies on such a topic.
The simplest approach to the problem at hand consists of regarding the formation as homogeneous with constant conductivity (a comprehensive overview of the existing analytical solutions can be found in Bruggeman 1999). In this case, flow is characterized by streamlines connecting the source with the sink. In particular, Koplik, Redner & Hinch (1994) have shown that the homogeneous advective field: u = u 0 r/r d−1 (where d = 2, 3 is the space-dimensionality) greatly influences transport causing, in particular, huge differences in the mass arrivals along different streamlines and ultimately producing a persistent tailing (modulated by the scaling parameter u 0 ) in the BTC.
However, natural porous formations are, as a rule, heterogeneous, with K varying in the space by several orders of magnitude (Rubin 2003). Especially in sedimentary formations, heterogeneity manifests as elongated inclusions (lenses) of different permeabilities, which result in a layered pattern. This sets up de facto 'shortcuts' connecting the lines through high conductivity zones, which lead to BTCs earlier than in an homogeneous medium (Kurowski et al. 1994). As a consequence, dispersion is tremendously augmented, as has been detected in a few field-scale transport experiments (Fernández-Garcia, Illangasekare & Rajaram 2004;Ptak, Piepenbrink & Martac 2004).
To account for its erratic variations and the associated uncertainty, it is customary to model the hydraulic conductivity K as a random field. Thus, the log-conductivity Y = ln K is regarded as stationary and normal, defined completely by the geometric mean K G = exp( Y ) (hereafter, the symbol will denote the ensemble average operator), variance σ 2 Y and two-point autocorrelation ρ Y . The latter is anisotropic, with the horizontal integral scale, I, larger than the vertical one, I v .
Transport in doublet-type flows, where the spatial variability of the hydraulic conductivity is accounted for, has been scarcely studied, its theoretical and practical importance notwithstanding (Kurowski et al. 1994). With the exception of a few analytical studies Zech et al. 2018), to our knowledge, there are only numerical studies simultaneously accounting for the heterogeneity of the aquifer and the non-uniformity of the flow pattern (see Bianchi et al. 2011, and references therein).
Q Q x 3 x 1 However, the strong coupling between the spatially variable hydraulic conductivity and the non-uniformity of the flow poses serious numerical issues. Very dense grids are required in the regions of the flow domain where streamlines converge, which therefore prevent the ability to achieve accurate solutions (especially for highly anisotropic formations). Instead, analytical tools lead to simple (i.e. closed form) solutions. These provide explicit relationships between the input parameters and the model output, therefore giving physical insight into the problem, without resorting to computationally heavy (sometimes prohibitive) numerical simulations. We aim to investigate the dispersion mechanism in a doublet flow by means of spatial moments. While the dispersion of both passive and reactive solutes in uniform mean flows has been studied extensively in the past (see, e.g. Dagan 1984;Cvetkovic & Dagan 1994), much less has been done at the conditions typical of source flows (a review of recent results can be found in Severino, Leveque & Toraldo 2019) and we are not aware of any theoretical derivation of spatial moments for a doublet-flow configuration. Before proceeding, we wish to clarify the difference between the present study and that of  (similar conclusions can be drawn for the study of Zech et al. 2018). In these studies, the aim was: (1) to compute the BTC at the sink given the input at the source; and (2) to assess how the heterogeneity of the formation impacts the shape of the BTC. However, such a methodology does not allow for quantification of the dispersion process in the intermediate zone between the two lines (figure 1). In fact, the approach of  relies on the quantification of the probability distribution function of the arrival times of solute particles at the source, and therefore it does not account for the entire history of the spatial dispersion between the release and recovery of the solute. From the standpoint of applications, an important implication concerns the identification of the heterogeneous nature of the conductivity. In fact, unlike the BTC, spatial moments are affected per se by the distribution of the advective velocity between the injecting and pumping lines. As a consequence, the match between the theoretical second-order spatial moments (derived in the present study) and the experimental ones leads to a robust identification of the conductivity field.
The paper is organized as follows: after recalling the general expression of spatial moments, we employ an analytical approximate solution, which ultimately leads to a simple (closed form) expression for the first-and the second-order spatial moments. Then, we discuss the main differences in the dispersion mechanism with transport taking place in other non-uniform flows and we end with our concluding remarks.

The transport problem
A passive solute is injected injected in Ω through the sink (figure 1) at a concentration C 0 ≡ C 0 (a). We assume that Ω is initially solute free. Spatial dispersion can be quantified by means of the first order spatial moments (where the constants M = ϑ da C 0 (a) and ϑ are the injected mass and the porosity, respectively). For a pulse of constant concentration C 0 , moments (2.1)-(2.2) are written as (Dagan 1989). These expressions rely upon two fundamental hypotheses.
(i) Pore scale dispersion D is neglected. To discuss the feasibility of this approximation, it is instrumental to recall that advection results in the fragmentation of the plume into spots moving quicker/slower than the mean velocity. This produces a concentration gradient between adjacent spots, which ultimately leads to a Fick-type mixing. While this mechanism determines a reduction (dilution) of the local concentration (Fiori & Dagan 2000), the matter here is whether D is also influential on second-order moments (2.2). This can be easily established by comparing the characteristic advection time, t a , relative to the characteristic mixing time, t D .
In particular, if t a is larger than t D , then local dispersion is expected to have a 'non-negligible' impact, otherwise it can be ignored. Because typical values of the above characteristic times are such that t a /t D 1 (see, e.g. Dagan 1989;Rubin 2003), one concludes that pore-scale dispersion can be neglected in the majority of the applications (see also Severino, Santini & Sommella 2011;Severino et al. 2012a;Zech et al. 2018); (ii) Ergodicity, the condition which de facto allows for replacing spatial averages with their statistical counterparts in (2.3a,b), is fulfilled within the present study. More precisely, to invoke the ergodic argument, the thickness D is required to be much larger than the vertical integral scale I v (for details, see Dagan 1989, ch. 1.10). Rubin 2003), ergodicity is met in most of the real-world situations. As indirect consequence, we can also neglect the effect of the lower/upper boundaries on flow and transport, therefore assuming Ω ∼ = R 3 .
Overall, determining the higher-order statistics (typically skewness and kurtosis) would be worthwhile to comprehensively capture the behaviour of the plume. However, logistic as well as economic limitations in the experiments (see, e.g. Fernández-Garcia et al. 2004) do not allow for acquiring the proper dataset needed to compute skewness and kurtosis. For this reason, stochastic tools for the analysis of transport experiments taking place in randomly heterogeneous porous formations generally employ second-order statistics. The usefulness of adopting moments (2.3a,b) to characterize the plume spreading is that they depend directly upon the mean X and the fluctuation X = X − X of the trajectory X of a fluid particle which, in turn, are determined by the kinematic equation Thus, central to quantifying transport is the velocity field u. The latter is not exactly solvable, even for the case of a mean uniform flow (an extensive treatment on this topic can be found in Dagan 1989). In the following, we adopt an approximation for the flow field leading to a simple solution for the transport problem, thereby avoiding heavy (sometimes prohibitive) numerical simulations. Although approximate, the solution for the flow field keeps the salient features of the problem.
2.1. Approximate solution of the velocity field The system of injecting/pumping wells is replaced by a point-like source/sink (dipole). Such an approximation is valid when the radius r w of the wells is much smaller than the distance between them (figure 1a). Because r w ∼ O(1-10 cm), whereas ∼ O(1-10 m), this is a reasonable approximation.
We adopt a first-order approximation in the fluctuation Y = Y − Y , which has been shown to be quite a robust hypothesis (see, e.g. Firmani, Fiori & Bellin 2006). Hence, the flow variables (specific energy, flux, etc.) can be expanded in an asymptotic series of Y and, for each of them, up to the σ 2 Y -order, one retains the leading-order term (mean) and the first-order (fluctuation) approximation (an approach similar to the 'frozen field' approximation in turbulent flows). Owing to its importance for the transport problem, we focus in the following upon the velocity field. In particular, at the leading order, one has (2.5) , where x r = x 2 1 + x 2 2 is the magnitude of the vectorial position x r ≡ (x 1 , x 2 ) in the horizontal plane. Thus, at the leading order, the velocity field u (0) is de facto two-dimensional, although transport is still three-dimensional (the fluctuation X obeys the three-dimensional (2.12)).
To simplify the computational burden, we deal with a formation having the anisotropy ratio λ = I v /I much less than one. This heterogeneous structure is typical of those (highly anisotropic) formations where the fluctuation Y along the vertical direction is much larger than that in the horizontal plane (see, e.g. Sudicky 1986;Zinn & Harvey 2003). This authorizes the replacement of the vertical autocorrelation with a 'white noise' signal, Fiori, Indelman & Dagan 1998). This approximation has two fundamental implications.
(i) The fluctuation of the velocity is written as (2.6) . In fact, adopting (2.6) was found to yield accurate results for the dispersion mechanism already for λ ≤ 0.2 ). More generally, it has been recently shown that the approximation (2.6) becomes an exact result for λ → 0 (stratified formation), even in unsteady non-uniform flows (Severino & Cuomo 2020). Because numerous (typically sedimentary) formations are such that I v ≤ I/10 (see e.g. tables 2.1-2.2 in Rubin 2003), the above approximation is applicable to many situations of practical interest. (ii) Transport variables become stationary along the vertical and, therefore, they depend only on x r .
Based on this, we proceed with the computation of moments and, in particular, with the quantification of the dispersion pattern in the zone between the source and the sink.

Longitudinal dispersion along the central strip delimited by the source/sink
We consider a solute injected at the source (figure 1a). For t > 0, the plume is advected outwards and it is distorted owing to the Y-heterogeneity (figure 1b). It is seen from (2.3a,b) that to characterize the migration of the plume, it suffices to compute the second-order statistics of the trajectory X ≡ (X 1 , X 2 , X 3 ) of a fluid particle. In particular, because the average velocity is identical to that in an homogeneous medium of conductivity K G , the mean trajectory X results from (2.4a,b) as follows (chain-rule of derivation): Then, accounting for the velocity field (2.5) leads to the nonlinear first-order equation, which is transformed into a linear one after introducing Y = X 1 2 and Z = ln X 2 , i.e.
The second part of (2.9) corresponds, for any C ∈ R, to a streamline originating at the source. It is easy to show that all the streamlines belong to a bundle of circles (Apollonius), passing by the sink and the source, with centres and radii given by (0, C) and C 2 +¯ 2 , respectively. In particular, the trajectories within a tiny strip between the source and the sink are such that X 2 X 1 (figure 1b) and, therefore, the longitudinal trajectory X 1 ≡ X 1 (t) can be obtained by expanding in Taylor series the right-hand side of (2.4a,b), with u (0) 1 given by the first part of (2.5). This leads to (2.10) Solving the above Cauchy problem provides Thus, the central normalized trajectoryX is an increasing function of the time, and it takes time t(1) = 4/3 t c to reach the sink.
We are now in position to compute the trajectory variance X m X n , which in turn is related to the fluctuation X by means of (Indelman & Rubin 1996). We restrict the analysis to the trajectories along the x 1 -direction, because they are of most interest for the spreading mechanism. In fact, the spreading of the plume is advection dominated (see, e.g. Severino, Cvetkovic & Coppola 2005) and therefore dispersion is enhanced along streamlines where the advection velocity is larger. To establish which is the streamline with higher velocity, in the spirit of the employed perturbation approach, one can look at the leading-order approximation u (0) of the flow field. Toward this aim, it is convenient to switch to the new variables (φ, θ ) that are related to x r as (where, in particular, θ is the attack angle of the outgoing trajectory). Thus, in the new framework, the magnitude of the velocity results from (2.13a,b) as |u (0) | ∼ cosh φ + cos θ, which shows that the streamline exhibiting the largest dispersion is the one along the x 1 direction (θ = 0). For this trajectory, the probability that X 2 /X 1 ∼ O(1) is quite small, and concurrently the fluctuation X 1 results from (2.12) as (2.14) The solution of (2.14) is , (2.15) and thus the variance X 11 (t) = X 2 1 (t) along the central streamline reads as , (2.16) where, for the sake of brevity, we have replaced u ] is the covariance of the velocity field u that is computed by the approximation (2.6) as ( 2.17) Then, substitution into (2.16) leads to It is worth noting that for a mean uniform flow, i.e. u (0) 1 = const, one recovers from (2.18) the same expression (see (3.10) in Dagan 1987) valid for a groundwater flow. For convenience, we switch to the coordinate X 1 ≡ X 1 (t) by virtue of the kinematic (

2.19)
Introducing the new variables u = α − β and v = β enables one to write (2.20) Thus, if the leading-order term u (0) 1 makes it possible to express in closed form the inner quadratures in (2.20), the computation of X 11 is reduced to the evaluation of a single integral. This latter is then carried out (either analytically or numerically) once the shape of the autocorrelation function ρ Y is selected. This is just the case for the problem at hand, and the final result reads as where we have set In what follows, we shall discuss the combined effect upon solute dispersion of the heterogeneity of the medium and the flow configuration.
We wish to illustrate and discuss the results achieved so far. To show, in a simple manner, the main features of the dispersion between the two lines, we regard the fluctuation Y as a white noise signal (in analogy with the δ-approximation of the energy spectrum in the theory of homogeneous turbulence). This simplifies tremendously computing the right-hand side of (2.21), and the final result is (for convenience, we have taken¯ as the characteristic length scale). It is seen that X 11 grows monotonically with the scaled travel distanceX ∈] − 1, 1[ (and ultimately with the time t). Such a behaviour has a straightforward mechanical explanation by recalling that the variance (3.1) is computed from the outset of the advective velocity field. In fact, close to the release zone, transport is affected by the boundary condition and, concurrently, dispersion is negligible there (it is reminded that the pore-scale dispersion is neglected). Instead, as the travel distanceX increases, the random fluctuations of the velocity field overtake over and over and this justifies the increase of X 11 . Ultimately, the latter becomes unbounded owing to the fact that the velocity at the sink becomes singular. The transitional behaviour from the near (i.e.X −1) to the far (i.e.X 1) field is in agreement with  and Severino (2011), who have investigated transport generated by source-type flows.
In figure 2, we have depicted the scaled longitudinal moment X 11 /(I¯ σ 2 Y ) obtained from (2.21) by considering the exponential ρ Y (x) = exp(−x/I), i.e. 2) (a similar result, although much more cumbersome, is obtained by adopting Gaussian ρ Y ) as a function of the distance R = tan[π( X 1 + )/(4 )]. The utility of the mapping X 1 ∈] − , [ → R ∈ [0, ∞[ relies upon the fact that, in this way, one can compare with similar results in source-/line-type flows. The most evident feature detected from figure 2 is the increasing dispersion (for a given R) with the smallerĨ, ultimately reaching the upper bound corresponding to a stratified formation (red line). In fact, to pass through a low conducting inclusion, a fluid particle has to cover a large distance, therefore increasing the correlation length, and hence the velocity covariance C u 11 . As a consequence, the second-order moment X 11 will also grow owing to its dependence (see (2.16)) upon C u 11 . Such an increase is enhanced by the smallestĨ values, because a formation with I implies a longer correlation distance. In addition to the computational simplification, the expression (3.1) lends itself as an upper bound for the dispersion mechanism, which in turn can be conveniently  Figure 2. The normalized trajectory variance X 11 /(I¯ σ 2 Y ) for exponential autocorrelation (black lines) as a function of R and a few values of the non-dimensional integral scaleĨ = I/ . For comparison purposes, the trajectory variances for a line (blue) and source (green) flow  have also been included. Finally, the red line depicts (3.1), which corresponds to a stratified formation.
accounted for when designing remediation strategies. More precisely, one can select the injecting/pumping flow rate Q as well as the distance such that dispersion at a certain distance x ∈] − , [ is less than a given (regulatory) threshold. It is therefore clear that dealing with a stratified formation, and thus by accounting for the simple expression (3.1), would lead to a conservative estimate.
At the other extreme, line (and a fortiori source) flow constitutes a lower bound for dispersion in a dipole. This is clearly observed in figure 2, where we have also depicted the trajectory variance for a line and source flow, i.e.
(3.5) . The different behaviour of (3.2) as compared with the latter is explained by noting that the velocity, i.e. u(R) = (2R) −1 , in a line-type flow and that, i.e. u(R) = (4R 2 ) −1 , in a source-type flow are decreasing with the scaled distance R (figure 3), whereas in a dipole flow, the velocity u(R) = (π/4) 2 arctan R (π/2 − arctan R) (3.6) (which is obtained by replacing X 1 → (4 /π) arctan R − in (2.10)) is unbounded at the source (i.e. R = 0) and at the sink (i.e. R → ∞). As a consequence, the fluctuation of the velocity (which, in the spirit of the perturbation expansion employed in the present study, slightly differs from the mean velocity) is unbounded both at the release and at the recovery of the solute, therefore producing a larger dispersion as compared with that detected in source/line flows.

Summary and conclusions
Solute transport in a doublet-type flow through a heterogeneous porous formation has received, in the past, very little attention, its importance in the applications notwithstanding. The main difficulty, which renders this problem very complex, is the coupling between the spatially variable hydraulic conductivity and the strong non-uniformity of the flow field. Even worse, the numerical approach does not seem computationally affordable for highly anisotropic formations of a three-dimensional structure.
In the present study, we have focused on the dispersion process occurring in the strip delimited by the source and the sink ( figure 1). An analytical (closed form) expression for the longitudinal spatial moment has been derived. It has been achieved by adopting a few approximations: (i) the log-conductivity is a stationary random field of axisymmetric anisotropy; (ii) a perturbation solution for the flow field is sought; and (iii) pore-scale dispersion is neglected. Dealing with a highly anisotropic formation simplifies considerably the computation of the statistics of the flow field and concurrently the evaluation of the longitudinal spatial moment. The latter is expressed in closed form for any autocorrelation function ρ Y , which can be straightforwardly evaluated once the shape of ρ Y is specified. The major difference between transport in a doublet and other similar (typically source and line type) flows is that dispersion in the former is larger than in the latter. This effect arises from the advective velocity (3.6) which, unlike that in source/line flows, is rapidly increasing in the far field owing to the presence there of the singularity.
In addition to the theoretical aspects, expression (2.21) is also useful for the applications. In fact, by taking the limit λ → 0 (stratified formation) provides an upper bound for the dispersion mechanism (figure 2). Thus, it can be adopted as an (conservative) estimate to select the strength Q and the distance of the doublet, when designing remediation strategies. In addition, it is shown that dispersion arising from a single injecting well ) constitutes a lower bound for the same mechanism. This is clearly owing to the fact that, unlike the doublet, the advective velocity decays like R −1 . Even such a lower bound finds use in the applications concerning the use of chemicals to neutralize dissolved pollutants. In this case, one has to maximize the dispersion to enhance reactions (Di Dato et al. 2018), and therefore the design accounting for the lower bound certainly provides another conservative estimate of the efficiency of the procedure.
Although the study has focused on dispersion of a passive scalar, it can be extended straightforwardly to reactive solutes along the lines developed by Cvetkovic & Dagan (1994). Even in this case, if one is interested in global quantities (like BTCs or spatial moments), then pore-scale dispersion can be neglected. In fact, the involved reaction equation very often overtakes local mixing (Severino, Dagan & van Duijn 2000;Severino et al. 2012b). Nevertheless, to assess the efficiency of certain remediation procedures (specifically those relying upon the concentration values), the above global quantities may not result in an appropriate tool. In this case, the computation of the probability distribution function of the concentration value(s) would be preferable and dispersion can not be neglected anymore. Finally, the present study lends itself as a benchmark for validating more involved numerical codes.