Macrodispersion and Recovery of Solutes and Heat in Heterogeneous Aquifers

The recovery efficiency of aquifer storage systems with radial flow fields are studied for heterogeneous aquifers. Macrodispersion, arising from spatially heterogeneous hydraulic conductivity, is modeled as a scale‐dependent mechanical dispersion process. Approximate solutions for the recovery efficiency as a function of local dispersion and macrodispersion parameters, the injection‐extraction rate Q $Q$ and duration T $T$ , and storage cycle count, are derived and validated against numerical simulations. If macrodispersion dominates and the macrodispersion coefficient scales linearly with distance, the recovery efficiency is independent of both Q,T $Q,T$ . For sublinear and superlinear scalings, recovery increases and decreases respectively if Q,T $Q,T$ increases. However, if local dispersion dominates, increasing Q,T $Q,T$ always increases recovery. As macrodispersion becomes increasingly dominant with scale, the recovery efficiency may be a nonmonotonic function of Q,T $Q,T$ , with a maximum. In homogeneous aquifers, nonmonotonicity does not occur for 1D and 2D radial flow, but occurs for 3D radial flow fields only as a function of T $T$ , not Q $Q$ . These methods may also be used for fitting local dispersion and macrodispersion parameters with push‐pull tests using recovery data, with advantages in scope of applicability and ease of data acquisition and interpretation, compared to existing push‐pull test methods, which fit to breakthrough curves and do not consider macrodispersion. Furthermore, characterizing macrodispersion with push‐pull tests may be advantageous over methods that use observation wells, as observation well placement may be challenging in highly heterogeneous aquifers. The results show that the macrodispersion parameters are not innate aquifer hydraulic properties, as their values vary with flow field geometry.


Problem Statement
We use solute transport terminology, but note that the analysis is mathematically analogous and fully applicable to heat transport (Lee, 1998): molecular diffusion in solute transport is equivalent to thermal diffusion in heat transport, while mechanical dispersion applies similarly to both heat and solute transport. Linear chemical and thermal retardation, that occur when solutes adsorb onto the pore matrix, and when the heat capacity of the solid and liquid phases differ, respectively, are implicitly considered as both can be captured through a linear rescaling of time.
Consider the geometry of flow fields originating from wells in homogeneous aquifers with negligible background flow. The flow fields investigated in this study are the -dimensional spheres. A commonly studied conceptualization of ASR and ATES is that of a fully penetrating well in a confined aquifer, corresponding to a cylindrical flow field. Since the well is fully penetrating, the cylindrical flow field is equivalent to a disc-shaped field, which is a 2D sphere. Figure 1 illustrates that other radially axisymmetric macroscopic flow field configurations are possible in a 3D space: unidirectional planar flow is a 1D radial flow field, while spherical flow from a point source is a 3D radial flow field. 10.1029/2021WR030920 3 of 21 In heterogeneous aquifers, flow around wells can be described as radially axisymmetric flow around a point source, perturbed by velocity field fluctuations at a smaller spatial scale . The flow velocity ( 1, …, ) in a -dimensional heterogeneous flow field can therefore be approximated as the sum of the macroscopically averaged axisymmetric velocity ( ) and a spatially dependent fluctuation ( 1, …, ) (Dagan, 1988;Indelman, 2001) 2 is the radial position and the 1, …, are the Cartesian coordinates. For small-scale fluctuations with zero mean, the macroscopic flow field is radially aligned, and the macroscopically averaged axisymmetric velocity obeys where is the injection rate, is the porosity of the aquifer, Γ( ) is the gamma function, and are arbitrary variables. is the spatial dimensionality of the macroscopically averaged radial flow field (e.g., = 1 for planar, = 2 for disc, and = 3 for spherical flow). is shape coefficient of the injected water front, whose component within square brackets is the surface area of a -dimensional sphere of unit radius. The governing advection-dispersion equation (ADE) that describes the transport of solutes and heat in homogeneous soils is (Tang & van der Zee, 2021) where is the dispersion coefficient, is the molecular diffusion coefficient, is the longitudinal mechanical dispersivity constant, and ℎ is the hydrodynamic dispersion coefficient. Since no transverse dispersion occurs under radial flow in homogeneous aquifers (Gelhar & Collins, 1971), we assume that the effects of transverse mechanical dispersion in heterogeneous aquifers are captured through macrodispersion, which will be introduced later. Lessoff and Indelman (2004) found that varying the transverse dispersivity has only a minor impact on overall dispersion. Hence, in the numerical simulations that follow, the transverse mechanical dispersivity constant is set as equal to to maximize computational efficiency.
We consider a problem in which at the well, a duration of injection rate alternates with the same duration of extraction rate -, in a step-cyclic manner in a -dimensional infinite domain. Assume that the flow field immediately achieves steady-state as the injection switches to extraction or vice versa. The injected volume in = and extracted volume ex = are equal. The water injected by the well contains a solute at some dimensionless reference concentration = 0 , while the ambient water surrounding the well initially contains solute at = 0 . The well injects a total mass = 0 = 0 in of solute, and extracts ex of water, recovering of solute. The recovery efficiency of solute is = . (6)

Numerical Model
The problem described above is solved with MODFLOW-2005(Harbaugh, 2005 and MT3DMS (Langevin & Guo, 2006) for flow and solute transport, respectively. The finite difference flow calculations were performed using the pre-conjugated gradient (PCG) package of MODFLOW in block-centered numerical grids. Solute transport in MT3D was solved with the third-order ULTIMATE TVD scheme. A uniform head initial condition was imposed across the numerical domain, and an imposed head boundary condition (magnitude equal to the initial condition) was imposed at the edges of the domain, to simulate the absence of background flow. Imposed flux boundary conditions were used at the well, with the flux uniformly distributed along the well screen, during both the injection and extraction phases. The initial and boundary conditions that describe solute transport are where subscripts and refer to the injection and extraction phases respectively and max is the r-coordinate at the edge of the numerical domain. The zero-gradient boundary condition at the well during the extraction phase implies that the concentration immediately inside the wellbore is equal to that immediately outside, and reflects the assumption that mixing within the wellbore is negligible.
The radial flow fields and aquifer geometries considered in the numerical simulations are planar flow ( = 1 ) in a two-dimensional domain and in a three-dimensional domain, and cylindrical flow fields ( = 2 ) in a three-dimensional domain. Solute transport in these flow fields with three types of heterogeneity structures were modeled: (a) spatially autocorrelated heterogeneity with exponential covariance; (b) white noise heterogeneity; and (c) 5 of 21 layer stratification. In the two-dimensional simulations, a 256 × 256 numerical grid with horizontal and vertical dimensions = 75, = 50 was used. In the three-dimensional simulations, a 128 × 128 × 64 numerical grid was used, with = 75, = 50 , = 10 . Grid cells along the edges of the square (two-dimensional) and cubic (three-dimensional) numerical grids were de-activated so that the active area of the numerical grid resembled a circle (two-dimensional) or sphere (three-dimensional) as closely as possible, to prevent edge effects caused by simulating radial flow in square or cubic domains. For cylindrical flow fields ( = 2 ), the well is fully penetrating along the z-axis. In both 2D and 3D models with planar flow ( = 1 ), the macroscopic flow direction is along the x-axis. Dispersion into confining layers of the aquifer is not explicitly considered in this study, because such scenarios are essentially specific cases of aquifer stratification with high conductivity variances: one can consider high conductivity layers to be aquifer layers, and low conductivity layers to be confining layers.
For the nonautocorrelated (white noise) fields, the log-conductivities were independently drawn at each cell. In the stratified aquifers, 16 layers were generated; the hydraulic conductivity is heterogeneous only along the y-axis for 2D simulations and z-axis for 3D simulations, with no autocorrelation between layers. As such, this model of stratification represents a limiting case of autocorrelated heterogeneity in which the correlation length is infinite within the strata and zero across the strata. For the exponentially autocorrelated fields Cov( ) = 2 exp(− ) , the simulated autocorrelation lengths along each axis are = ∕10 , where i refers to the index of each spatial dimension, to ensure non-trivial structures of heterogeneity within the confines of the numerical grid. No assumption of ergodicity is made, as in practice, for managed aquifer recharge applications such as ATES and ASR, the storage radius of the system may be even smaller than a single integral scale of heterogeneity (e.g., Sommer et al., 2013). For example, in the Netherlands, where most presently operating ATES systems worldwide are located, the storage radius of wells never exceeds 100m, with most being under 50 m (Bloemendal & Hartog, 2018). Real aquifers reported in many studies in the literature (e.g., Sommer et al., 2013;Vereecken et al., 2000;and Fernàndez-Garcia et al., 2005) have horizontal correlation lengths of between 2.5 and 100 m. This means that many aquifer storage systems in practice may not sample a sufficient volume of the aquifer they are located in to approach ergodic conditions. Therefore, our choice of a numerical domain with 10 integral scales appears to be a good middle-ofthe-road choice to show that non-ergodicity does not preclude the applicability of our findings. The exponentially autocorrelated conductivity fields were generated using the RandomFields R package (Schlather et al., 2015). The natural logarithms of the hydraulic conductivity at each cell was drawn from a normal distribution with mean 1 and four different standard deviations = {0.5, 1, 1.5, 2} . The modeled values of are field-realistic values that are spread around those of three experimentally studied aquifers, Cape Cod (

Scale-Dependent Macrodispersion
We define the scale-dependent macrodispersion coefficient as = (e.g., Dagan, 1988), where is the macrodispersivity constant and is a positional scaling exponent. Such macrodispersion scaling (i.e., log ∝ log ) is frequently reported in the literature (e.g., Wheatcraft & Tyler, 1988;Zech et al., 2015). By the principle of conservation of mass, due to the assumed incompressibility of water, the macroscopically averaged radial position ′ of the injected hydraulic front is related to the injection rate ∝ by Thus, the Macrodispersion coefficient at ′ is We account for macrodispersion in the dispersion coefficient by adding the macrodispersion coefficient to the hydrodynamic dispersion coefficient ℎ : Adding the macrodispersion coefficient to the hydrodynamic dispersion coefficient, and using the macroscopically averaged velocity ( ) to model velocity-dependent dispersion terms, is sufficient to successfully describe

Water Resources Research
TANG AND VAN DER ZEE 10.1029/2021WR030920 6 of 21 macroscopic spatial moments of plume spreading (Dagan, 1988). The ability to additively combine the local dispersion and macrodispersion terms into an overall dispersion term stems from the observation that in typical situations where the heterogeneity autocorrelation scale is larger than , the macrodispersion term and the local dispersion term are independent (Gelhar, 1986).

Recovery Efficiency Under Scale-Dependent Macrodispersion in Planar Flow Fields
The recovery efficiency in a planar flow field in a homogeneous aquifer subject to only scale-independent mechanical dispersion (i.e., = ) is (Tang & van der Zee, 2021) Here, the recovery efficiency increases as increases.
In the case of Equation 4 subject to = 0 , = 0 , and = 1 , in a planar flow field, under purely scale-dependent dispersion = , the exact solution for the concentration profile under continuous injection is (Pang & Hunt, 2001) where Γ(,) is the upper incomplete gamma function and are arbitrary variables. The mass of solute (i.e., the zeroth spatial moment) outside the hydraulic front during the injection phase can be found by integrating the concentration profile from the front to infinity and multiplying by the porosity : The recovery efficiency of a complete injection-extraction cycle is then approximately given by the following expression from Tang and van der Zee (2021): which is accurate only when dispersion is not too strong relative to advection (i.e., when ≥ 0.7 ) because it uses a boundary layer solution to the ADE.
Remarkably, in contrast to the position-independent mechanical dispersion (Equation 11), the recovery efficiency is independent of the porosity and the well operational parameters if the dispersion coefficient increases linearly with distance from the source, that is, = 1 . In summary, the recovery efficiency is independent of for = 1 , and the recovery efficiency increases as increases if = 0 as implied by Equation 11. Accordingly, assuming that the function mapping to the recovery efficiency is continuous, the recovery efficiency should decrease as increases if 1 . With multiple coexisting dispersion processes (e.g., = + ), whether the recovery efficiency increases or decreases as increase should then depend on whichever process dominates. This is because mechanical dispersion is velocity-dependent, whereas molecular diffusion is not, thus the 7 of 21 -dependence of flow velocity in radial flow fields causes the relative contribution of mechanical dispersion and molecular diffusion to solute losses to change with . For example, in three-dimensional radial flow in a homogeneous aquifer, increases as increases if mechanical dispersion dominates, but the opposite occurs if molecular diffusion dominates (Tang & van der Zee, 2021). Unfortunately, the solution to the ADE for values of other than 0 or 1 (Guerrero & Skaggs, 2010), or for higher dimensional flow field geometries (Hunt, 1998), are too complex to obtain a closed-form solution for the recovery efficiency. Therefore, we seek to obtain an approximate solution for the recovery efficiency that applies to a wide range of scenarios.

Approximate General Solution for Recovery Efficiency
An approximation by Gelhar and Collins (1971) states that the radial concentration profile ( ) in -dimensional flow in a homogeneous aquifer is The approximation is accurate when the average hydraulic front radius ′ ( ) is significantly larger than the width of the dispersed zone, that is, It was previously shown that in radial flow fields, the first spatial moment obtained from the solution of the radial ADE (Equation 4) is able to characterize the first spatial moment of solute concentration in heterogeneous aquifers, when the macrodispersion coefficient is substituted for the hydrodynamic dispersion coefficient (Adams & Gelhar, 1992). Since the first spatial moment of concentration essentially quantifies the total mass of solutes within a certain volume of space, this allows us to approximate the recovery efficiency, by first substituting (Equation 10) for ℎ in Equation 18, and also expressing ′ using Equation 8, which yields for the macroscopically averaged radial concentration profile of a heterogeneous aquifer.
The storage radius is the maximum extent of the hydraulic front achieved during a cycle, and can be found by substituting the injection duration into Equation 8, which yields = ( ) 1 . The mass of solute dispersed beyond the stored volume may be obtained by spatially integrating (20) beyond the storage volume: where ( ) is the volume of a -dimensional sphere of radius .
Zee, 2021), applicable when ≥ 0.7 , may computed as In the case of = 0 , = 0 , = 1 , and = 1 , as previously discussed in Section 3.2, the recovery efficiency (Equation 23) becomes which is independent of ( ) and , in agreement with Equation 16. Comparing obtained from the approximate concentration profile (Equation 24) to that obtained from the exact profile (Equation 16) reveals that the two approaches of calculating are effectively identical for low (macro) dispersion and high recovery scenarios (regions of low in Figure 2).

The implications of Equation 16 and Equation 24
that the recovery efficiency in this case is independent of and is extremely remarkable, therefore we numerically solved this problem (Equation 4 and Equation 10 with = 1 , = 0 , = 0 , and = 1 ) with the pdepe routine of MATLAB for the purpose of additional verification. The obtained with Equation 24 agrees exactly with those numerically obtained for ≤ 0.01 . More importantly, the numerical results confirm that varying and does not affect the recovery efficiency at all for ≤ 0.01 . As increases beyond 0.01, Equation 16 and Equation 24 begin to underestimate the numerical result, to a larger extent as becomes larger (i.e., the smaller becomes) due to the use of boundary layer approximations in the derivations, though the relative underestimation remains well below 10% at ∼ 0.1 (Figure 2). Laboratory and field evidence suggests that in practice is mostly on the order of 0.001-0.01 (discussed further in Section 4.3), with 0.1 being a soft upper bound. Therefore, Equations 16 and 24 are essentially accurate for most situations of practical relevance. Figure 2 also shows that for larger , varying and affects the numerically computed recovery efficiency slightly. Simulations with various combinations of and reveal that as the storage volume in = increases, Equation 16 underestimates the numerically computed recovery efficiency to a larger but still reasonable extent. In summary, the results show that Equation 16 is always larger than 24, and that Equations 16 and 24 closely approximate the numerical results for 0.1 . Furthermore, Equations 16 and 24 lower bound the numerically computed recovery efficiency for cases with atypically large or in . From Equation 23, it follows that if = 0 , = 0 , and 1 , then the recovery efficiency (Equation 23) decreases as increases whereas the opposite is true if 1 , in agreement with our hypothesis in Section 3.2.
Furthermore, from Equation 23 it can be deduced that if = , then for all combinations of = 1 the recovery efficiency is independent of , , and the scenario analyzed above ( = 1 , = 1 ) is only one of many possibilities. Equation 23 is also reducible to the solution for the recovery efficiency in homogeneous aquifers (Tang & van der Zee, 2021), by setting = 0 . Thus, Equation 23 captures the complexity relating to the two empirical macrodispersion parameters and , which have to be fitted from observational data and are somehow related to the aquifer properties and .

Relation Between Aquifer Structure and Macrodispersion Parameters
Both in 3D and 2D domains, when 2 < 2 , solute moments can be effectively described with only two parameters, namely and , and other characteristics of the structure of heterogeneity are mostly inconsequential (Di Dato et al., 2019;Jankovic et al., 2017;Zinn & Harvey, 2003). Therefore, here we relate and to the macrodispersion parameters and .
The macrodispersivity constant is related to the heterogeneity of the aquifer, but is not simple to derive from the statistics of the conductivity distribution (Fernàndez-Garcia et al., 2005), and its validity in representing true mixing is often disregarded. Expressions have been derived for (Dagan, 1988), including expressions specifically concerning cylindrical flow Indelman, 2004;Indelman & Dagan, 1999;Neuweiler et al., 2001) and spherical flow (Indelman & Dagan, 1999). Mostly, these expressions pertain to the two limiting cases and of autocorrelated heterogeneity. In these cited studies, the macrodispersion coeffi-

Water Resources Research
TANG AND VAN DER ZEE 10.1029/2021WR030920 9 of 21 cient is shown to be proportional to the flow velocity and the variance of the log-conductivity: Accordingly, decreases to zero if the aquifer is homogeneous ( = 0 ), and no macrodispersion occurs.
For autocorrelated structures, preferential flow paths may be more persistent than under the absence of autocorrelation, and this affects , which describes the positional scaling of the macrodispersion. Therefore, is a function of the autocorrelation length , which encodes how disordered or connected the structure of the aquifer heterogeneity is (Bianchi & Pedretti, 2017;Trinchero et al., 2008). In aquifers with longer autocorrelation ranges or less disorder, such as perfectly stratified aquifers, should be larger, as preferential flow channels form and persist over a larger spatial expanse. In aquifers with preferential flow channels, the amount of ambient water that the solute-rich injected water encounters increases significantly as the volume injected increases, because the surface area of the plume expands (Figures 3a and 3b), allowing more solutes to disperse into relatively immobile channels (Güven et al., 1985). Furthermore, as stratification represents the extreme case of infinite along the direction of flow, is expected to be larger in stratified than in exponentially autocorrelated aquifers.
In contrast, nonautocorrelated heterogeneity gives rise to some degree of "surface roughness" at the hydraulic front without persistent preferential flow (Figure 3c), with a relatively small surface area of interaction with the ambient water. Setting = 0 in Equation 23 reduces the macrodispersion term to a scale-independent mechanical dispersion term. Since macrodispersion under non-autocorrelated heterogeneity gives rise to apparently radially axisymmetric concentration fronts which are macroscopically similar to those created by mechanical dispersion in homogeneous aquifers aside from some jaggedness at the circumference, it is to be expected that is small or close to zero under non-autocorrelated heterogeneity, and that increases as increases. The observation that the plume spreading arising from nonautocorrelated heterogeneity behaves similarly as from local mechanical dispersion has also been made previously (Luther and Haitjema, 1988).

Recovery Efficiency
The relationship between and are illustrated in Figures 4 and 5, and the fitted macrodispersion parameters ( 0, ) are listed in Table 1, where the reported 0 is the of the scenario with = 0.5 . Fitting of ( 0, ) in each individual scenario (e.g., Figure 4a for planar flow in a non-autocorrelated heterogeneous 2D domain) was accomplished by substituting = 0 2 0.5 2 into Equation 23 and minimizing the total RMSE of all four sets of analytical and numerical curves corresponding to all four values of . In other words, is minimized to fit ( 0, ), under the condition that = 0 2 0.5 2 . The following discussion, and the good agreement between the analytical and numerical outcomes, also apply to the relationship between and (not shown).

Figures 4 and 5 show that Equation 23
characterizes well the numerical recovery efficiency of an injection-extraction system, and its sensitivity to    Recall the approximate theoretical relationship ∝ 2 that was previously discussed. In non-autocorrelated aquifers, ∝ 2 is remarkably accurate (Figures 4a, 4b, and 5a). Figures 4 and 5 also shows that ∝ 2 is valid but sometimes slightly underpredicts (i.e., slightly overpredicts ) for stratified and exponentially autocorrelated heterogeneity, in agreement with Sommer et al. (2013) and Gelhar and Axness (1983). Despite the slight overprediction with ∝ 2 , the relative errors Figures 4 and 5 do not exceed 10% even though and 2 are varied over an order of magnitude.
In some cases, the recovery efficiency varies non-monotonically with injection duration, and attains a maximum. For a fixed value of , the value of at which the maximum point is located becomes smaller as increases, as can be seen in Figures 4e and 5c. While non-monotonicity in the recovery efficiency is only possible under 3D spherical flow in homogeneous aquifers (Equation 23; also see Tang & van der Zee, 2021), and only for ( ) but not ( ) , we show here that it is possible in 1D and 2D radial flow in heterogeneous aquifers for both ( ) and ( ) due to the positional scaling of macrodispersion.
Macrodispersion is weak compared to local dispersion at small plume travel distances (Molz et al., 1983), thus the recovery efficiency increases with (see Equation 11). In contrast, when macrodispersion dominates the spreading at large travel distances, the effects of local dispersion may be disregarded (Zhang & Neuman, 1996) and the recovery efficiency may decrease as increases for some combinations of . Depending on the values of , the dominance of macrodispersion over local dispersion increases with the distance from the inlet. Therefore, in scenarios where local dispersion dominates at small distances from the inlet, the recovery efficiency maximizes at some intermediate storage volume. The corresponding to maximum recovery can be found by equating the derivative of Equation 23 to zero.
The area-to-volume ratio A/V under the assumption of aquifer homogeneity, which is the area-to-volume ratio of a uniform d-dimensional sphere calculated as = ( ) 1 (Tang & van der Zee, 2021), is a popular metric for the estimation of the recovery efficiency of radial transport systems. Note that is not the actual area-to-volume ratio that takes into account heterogeneous flow paths. In homogeneous aquifers, an increase in the injected volume necessarily leads to a decrease in the actual A/V (which is identical to ) for geometrical reasons. Many studies therefore suggest that an increased injected volume leads to an increase in the recovery efficiency (e.g., Bloemendal & Hartog, 2018;Forkel & Daniels., 1995;Novo et al., 2010;Schout et al., 2014;Sommer et al., 2015). In aquifers with heterogeneous hydraulic conductivity and scale-dependent macrodispersion, this rule might no longer apply, as increasing would necessarily decrease but possibly lead to an increase in the actual A/V, due to the large interfaces that may develop between fast and slow flow zones. Accordingly, we found that increasing in beyond some threshold will decrease the recovery efficiency if 1 , even in scenarios where the recovery efficiency would have increased monotonically as in increased if the aquifer were to be homogeneous.

Effects of Heterogeneity Structure and Flow Field Geometry on Macrodispersion Parameters
The fitted values support the hypotheses presented in Section 3.4: is smallest under non-autocorrelated heterogeneity, intermediate under exponentially autocorrelated heterogeneity, and largest under layer stratification. Indeed, is smallest under non-autocorrelated heterogeneity for all tested scenarios. Preferential flow paths in stratified aquifers do not collide or fuse, therefore allowing preferential flow paths the greatest extent of persistence. Thus, the positional scaling of macrodispersion is larger in stratified than in heterogeneous aquifers with a finite autocorrelation length (Indelman & Dagan, 1999). Furthermore, the macrodispersion coefficient appears to converge toward some asymptotic value at large advective scales in practice (Zech et al., 2015), which agrees with our observation that is small or close to zero for non-autocorrelated heterogeneous aquifers, where the advective scale is much larger than the correlation scale of heterogeneity.
For planar flow in exponentially autocorrelated aquifers, if the storage radius is heavily increased such that , the system will eventually approach the limiting behavior of small associated with non-autocorrelated aquifers. Thus, is a decreasing function of when is neither infinite (stratified aquifer) nor zero (nonautocorrelated aquifer). However, at least within the parameter ranges investigated in this study, we observe the following. As seen in Figures 4c, 4d, and 5b even when is varied by one order of magnitude, modeling the recovery efficiency with a constant yields a good fit with the numerical simulations. This means that any variation of

Water Resources Research
TANG AND VAN DER ZEE 10.1029/2021WR030920 13 of 21 with is a higher-order effect that may be neglected in the design of aquifer storage systems. For stratified aquifers (i.e., infinite ) with flow parallel to the layers, the macrodispersion coefficient never converges to an asymptotic value (Matheron & De Marsily, 1980). Thus, grows unboundedly with in stratified aquifers, and never converges to zero.
The results listed in Table 1 agree with the mathematical analyses of Indelman and Dagan (1999). In all simulated heterogeneous aquifers, is larger under planar than cylindrical flow having otherwise identical circumstances. The scale-dependence of macrodispersion is larger for planar than for cylindrical flow in stratified and exponentially autocorrelated aquifers. Table 1 shows that for non-autocorrelated heterogeneity, is larger under cylindrical flow than planar flow. This agrees with Indelman (2004), who showed that when , as is the case for non-autocorrelated heterogeneity, the macrodispersion coefficient grows slowly with distance for cylindrical flow but not for planar flow, although we found that is small but not exactly zero under planar flow. In summary, for identical heterogeneous aquifers, imposing a different flow field results in different fitted macrodispersion parameters . Therefore, the macrodispersion parameters are emergent properties that arise from the dynamic interactions between flow field geometry and aquifer properties , . In contrast, the static properties , are intrinsic to the structure of aquifer heterogeneity. Hence, the presence of additional factors affecting the flow field, such as background flow, multiple interacting wells, or an unconfined upper boundary, can be expected to alter .

Validation of Fitted Macrodispersion Parameters
Laboratory experiments by Pang and Hunt (2001) suggest that is around 0.005. Review studies have shown that fitted to field tracer tests mostly range between 0.001 and 0.1 (Gelhar, 1986;Gelhar et al., 1992), though a more recent evaluation shows that clusters around 0.01 and below (Zech et al., 2015). The fitted to the simulation results (Table 1) range between 0.001 and 0.008 for = 0.5 , leading accordingly to a range between 0.016 and 0.128 for = 2 based on ∝ 2 . Therefore, the range of explored in this study covers the entire range of encountered in practice. Recall that = ( ) ∝ ′ . Table 1 shows that is different for different flow field geometries, and that for the simulated scenarios, lies within the range 1.04-1.4 (excluding non-autocorrelated aquifers). Zech et al.'s (2015) review, and Pickens and Grisak's (1981) numerical study, suggests that appears to vary across different geological formations and is thus non-universal. Various studies and reviews of field, laboratory, and numerical studies (Arya et al., 1988;Neuman, 1990;Schulze-Makuch, 2005) have found a range of 0.755-1.46 for in heterogeneous aquifers, which agrees well with the range of we found.
The review studies referenced in the previous paragraph obtained and by fitting a regression line between observed longitudinal dispersivity values and aquifer scale from an extensive database of numerical and experimental breakthrough curves. The agreement between the range of and we find and the range found in the cited studies indicate that the and we fit to recovery data, and the and those studies have obtained, correspond to the same physical process. Furthermore, our findings that ∝ 2 approximately, which as previously discussed agrees with the literature, also emphasize that the macrodispersion described by our model of recovery efficiency corresponds mathematically and physically to the macrodispersion phenomenon observed in breakthrough curves. These agreements thus validate the methodology that we have developed for predicting and interpreting the effects of aquifer heterogeneity and macrodispersion on recovery efficiency.
For the Cape Cod experiment, Fernàndez-Garcia et al. (2005) report 0.14 ≤ 2 ≤ 0.24 , = 1 , and 45 ≤ ′ ≤ 54 . Based on our findings for planar flow in an exponentially autocorrelated aquifer that = 0.004 when = 0.5 , and using the relationship ∝ 2 , we estimate that the macrodispersion coefficient ′ lies between 0.1 and 0.21 m. Field evidence suggests a value of 0.96 m for ′ (Garabedian et al., 1991). In contrast, Hess et al. (1992) estimated ′ in Cape Cod to be between 0.35 and 0.78 m based on the significantly more exact and complex methodology of Gelhar and Axness (1983). For the Borden experiment, also using data reported by Fernàndez-Garcia et al. (2005), we estimate ′ to be 0.18-0.42 m, whereas the field observed ′ is 0.45-0.49 m, and the estimate based on Gelhar and Axness' method is 0.61 m (Sudicky, 1986).
Although our model is somewhat less accurate for the Cape Cod experiment than for the Borden experiment, our model nevertheless provides reasonable estimates for both considering the similarly large margins of error of Gelhar and Axness' (1983) more exact and complex method. Other complex methods of estimating ′ , such 14 of 21 as those of Attinger et al. (2001), Gelhar (1993), and Chang and Yeh (2012), also result in similarly large or even larger margins of error, even when compared against numerical simulation data  that preclude the influence of measurement errors. We thus emphasize that our easily evaluable model provides utility for the characterization of macrodispersion despite being approximate.

Push-Pull Tests
Aquifer characterization by parameter inversion from single source injection-extraction cycles is known as the push-pull test method, and has thus far been successfully employed to determine a variety of aquifer properties (Ptak et al., 2004), such as the porosity (Hall et al., 1991), in situ microbial activity (Istok et al., 1997), fracture geometry (Klepikova et al., 2016), and geochemical reaction rates (Haggerty et al., 1998). The currently established methodology for obtaining field-scale data on dispersion and macrodispersion parameters is to measure breakthrough curves downstream of a solute source . For intentional solute introduction (with a cylindrical flow field), at least two wells are required, one for the injection and another downstream for the observations (e.g., Ptak & Teutsch, 1994). Rows of recharge wells can be lined up against rows of pumping wells to simulate forced planar flow (Molinari & Peaudeceff, 1977;Sauty, 1977). Furthermore, in heterogeneous aquifers, choosing suitable locations to place injection and observation wells is challenging due to possible bypass flow, and the scale-dependence of macrodispersion may require more wells to resolve. Gelhar and Collins' (1971) original model based on Equation 17, which does not account for macrodispersion, has been widely used to interpret the results of push-pull tests (e.g., Schroth & Istok, 2005;Schroth et al., 2000;Wang et al., 2018). For push-pull tests in aquifers with significant heterogeneity and macrodispersion, performing parameter estimation with such local-scale dispersion models will result in significant inaccuracy (Schroth et al., 2000). Our methods, which extend Gelhar and Collins' model to explicitly consider scale-dependent macrodispersion may therefore be more appropriate for interpreting the results of push-pull tests in heterogeneous aquifers.
Our findings that it is possible to fit these macrodispersion parameters with recovery data suggest the possibility of performing tracer tests with recovery data, in addition to solely breakthrough curve data. If Equation 23 is employed to estimate dispersion parameters, only a single well is required for cylindrical flow, and only a single row is required for planar flow, halving the number of required wells. While such "single-source" methods do not solve all difficulties related to aquifer characterization, such as an inability to resolve changes in the structure of heterogeneity across space, it can simplify the problem. If regional flow is negligible, the streamlines during the extraction phase are identical to those of the injection phase, but reversed. This ensures that the experiment is less sensitive to well placement, reducing the technical complexity and reducing the minimum aquifer size necessary to simulate a (semi-)infinite domain.
Though using breakthrough curves in push-pull tests imply more flexibility in characterizing observational data, in practice the shape of the breakthrough curve is often fixed a-priori by researchers (e.g., Davis et al., 2002;Schroth et al., 2000;Snodgrass & Kitanidis, 1998). With only dispersivity to be fitted for a designated flow geometry, breakthrough curves lose this advantage over recovery efficiency. While measuring the total mass of solutes in the recovered water to calculate the recovery efficiency is straightforward, the measurement of concentrations, especially flux-averaged concentrations (Bloem et al., 2008), to construct breakthrough curves is more complicated, prone to measurement errors (Moench, 1989;Novakowski, 1992;Palmer, 1988), and often involves uncertainty regarding field-appropriate boundary conditions (Ellsworth et al., 1996;Kreft & Zuber, 1978;Parker & van Genuchten, 1984). For push-pull tests with planar flow fields (Molinari & Peaudeceff, 1977;Sauty, 1977), using the recovery efficiency instead of breakthrough curves allows the challenging task of temporally synchronizing breakthrough curves recorded by a row of observation wells to be avoided. Breakthrough curve data is also vulnerable to random errors and limited measurement precision, whereas those would be smoothed out in recovery data, as it is essentially the cumulative breakthrough curve. This suggests that using the recovery efficiency data of single-source setups may be an attractive method of characterizing aquifers.
The main limitation of using the recovery efficiency in a push-pull test is that to resolve and , one has to obtain a function of or . Thus, multiple cycles with varying or are required, but the plumes of the second and subsequent cycles will interact with the residual solutes from the previous cycles, convoluting the data. Methods to deal with this disadvantage is a potential avenue for further research, although simply waiting between each cycle to allow the concentration or temperature to return to background levels is a possible solution. An alter-native is to take into account the effects of multiple consecutive cycles in the calculations when estimating the macrodispersion parameters. In addition, as the numerical scenarios simulated were meant to represent typical ATES and ASR applications, in which the storage radius may not be large enough to sample the aquifer heterogeneity under ergodic conditions, we emphasize that the macrodispersion parameters obtained from push-pull tests with the described method characterize only the local macrodispersion. The applicability of this method in obtaining ergodic macrodispersion parameters may perhaps be clarified through stochastic testing.

Multiple Storage Cycles
In scenarios with multiple cycles, the recovery efficiency of the -th cycle in homogeneous aquifers is given by (Tang & van der Zee, 2021) where 0 is the recovery efficiency of the first cycle, and is a constant with analytical bounds 1 2 < ≤ 1 . In this study, our numerical results show that this relation applies also to heterogeneous aquifers. Figure 6a shows that z increases as hydrodynamic dispersion increases, but does not depend on whether a 2D or 3D aquifer is simulated, assuming the conductivity distribution and flow field geometry are identical. Figure 6b shows that when comparing scenarios with similar 0 but different flow field geometries, z is not sensitive to the flow field geometry. Results show that z increases when increases if macrodispersion is dominant (i.e., large T) (Figure 6c), but not for small T where macrodispersion is weak (Figure 6d). Figures 6e and 6f show that z increases as the correlation length of heterogeneity increases (i.e., z is larger in stratified than in exponentially autocorrelated aquifers) for planar and cylindrical flow, respectively. Altogether, the various plots in Figure 6 show that Equation 26 excellently describes the recovery efficiency under multiple cycles.
Even when additional complexity such as dispersion to confining layers, multiple wells in a densely packed configuration, storage periods without flow, density-driven flow, and regional flow are present, the general form of Equation 26 agrees with experimental data (see Figure 7 of Bakr et al. (2013) and Figure 10b of Kastner et al., 2017), and numerical studies (see Figure 4 of Sommer et al., 2013 andFigure 5 of Zeghici et al., 2015). In particular, the heat or solute mass lost per cycle decreases rapidly, and the recovery efficiency increases sublinearly with an apparent power-law relation, as the number of cycles increases.
The exponent has to be empirically fitted and has a clear physical meaning, as is a measure of how much the recovery efficiency of a cycle is affected by all the preceding cycles. The solutes left behind in the aquifer after each cycle decreases the concentration gradient between injected and ambient water in subsequent cycles, thus decreasing mixing, and increasing the recovery efficiency with each new cycles. If = 1 , then the recovery efficiency of every cycle is identical, meaning that solutes dissipate away from the hydraulic front infinitely quickly. Hence, the larger is, the smaller this memory effect, and it is thus physically necessary that ≤ 1 . In the numerical results (Figure 6), increases as dispersion increases because in each subsequent cycle, the ambient chemical gradient enveloping the new injection water front dissipates more rapidly, which implies weakening the memory effect. It can also be observed in Figure 6 that in general, is larger when 0 is smaller, and vice-versa, though there does not appear to be a unique relationship between and 0 . In our simulations, highly efficient systems with 0.75 < 0 < 1 , which are of most interest for actual implementation of aquifer storage systems, correspond to 0.5 < < 0.65 . Additional numerical simulations (not shown) suggest that increases toward ≈ 0.8 as 0 ≈ 0.5 , then ≈ 1 as 0 ≈ 0 . This suggests the potential for future research to possibly uncover simple analytical or empirical relationships between 0 and . Ward et al. (2008) found that if stratified aquifers are simplified and modeled as homogeneous aquifers using the harmonic mean conductivity, the recovery efficiency of the simplified model converges to that of the actual stratified aquifer after 10 cycles. Figure 6 indeed reveals that after 10 cycles, the recovery efficiency of all the simulated scenarios converge to approximately 0.9 regardless of the heterogeneity structure and flow field geometry. After 10 cycles, the local dispersion parameters become more important than structural and geometrical factors in determining . This can be explained by the fact that the large surface area of solute plumes associated with stratified and autocorrelated heterogeneity become inconsequential after many cycles, as the ambient water's concentration between preferential flow channels gradually equilibrates with that of the injected water. Note, however, that Ward et al. simulated only one combination of well operational parameters. As we have shown, the 10.1029/2021WR030920 16 of 21 sensitivity of the recovery efficiency to may be directionally opposite between aquifers with and without scale-dependent macrodispersion. Hence, the simplification of stratified aquifers to a homogeneous one using a harmonic mean conductivity, or analogously the simplification of exponentially autocorrelated aquifers using a geometric mean conductivity (Zinn & Harvey, 2003), should not be applied to studies where the effect of varying the injection front position is to be investigated.

Limitations
A limitation of our methods is that the model of the local mechanical dispersion coefficient used in the numerical simulations is ( ) , but that used in the analytical solution is ⟨ ( )⟩ , where ⟨ ( )⟩ refers to the ensemble arithmetic mean velocity. When the spatial variance in flow velocity is small (weak heterogeneity), ⟨ ( )⟩ approximates ( ) , but the difference between the two models increases as the spatial variance in flow velocity increases (e.g., aquifers with strong heterogeneity and/or stratification). This may partly explain why the curve fitting under the imposed condition that ∝ 2 leads to the largest error for the stratified aquifers, because as increases the analytical model becomes insufficient to account for the effects of heterogeneity on local mechanical dispersion. This error appears to become large only for ≫ 1 (Boso et al., 2013;Jankovic et al., 2006). An analysis of the numerical results shows that the fit between the analytical and numerical curves, and the theoretical relationship ∝ 2 , are both stronger in stratified aquifers when scenarios with > 1 are excluded from the fitting. Nevertheless, the results of Jankovic et al. (2006) suggest that even in aquifers with ≫ 1 , applying an appropriate mean (i.e., fitted in the case of this study) and average velocity ( ) to Equation 23 will accurately account for over 95% of the transported solute mass.
Another limitation of our study is that we were unable to account for the magnitude of transverse dispersion in the mathematical analysis. Nevertheless, the recovery curves of numerical scenarios with both extremes of transverse dispersion magnitude: relatively small (nonautocorrelated heterogeneity) and relatively large (stratified aquifers), were successfully fitted to Equation 24. Furthermore, the extrapolated macrodispersion parameter values presented in Section 4.3 were reasonable. Therefore, we do not expect the omission of transverse dispersion from the mathematical analysis to invalidate our findings, at least under the common assumption that the transverse mechanical dispersivity is not larger than the longitudinal mechanical dispersivity. We reiterate that the purpose of this study is to argue that it is possible to fit the macrodispersivity parameters from recovery data using our methods, and not to find the parameter values corresponding to the heterogeneous conductivity distribution we simulated or to the field experiments discussed in this section using the most accurate possible method. In real-life cases, the intrinsic conductivity distribution and autocorrelation is seldom known, after all, thus it is arguably more important to understand the effects of macrodispersion on the recovery efficiency of aquifer storage systems, than to mathematically relate one abstract notion (conductivity distribution) to another (dispersion parameters).

Conclusion
The methods of Tang and van der Zee (2021) are extended to calculate the recovery efficiency of an injection-extraction well under scale-dependent dispersion, such as macrodispersion in a heterogeneous aquifer. An effective dispersion constant is constructed by adding a scale-dependent macrodispersion coefficient = to the local hydrodynamic dispersion coefficient ℎ = + . This general approach is able to characterize the recovery efficiency of aquifer storage systems and how they respond to variations in well operational parameters, number of recovery cycles, dispersion and macrodispersion parameters, aquifer heterogeneity structure, and flow field geometry.
Remarkably, if macrodispersion dominates and = 1 , the recovery efficiency becomes independent of both and . As increases, if 1 then increases, and if 1 then decreases. Hence, linear scaling marks an important threshold in the effects of macrodispersion on recovery. Since macrodispersion becomes stronger with displacement, whereas pore-scale dispersion is either scale-independent (molecular diffusion) or weakens with displacement (mechanical dispersion), the recovery efficiency may be a non-monotonic function of the operational parameters and in heterogeneous aquifers in 1D and 2D radial flow. In contrast, in homogeneous aquifers, nonmonotonicity of the recovery efficiency occurs only for ( ) under 3D radial flow, and never for ( ) . The macrodispersion parameters and may be fitted from recovery efficiency data, and are parameters that are not intrinsic to the hydrogeological structure and properties of the subsurface, as they depend also on the geometry of the flow field. The expression ∝ 2 , derived in several prior theoretical studies, agrees quite well with the values fitted to our numerical simulations. In aquifers where the heterogeneity has a longer autocorrelation range, such as perfectly stratified aquifers, tends to be larger.

10.1029/2021WR030920
18 of 21 Using the recovery efficiency as observational data in single-well tracer tests has several advantages over the classical downstream sampling of breakthrough curves. The necessary data are easier to measure and interpret, less vulnerable to erratic data caused by fluctuations in environmental conditions, less affected by bypass flow, necessitates fewer well installations, and can be performed in smaller aquifers. As far as we are aware, transport data from single well push-pull dispersivity tests have thus far been fitted only to local-scale mechanical dispersion models, and not to macrodispersion models in particular. Single well push-pull tests based on the recovery efficiency instead of the breakthrough curve have also not yet been attempted. The findings of this paper suggest that these two endeavors are possible, especially in combination, due to their mutual complementarity.
Our methods for interpreting single well push-pull tests in heterogeneous aquifers readily reveal key insights that are complex to uncover with classical methods involving breakthrough curve fitting. Under planar or cylindrical flow, if and only if the recovery efficiency is a non-monotonic function of the injection rate, duration, or volume, then the macrodispersion coefficient of the tested system increases superlinearly with distance from the well, implying the presence of well-connected high conductivity flow paths in the aquifer. The possible universality of macrodispersion scaling, as well as whether it is linear, sublinear, or superlinear with distance, remains a longstanding unresolved topic of ongoing research. While we did not fully resolve this question, we have introduced for the first time important practical implications for aquifer storage systems and push-pull tests that arise at the threshold separating sublinear, linear, and superlinear macrodispersion scaling.

Data Availability Statement
The MATLAB script used to solve the ADE to verify Equation 16 and Equation 24 is available at https://doi. org/10.5281/zenodo.5703348.