New Model of Reactive Transport in Single-Well Injection- Withdrawal Test with Aquitard Effect

The model of single-well injection-withdrawal (SWIW) test has been widely used to investigate reactive radial dispersion in remediation or parameter estimation of the in situ aquifers. Previous analytical solutions only focused on a completely isolated aquifer for the SWIW test, excluding any influence of aquitards bounding the tested aquifer. This simplification 10 might be questionable in field applications when test durations are relatively long, because solute transport in or out of the bounding aquitards is inevitable due to molecular diffusion and cross-formational advective transport. Here, a new SWIW model is developed in an aquifer-aquitard system, and the analytical solution in the Laplace domain is derived. Four phases of the test are included: the injection phase, the chaser phase, the rest phase and the extraction phase. The Green‟s function method is employed for the solution in the extraction phase. As the permeability of aquitard is much smaller than the 15 permeability of the aquifer, the flow is assumed to be perpendicular to the aquitard, thus only vertical dispersive and advective transports are considered for aquitard. The validity of this treatment is tested by a numerical solution. The sensitivity analysis demonstrates that the influence of vertical flow velocity and porosity in the aquitards, and radial dispersion of the aquifer is more sensitive to the SWIW test than other parameters. In the injection phase, the larger radial dispersivity of the aquifer could result in the smaller values of breakthrough curves (BTCs), while greater values of BTCs of 20 the chaser and rest phases. In the extraction phase, it could lead to the smaller peak values of BTCs. The new model of this study performs better than previous studies excluding the aquitard effect for interpreting data of the field SWIW test.


Introduction 25
A single-well injection-withdrawal (SWIW) test could be applied for investigating aquifer properties related to reactive transport in subsurface instead of the inter-well tracer test, due to its advantages of efficiency, low cost, and easy implementation. The SWIW test is sometimes called the single-well push-pull test, or single-well huff-puff test, or singlewell injection-backflow test (Jung and Pruess, 2012). A complete SWIW test includes the injection, the chaser, the rest, and the extraction phase. The second and third phases are generally ignored in the analytical solutions, but recommended in the 30 field applications, since they could increase the reaction time for the injected chemicals with the porous media (Phanikumar and McGuire, 2010;Wang and Zhan, 2019).
Similar to other aquifer tests, the SWIW test is a forced-gradient groundwater tracer test, and analytical solutions are often preferred to determine the in situ aquifer properties, due to the computational efficiency. Currently, many analytical models were available for various scenarios of the SWIW tests (Gelhar and Collins, 1971;Huang et al., 2010;Chen et al., 2017;35 Schroth and Istok, 2005;Wang et al., 2018). However, these studies were based on a common underlying assumption, that the studied aquifer was isolated from adjacent aquitards. In another word, the aquitards bounding the aquifer are taken as two completely impermeable barriers for solute transport. To date, numerous studies demonstrated that such an assumption might cause errors for groundwater flow (Zlotnik and Zhan, 2005;Hantush, 1967), and for reactive transport (Zhan et al., 2009;Chowdhury et al., 2017;Li et al., 2019). This is because even without any flow in the aquitards, molecular diffusion is 40 inevitable to occur when solute injected to the aquifer is close to the aquitard-aquifer interface. This is particularly true when a fully penetrating well is used for injection, thus a portion of injected solute is very close to the aquitard-aquifer interface and the SWIW test duration is relatively long so the effect of molecular diffusion can be materialized. Another important point to note is that the materials of aquitard are usually clay and silt which have strong absorbing capability for chemicals and great mass storage capacities (Chowdhury et al., 2017). Actually, the influence of aquitard on reactive transport in 45 aquifers has attracted attentions for several decades. As for radial dispersion, Chen (1985), Zhan (2013) andZhou et al. (2017) presented analytical solutions for radial dispersion around an injection well in an aquifer-aquitard system. However, these models only focus on the first phase of the SWIW test (injection).
Another assumption included in many previous models of radial dispersion is that the concentration of the mixing water with the injected tracer is equal to the injected tracer concentration during the injection phase. The examples of employing such 50 an assumption include Gelhar and Collins (1971), Chen (1985Chen ( , 1987, Moench (1989), Chen et al. (2007Chen et al. ( , 2012, Schroth et al. (2001), Tang and Babu (1979), Chen et al. (2017), Huang et al. (2010), Chen et al. (2012), andZhou et al. (2017). This assumption implies that the mixing effect in the wellbore is not considered, where the mixing effect refers to the mixture between the original (or native) water and the injected tracer in the well. Such effect is excluded in almost all previous studies except Wang et al. (2018) for the SWIW test, who developed two-phase (injection and extraction) models with 55 specific considerations of the mixing effect. In many field applications, the chaser and rest phases are generally involved and the mixing effect also happens in these two phases in the SWIW test.
Besides above-mentioned issues in previous studies, another issue is that the advection-dispersion equation (ADE) was used to govern the reactive transport of SWIW tests (Gelhar and Collins,1971;2018;Jung and Pruess, 2012). The validity of ADE was challenged by numerous laboratory and field experimental studies in heterogeneous media, mostly 60 because ADE could not adequately interpret anomalous reactive transport, e.g. the early arrivals and/or heavy late-time tails of the breakthrough curves (BTCs). Alternatively, the multi-rate mass transfer (MMT) model was proposed to interpret the data of SWIW test (Huang et al., 2010;Chen et al., 2017). In the MMT model, the porous media is divided into two domains: One mobile domain where both dispersion and advection happen, and other immobile domain that only diffusion occurs (Haggerty et al., 2000;Haggerty and Gorelick, 1995). A subset of MMT is the mobile-immobile model (MIM) in which the 65 mass transfer between two domains becomes a single parameter instead of a function. The MIM model can grasp most characteristics of MMT and is mathematically simpler than MMT. Besides the MMT model, the continuous time random walk (CTRW) model and the fractional advection-dispersion equation (FADE) model were also applied for anomalous reactive transport in SWIW tests (Hansen et al., 2017;Chen et al., 2017). However, due to the complexity of the mathematic models of CTRW and FADE, it is very difficult, or even not possible to derive analytical solutions for those two models, 70 although both methods perform well in a numerical framework.
In this study, a new model of SWIW tests will be established by including both mixing effect in the wellbore and the aquitard effect under the MIM framework. The reason to choose MIM as the working framework is to capture the possible anonymous transport characteristics that cannot be described by ADE but at the same time to make the analytical treatment of the problem possible. Four stages of a SWIW test will be considered. The model of the mixing effect will be developed 75 using a mass balance principle in the chaser and rest phases. Analytical solution will be derived to facilitate the data interpretation for SWIW tests. The newly developed model will be checked against numerical solutions and field experimental data.

Model statement of the SWIW test
A single test well is assumed to fully penetrate an aquifer with uniform thickness. Both the aquifer and aquitards are 80 homogeneous and extend laterally to infinity. The concept of homogeneity here deserves some clarification. First, despite the fact that the homogeneity assumption is commonly used in developing analytical models of subsurface flow and transport, one should be aware that a rigorous sense of homogeneity probably never exists in a real-world setting (unless the media are composed of idealized glass balls as in some laboratory experiments). Therefore, the homogeneity concept here should be envisaged as a media whose hydraulic parameters vary within relatively narrow ranges, or the so-called weak heterogeneity. 85 Some examples of weak heterogeneity include the Borden Site of Canada (Sudicky, 1988). Wang et al. (2018) employed a stochastic modeling technique to test the assumption of homogeneity associated with the SWIW test, and found that such an assumption could be used to approximate a heterogeneous aquifer when the variance of spatial hydraulic conductivity was small. Second, for moderate or even strong heterogeneous media such as Cape Code site (Hess, 1989) or MADE site https://doi.org/10.5194/hess-2019-699 Preprint. Discussion started: 14 February 2020 c Author(s) 2020. CC BY 4.0 License. (Bohling et al., 2012), the analytical model developed under the homogeneity assumption is also valuable, but in a statistical 90 sense, as long as the media heterogeneity can be regarded as spatially stationary, meaning that the statistical structure of the media heterogeneity does not vary in space. In this setting, the analytical model developed under the homogeneity assumption is used to describe the (ensemble) average characteristics of an ensemble of heterogeneous media which are statistically identical but individually different. In another word, such an analytical model will provide a statistically average description of many realizations (an ensemble) which are similar to the heterogeneous media of concern, but it cannot 95 provide an exact description for the particular heterogeneous media under investigation. A cylindrical coordinate system is employed in this study, and the origin is located at the well center. The z-axis and the raxis are vertical and horizontal, respectively. A schematic diagram of the model investigated by this study is similar to

Reactive transport model 100
Considering advective effect, dispersive effect and first-order chemical reaction in describing solute transport under the MIM framework, the governing equations the SWIW test are: where subscript "" "" refers to parameters in the upper aquitard; subscript " " refers to parameters in the lower aquitard; 110 subscript "" "" refers to parameters in the mobile domain; subscript " " refers to parameters in the immobile domains; The symbol of the advection term is positive in the extraction phase in above equations, while it is negative before that. The dispersions are assumed to be linearly changing with the flow velocity, and one has: Initial conditions are: The boundary conditions at infinity are: Due to the concentration continuity at the aquifer-aquitard interface, one has: The flux concentration continuity (FCC) is applied on the surface of wellbore, and one has: 135 , where , , and are the end moments [T] of the injection phase, the chaser phase, the rest phase and the 140 extraction phase, respectively; ( ), ( ), ( ) and ( ) represent the wellbore concentrations [ML -3 ] of tracer in the injection phase, the chaser phase, the rest phase and the extraction phase, respectively.
The variation of the concentration with mixing effect in the injection phase could be described by (Wang et al., 2018): , where is the wellbore water depth [L] in the injection phase.
As for the chaser phase, the models describing the concentration variation in the wellbore could be obtained using mass balance principle: 150 , where is the wellbore water depth [L] in the chaser phase.
In the extraction phase, the boundary condition is (Wang et al., 2018): 155 , where is the wellbore water depth [L] in the extraction phase.

Flow field model 160
The flow problem must be solved first before investigating the transport problem of the SWIW test. The velocity involved in the advection and dispersion terms of the governing equations (1a) and (1b) is: where is the pumping rate [L 3 T -1 ], and it is negative for injection and positive for pumping. The use of Eq. (15) implies that quasi-steady state flow can be established very quickly near the injection/pumping well, thus the flow velocity becomes 165 independent of time. This approximation is generally acceptable given the very limited spatial range of influence of most SWIW tests. For instance, if the characteristic length of SWIW test is l and the aquifer hydraulic diffusivity is D=K a /S a , where K a are S a are the radial hydraulic conductivity and specific storage, then the typical characteristic time of unsteadystate flow is around . For instance, for a typical l=10 m, K a =10 m/day and S a =10 -5 (m -1 ) (which are representative of an aquifer consisting of medium sands), the value of t c is found to be day. 170 The water levels in the wellbore in Eqs. (12) -(14) could be calculated by the models of Moench (1985): where is Laplace transform variable; represents the inverse Laplace transform; the over bar represents the Laplacedomain variable, and . / , . / , .

New solution of reactive transport in the SWIW test 185
In this study, the Laplace transform and Green"s function methods will be employed to derive the analytical solution of the new SWIW test models described in Section 2. The dimensionless parameters are defined as follows: , , detailed derivation of the new solution is listed in Section S1 of Supplementary Materials.

Solutions in Laplace domain
As for the injection phase of the SWIW test, the solutions in Laplace domain are: where s represents the Laplace transform parameter for (which is proportional to p); ( ) is the Airy function ( ) is the 200 derivative of the Airy function; the expressions for , , , , listed in Table 1.
In the chaser phase, the solutions of the SWIW test in Laplace domain are:  Table 2.
For the rest phase, the solutions of the SWIW test in Laplace domain are: 220  As for the extraction phase of the SWIW test, the solutions in Laplace domain are: well mixing, and reducing the four-phase SWIW test into a two-phase SWIW test, the new model this study becomes Gelhar and Collins (1971).

Solutions from Laplace domain to real-time domain
Because the analytical solutions in Laplace domain are too complex, it seems impossible to transform it into the real time 255 domain analytically. Alternatively, a numerical method will be introduced for the invers Laplace transform. Currently, several methods are available, like the Stehfest model, Zakian model, Fourier series model, de Hoog model, and Schapery model (Wang and Zhan, 2015). Here, the de Hoog method will be applied to conduct the inverse Laplace transform, since it performed well for radial-dispersion problems (Wang et al., 2018;Wang and Zhan, 2013).

Assumptions included in the new SWIW test model 260
Although the new SWIW test model is a generalization of many previous studies, three assumptions still remain. First, the flow is in the steady state, e.g. Eq. (15). Second, the groundwater flow is horizontal in the aquifer, and is vertical in the aquitard. This treatment relies on the basis that the permeability of the aquitard is smaller than the permeability of the aquifer (Moench, 1985). Third, the model is simplified for the solute transport. For example, only vertical dispersion and advection effects are considered in the aquitard, and only radial dispersion and advection effects are considered in the aquifer. The 265 validation of these assumptions will be discussed in the Section 4.2.

Test of the new solution with previous solutions
To test the new solutions, the model of Chen et al. (2017) serves as a benchmark, who ignored the aquitard effect and wellbore mixing effect in the SWIW test. Figure 1 shows  =0" represent =0, =0 and =0, and imply that the mixing effect in the wellbore is neglected. The values of = =0 m/d mean that aquitards are neglected. As shown in Figure 1, both solutions agree well 275 for the mobile and immobile domains.

Test of assumptions involved in the analytical solution
To test the three assumptions outlined in Section 3.3, a numerical model will be established, where general three- As the first assumption in Section 3.3 has been elaborated in Section 2.2, the following discussion will only focus on the second and third assumptions. Therefore, the above-mentioned third assumption in Section 3.3 is generally unacceptable in describing solute transport in 305 the aquitard in the SWIW test, but works well when the aquifer is of the primary concern. 5 Discussions

Model applications
As mentioned in Section 3.1, the new model is a generalization of many previous models, and the conceptual model is more close to reality. However, there are many parameters involved in this new model that have to be determined first for applying this model. For instance, the involved parameters for the aquitards include dispersivity ( and ), first-order mass transfer coefficient ( and ), retardation factor ( , , , and ), porosity ( , , and ), reaction rate ( , , and ), and velocity ( and ). The involved parameters for the aquifer include , , , , , , and . Generally, these parameters could be measured directly. Otherwise, they could be obtained by fitting the experimental data using the forward model. 315 Parameter estimation is an inverse problem, and it is generally conducted by an optimization model, such as genetic algorithm, simulated annealing, and so on. Due to the ill-posedness of many inverse problems or insufficient observation data, the initial guess values of unknown parameters of interest are critical for finding the best values or real values of those parameters in the optimization model. Here, we recommend using values of parameters from literatures as the initial guesses for similar lithology. Table 4 lists some parameter values for sandy and clay aquifers in previous studies. When result is not 320 sensitive to a particular parameter of concern, the value from previous publications for similar lithology and/or situations could be taken as estimated value of that parameter, if there is no direct measurement of that particular parameter of concern.
To prioritize the sensitivity of parameters involved the new model, a sensitivity analysis is conducted in Section 5.2.

Sensitivity analysis
From the analytical solutions of Eqs. (26) -(28), one may find that BTCs are affected by several parameters, like , , 325 , , , and . In this section, the sensitivity analysis is conducted, and the model is (Kabala, 2001;Yang and Yeh, 2009): where is a small increment of I j ; is the sensitivity coefficient. could be any individual parameter of interest. A larger | | value represents that the result is more sensitive to that particular parameter. 330 Fig. 4 represents the sensitivity coefficients of BTCs. One may find that the influence of , , and on the results is more obvious than others. As the values of is genrouslly small, we mainly focus on the discussion of , and in the following sections.

Effect of the aquitard
As shown in Section 4.2, the new analytical solution is a good approximation for the numerical model in the aquifer when 335 and are smaller then 0.01. In this section, we try to figure out how the aquitards will affect BTCs of the SWIW tests. Since the porosity is an important factor of concern, three sets of porosity values are used for the aquitards: , 0.1, and 0.25. The other parameters are from the case in Figure 3. Therefore, the aquitard effect on transport in the aquifer is quite obvious and should not be ignored in general.

Effect of the aquifer radial dispersion
Another important parameter is the radial dispersion in the aquifer. In this section, three sets of the radial dispersivity values 345 will be used to analyze the influence: 1.25m, 2.50m, and 5.00m.

Data interpretation: Field SWIW test
To test the performance of the new model, the field data reported in Chen et al. (2017) will be employed. Specifically, the experimental data of S1 conducted in the borehole TW3 will be analysed. The reason choosing this dataset is because this borehole penetrated several layers, and it had been interpreted by Chen et al. (2017) before (using a model without 355 considering the aquitard effect and the mixing effect). Detailed information of experimental data could be seen in the references of Assayag et al. (2009) andYang et al. (2014).  Apparently, the fitness by the new solution is better than the model of Chen et al. (2017). As for the error between the observed and computed BTCs, the new solution is also smaller than that of Chen et al. (2017) as well, where the error is defined as where and are the observed and computed concentrations, respectively, and is the number of sampling points.

Summary and conclusions
The single-well injection-withdrawal (SWIW) test could be applied to estimate the dispersivity, porosity, chemical reaction rates of the in situ aquifers. However, previous studies mainly focused on an isolated aquifer, excluding all the possible effect of aquitards bounding the aquifer. In another word, the adjacent layers are assumed to be non-permeable, which is not exactly true in reality. In this study, a new analytical model is established and its associate solutions derived to inspect the effect of overlying and underlying aquitards. Meanwhile, four stages are considered in the new model, including the injection phase, the chaser phase, the rest phase and the extraction phase. The anomalous behaviors of reactive transport in the test were described by a mobile-immobile framework. The mixing effect was considered in the wellbore.
To derive the analytical solution of the new model, some assumptions are inevitable. For instance, only vertical advection 375 and dispersion are considered in the aquitard and only horizontal advection and dispersion are considered in the aquifer, and the flow is quasi-steady state. Although these assumptions have been widely used to describe the radial dispersion in previous studies, the influences on reactive transport have not been discussed in a rigorous sense before. In this study, numerical modelling exercises will be introduced to test the above-mentioned assumptions of the new model. Based on this study, the several conclusions could be obtained. 380 1. A new model of the SWIW test is a generalizing of many previous models by considering the aquitard effect, the wellbore mixing effect, and the mass transfer rate in both aquifer and aquitards. The sub-model of the wellbore mixing effect is developed.
2. Assumption of vertical advection and dispersion on the aquitard and horizontal advection and dispersion in the aquifer is tested by specially designed finite-element numerical models using COMSOL and the result shows that this assumption is 385 acceptable when the aquifer is of primary concern, provided that the ratios of the aquitard/aquifer permeability are less than 0.01; while such an assumption is generally unacceptable when the aquitards are of concern, regardless of the ratios of the aquitard/aquifer permeability.
3. The new model is most sensitive to the aquitard porosity and aquifer radial dispersivity after a comprehensive sensitivity analysis. A larger aquifer radial dispersivity decreases BTCs in the injection stage, increases BTCs in the chaser and rest 390 stages. It decreases BTC peak values in the extraction stage.