New model of reactive transport in a single-well push–pull test with aquitard effect and wellbore storage

The model of single-well push–pull (SWPP) test has been widely used to investigate reactive radial dispersion in remediation or parameter estimation of in situ aquifers. Previous analytical solutions only focused on a completely isolated aquifer for the SWPP test, excluding any influence of aquitards bounding the tested aquifer, and ignored the wellbore storage of the chaser and rest phases in the SWPP test. Such simplification 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 SWPP model is developed in an aquifer– aquitard system with wellbore storage, 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. As the permeability of the aquitard is much smaller than the permeability of the aquifer, the flow is assumed to be perpendicular to the aquitard; thus only vertical dispersive and advective transports are considered for the aquitard. The validity of this treatment is tested against results grounded in numerical simulations. The global sensitivity analysis indicates that the results of the SWPP test are largely sensitive (i.e., influenced by) to the parameters of porosity and radial dispersion of the aquifer, whereas the influence of the aquitard on results could not be ignored. In the injection phase, the larger radial dispersivity of the aquifer could result in the smaller values of breakthrough curves (BTCs), while there are greater BTC values in 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 is a generalization of several previous studies, and it performs better than previous studies ignoring the aquitard effect and wellbore storage for interpreting data of the field SWPP test reported by Yang et al. (2014).


Introduction
A single-well push-pull (SWPP) test could be applied for investigating aquifer properties related to reactive transport in the subsurface instead of the inter-well tracer test, due to its advantages of efficiency, low cost and easy implementation. The SWPP test is sometimes called the singlewell injection-withdrawal test, single-well huff-puff test or single-well injection-backflow test (Jung and Pruess, 2012). A complete SWPP 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 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 SWPP test is a forcedgradient 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 are available for various scenarios of the SWPP tests (Gelhar and Collins, 1971;Huang et al., 2010;Chen et al., 2017;Schroth and Istok, 2005;Wang et al., 2018). However, these studies are based on a common underlying assumption, i.e., that the studied aquifer was isolated from adjacent aquitards. In other words, the aquitards bounding the aquifer are taken as two completely impermeable barriers for solute transport. To date, numerous studies have 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 inevitably occurs when solute injected to the aquifer is close to the aquitardaquifer 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 SWPP test duration is relatively long, so the effect of molecular diffusion can materialize. Another important point to note is that the materials of the aquitard are usually clay and silt, which have strong absorbing capability for chemicals and great mass storage capacities (Chowdhury et al., 2017). To date, the influence of the aquitard on reactive transport in aquifers has attracted attention for several decades. As for radial dispersion, Chen (1985), Wang and Zhan (2013), and Zhou 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 SWPP test (injection).
Another assumption included in many previous models of radial dispersion is that the wellbore storage is ignored for the solute transport. In the injection phase of the SWPP test, the wellbore storage refers to the mixing processes between the prepared tracer injected into the wellbore and the original (or native) water in the wellbore. As a result of the wellbore storage, the concentration inside the wellbore varies with time until reaching the same value as the injected concentration, as shown in Fig. 1a. When ignoring it, the concentration inside the wellbore is constant during the entire inject phase, which is certainly not true. Similarly, the wellbore storage in the chaser, rest and extraction phases refers to the concentration variation caused by mixing processes between the original solute in the wellbore and the tracer moving in or out the wellbore. The examples of ignoring wellbore storage include Gelhar and Collins (1971), Chen (1985Chen ( , 1987, Moench (1989), Chen et al. (2007Chen et al. ( , 2012Chen et al. ( , 2017, Schroth et al. (2001), Tang and Babu (1979), Huang et al. (2010) and Zhou et al. (2017). Recently, Wang et al. (2018) developed a two-phase (injection and extraction) model for the SWPP test with specific considerations of the wellbore storage. In many field applications, the chaser and rest phases are generally involved and the mixing effect also happens in these two phases in the SWPP test, which will be investigated in this study.
Besides the abovementioned issues in previous studies, another issue is that the advection-dispersion equation (ADE) was used to govern the reactive transport of SWPP tests (Gelhar and Collins, 1971;Wang et al., 2018;Jung and Pruess, 2012). The validity of ADE has been challenged by numerous laboratory and field experimental studies before, when using a single representative value of advection, dispersion and reaction to characterize the whole system. In a hypothetical case, if great details of heterogeneity are known, one may employ a sufficiently fine mesh to discretize the porous media of concern and use ADE to capture anomalous transport characteristics fairly well (e.g., the early arrivals and/or heavy late-time tails of the breakthrough curves, BTCs). However, such a hypothetical case has rarely materialized in real applications, especially for field-scale problems. To remedy the situation (at least to some degree), the multi-rate mass transfer (MMT) model was proposed as an alternative to interpret the data of the SWPP test (Huang et al., 2010;Chen et al., 2017). In the MMT model, the porous media is divided into many overlapping continuums (Haggerty et al., 2000;Haggerty and Gorelick, 1995). A subset of MMT is the two overlapping continuums or the mobileimmobile model (MIM) in which the mass transfer between two domains (mobile and immobile) 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 advectiondispersion equation (FADE) model were also applied for anomalous reactive transport in SWPP tests (Hansen et al., 2017;Chen et al., 2017). Due to the complexity of the mathematic models of CTRW and FADE, it is very difficult or even impossible to derive analytical solutions for those two models, although both methods perform well in a numerical framework.
In this study, a new model of SWPP tests will be established by including both wellbore storage and the aquitard effect under the MIM framework. The reason for choosing MIM as the working framework is to capture the possible anomalous 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 SWPP test will be considered. The model of the wellbore storage will be developed using a mass balance principle in the chaser and rest phases. It does not seem difficult to solve this model of this study using numerical packages like MODFLOW-MT3DMS, TOUGH and TOUGHREACT, and FEFLOW. However, the numerical solutions may cause errors in treating the wellbore storage, since the volume of water in the wellbore was assumed to be constant (Wang et al., 2018), while in reality it changes with time and well discharge. Meanwhile, the numerical errors (like numerical dispersion and numerical oscillation) have to be considered in solving the ADE equation, especially for advection-dominated transport. In this study, an analytical solution will be derived to facilitate the data interpretation. Due to the format of analytical solutions, it is much easier to couple such solutions with a proper optimization algorithm (like a genetic algorithm). The analytical solution could serve as a benchmark to test the numerical solutions as well.

Model statement of the SWPP test
A single test well is assumed to fully penetrate an aquifer with uniform thickness. Both the aquifer and aquitards are homogeneous and extend laterally to infinity. Linear sorption and first-order degradation are included in the mathematic model of the SWPP test. Such assumptions might be oversimplified for cases in reality, while they are inevitable for the derivation of the analytical solution, especially for the aquifer homogeneity. For a heterogeneity aquifer, the solution presented here may be regarded as an ensemble-averaged approximation if the heterogeneity is spatially stationary. If the heterogeneity is spatially non-stationary, then one can apply non-stationary stochastic approach and/or Monte Carlo simulations to deal with the issue, which is out of the scope of this investigation.
The concept of homogeneity here deserves clarification. Despite the fact that the homogeneity assumption is commonly used in developing analytical and numerical 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. The Borden site of Canada (Sudicky, 1988) is one example of weak aquifer heterogeneity. Wang et al. (2018) employed a stochastic modeling technique to test the assumption of homogeneity associated with the SWPP test and found that such an assumption could be used to approximate a heterogeneous aquifer when the variance of spatial hydraulic conductivity was small. A cylindrical coordinate system is employed in this study, and the origin is located at the well center, as shown in Fig. 1c. The z axis and the r axis are vertical and horizontal, respectively. Figure 1 is a schematic diagram of the model investigated by this study.

Reactive transport model
Considering advective effect, dispersive effect and first-order chemical reaction in describing solute transport under the MIM framework, the governing equations for the SWPP test are

3986
Q. Wang et al.: New model of reactive transport in a single-well push-pull test with aquitard effect where subscripts "u" and "l" refer to parameters in the upper and lower aquitards, respectively; subscripts "m" and "im" refer to parameters in the mobile and immobile domains, respectively; C m and C im are the concentrations (ML −3 ) of the aquifer; C um , C uim , C lm and C lim are concentrations (ML −3 ) of the aquitards; t is the time (T ); B is half of the aquifer thickness (L); r is the radial distance (L); z represents the vertical distance (L); r w is the well radius (L); D r is aquifer dispersion coefficient (L 2 T −1 ); D u and D l are vertical dispersion coefficients (L 2 T −1 ) of the upper aquitard and lower aquitard, respectively; v a represents the average velocity (LT −1 ) in the aquifer and v a = u a θ m ; u a is Darcian velocity (LT −1 ); v um and v lm are vertical velocities (LT −1 ) in the aquitards; µ m , µ im , µ um , µ uim , µ lm and µ lim are reaction rates; θ m , θ im , θ um , θ uim , θ lm and θ lim are the porosities (dimensionless); θ lim are the retardation factors (dimensionless); K d is the equilibrium distribution coefficient (M −1 L 3 ); ρ b is the bulk density (ML −3 ); ω a , ω u and ω l are the first-order mass transfer coefficients (T −1 ).
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 where α r , α u and α l are dispersivities (L) of the aquifer, upper aquitard and lower aquitard, respectively; D * r , D * u and D * l are the diffusion coefficients (L 2 T −1 ). 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 where t inj , t cha , t res and t ext are the end moments (T ) of the injection phase, the chaser phase, the rest phase and the extraction phase, respectively; C inj, m (t), C cha, m (t), C res, m (t) and C ext, m (t) represent the wellbore concentrations (ML −3 ) of tracer in the injection phase, the chaser phase, the rest phase and the extraction phase, respectively. Eqs. (8)-(11) indicate that the flux continuity across the interface between well and the formation is only considered for the mobile continuum (or mobile domain).
The variation in the concentration with mixing effect in the injection phase could be described by Wang et al. (2018): where h w, inj is the wellbore water depth (L) in the injection phase, C 0 is concentration (ML −3 ) of prepared tracer. As for the chaser phase, the models describing the concentration variation in the wellbore could be obtained using mass balance principle: where h w, cha is the wellbore water depth (L) in the chaser phase.
In the extraction phase, the boundary condition is (Wang et al., 2018) V w, ext dC ext, m dt where h w, ext is the wellbore water depth (L) in the extraction phase.

Flow field model
The flow problem must be solved first before investigating the transport problem of the SWPP test. The velocity involved in the advection and dispersion terms of the governing Eqs. (1a) and (1b) is where Q 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 independent of time. This approximation is generally acceptable given the very limited spatial range of influence of most SWPP tests. For instance, if the characteristic length of SWPP test is l and the aquifer hydraulic diffusivity is D = K a /S a , where K a are S a are, respectively, the radial hydraulic conductivity and specific storage, then the typical characteristic time of unsteady-state flow is around t c ≈ l 2 2D . The typical characteristic time refers to the time of the flow changing from transient state to quasi-steady state, where the spatial distribution of flow velocity does not change while the drawdown varies with time. This model is similar to the model used to calculate the typical characteristic length of the tide-induced head fluctuation in a coastal aquifer system (Guarracino et al., 2012). For K a =1 m d −1 , S a = 10 −5 m −1 and l = 10 m (which are representative of an aquifer consisting of medium sands), one has t c ≈ l 2 2D = 5.0 × 10 −3 d, which is a very small value. To test the model in computing t c , the numerical simulation has been conducted, where the other parameters used in the model are the same as ones used in Figs. 2 and 3. Figure S2 shows the flow is in quasisteady state when time is greater than t c , since two curves of t = 5.0 × 10 −3 d and t = 10.0 × 10 −3 d overlap. As for the typical characteristic length, if the values of K a , S a and B have been estimated by the pumping tests before the SWPP test, it could be calculated by numerical modeling exercises using different simulation times.
The water levels in the wellbore in Eqs. (12)-(14) could be calculated by the models of Moench (1985): where p is Laplace transform variable; L −1 represents the inverse Laplace transform; the over bar represents the Laplacedomain variable, and where K u and K l are hydraulic conductivities (LT −1 ); S u and S l are specific storages (L −1 ); S w is the wellbore skin factor (dimensionless); B u and B l are thicknesses (L); and K 0 (·) and K 1 (·) are the modified Bessel functions.
3 New solution of reactive transport in the SWPP test In this study, the Laplace transform and Green's function methods will be employed to derive the analytical solution of the new SWPP test models described in Sect. 2. The dimensionless parameters are defined as follows: The detailed derivation of the new solution is listed in Sect. S1 in the Supplement.

Solutions in the Laplace domain
As for the injection phase of the SWPP test, the solutions in the Laplace domain are where s represents the Laplace transform parameter for t D (which is proportional to p); A i (·) is the Airy function; A i (·) is the derivative of the Airy function; and the expressions for a 2 , b 1 , E, y inj , y inj, w , ε m , ε im , ε um , ε uim , ε lm , ε lim , β inj and φ 1 are listed in Table 1.
In the chaser phase, the solutions of the SWPP test in the Laplace domain are where η varies between r wD and ∞ -e.g., r wD ≤ η ≤ ∞; η u varies between 1 and ∞; η l varies between −1 and −∞; C mD r D , t inj, D and C imD r D , t inj, D are the concentrations (ML −3 ) of the aquifer at the end of injection stage, which could be calculated by Eq. (25a) and (25b) after applying the inverse Laplace transform; C umD r D , z D , t inj, D and C uimD r D , z D , t cha, D represent the concentrations (ML −3 ) of the upper aquitard at the end of the injection phase, which could be calculated by Eq. (25c) and (25d) after applying the inverse Laplace transform; C lmD r D , z D , t inj, D and C limD r D , z D , t inj, D are the concentrations (ML −3 ) of the lower aquitard at the end of the injection phase, which could be calculated by Eq. (25e) and (25f) after applying the inverse Laplace transform; g (r D , E a ; η), g u (z D , E u ; η u ) and g l (z D , E l ; η l ) are Green's functions; and the expressions for g (r D , E a ; η), g u (z D , E u ; η u ), g l (z D , E l ; η l ), δ 1 , δ 2 , s 1 , s 2 , E a , E u , E l , y cha , y cha, w , F , ϕ (r D ), f u (η u ), X, M 1 , M 2 , M 3 , M 4 , N 1 , N 2 , N 3 , N 4 , T 1 , T 2 , T 3 , T 4 and β cha, D are listed in Table 2.
For the rest phase, the solutions of the SWPP test in the Laplace domain are where C mD r D , t cha, D and C imD r D , t cha, D are the concentrations (ML −3 ) of the aquifer at the end of the chaser phase, which could be calculated by Eq. (26a) and (26b)

3990
Q. Wang et al.: New model of reactive transport in a single-well push-pull test with aquitard effect where C mD r D , t res, D and C imD r D , t res, D are the concentrations (ML −3 ) of the aquifer at the end of the rest phase, which could be calculated by Eq. (27a) and (27b) after applying the inverse Laplace transform; C umD r D , z D , t res, D and C uimD r D , z D , t res, D are the concentrations (ML −3 ) of the upper aquitard at the end of the rest phase, which could be calculated by Eq. (27c) and (27d) after applying the inverse Laplace transform; C lmD r D , z D , t res, D and C limD r D , z D , t res, D are the concentrations (ML −3 ) of the lower aquitard at the end of the rest phase, which could be calculated by Eq. (27e) and (27f) after applying the inverse Laplace transform; b u varies between 1 and ∞; b l varies between −1 and −∞; ε varies between r wD and ∞ (e.g., r wD ≤ ε ≤ ∞); g (r D , ζ ; ε), g u (z D , E u ; b u ) and g l (z D , E l ; b l ) are Green's functions; and the expressions for , H 1 ∼ H 4 , I 1 ∼ I 4 , m 1 ∼ m 2 , n 1 ∼ n 2 , P 1 ∼ P 4 , W, y ext , y ext, w and β ext, D are listed in Table 3.

Solutions from the Laplace domain to the real-time domain
Because the analytical solutions in the Laplace domain are too complex, it seems impossible to transform it into the realtime domain analytically. Alternatively, a numerical method will be introduced for the inverse 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 SWPP test model
The new SWPP test model is a generalization of several previous studies; for instance, the new solution reduces to the solution of Gelhar and Collins (1971) when Chen et al. (2017) when ω u = ω l = D u = D l = v um = v lm = V w, inj = V w, cha = V w, ext = 0 and that of Wang et al. (2018) when ω a = ω u = ω l = D u = D l = v um = v lm = t cha = t res = 0. t cha = t res = 0 represents the four-phase SWPP test becoming the two-phase SWPP test, where the chaser and rest phases are excluded. Actually, all values of ω a , ω u , ω l , D u , D l , v um , v lm , V w, inj , V w, cha , and V w, ext are not zero in reality, which has been considered in the new solutions of this study. However, three assumptions still remain. First, the flow is in quasi-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 idea 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 validation of these assumptions will be discussed in the Sect. 4.2.

Verification of the new model
In this section, the newly derived analytical solutions will be tested from two aspects. Firstly, the new solution of this study could reduce to previous solutions under special cases, as the model established in this study is an extension of previous ones, and comparisons between them will be shown in Sect. 4.1. Secondly, although some assumptions included in previous models have been relaxed in the new model, some other processes of the reactive transport in the SWPP test have to be simplified in analytical solutions. Assumptions included in the new model have been discussed and their applicability is elaborated in Sect. 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; they ignored the aquitard effect and wellbore storage in the SWPP test. Figure 2 shows the comparison of BTCs between them, and the parameters used in such a comparison are R m = R im = R um = R uim = R lm = R lim =1, θ um = θ lm = 0.1, α r = α u = α l = 0.1 m, µ m = µ im = µ um = µ uim = µ lm = µ lim = 10 −6 d −1 , , V w, cha = 0 and V w, ext = 0 and imply that wellbore storage is neglected. The values of v um = v lm = 0 m d −1 mean that aquitards are neglected. As shown in Fig. 2, both solutions agree well for the mobile and immobile domains.

Test of assumptions involved in the analytical solution
To test three assumptions outlined in Sect. 3.3, a numerical model will be established, where general three-dimensional transient flow and solute transport are considered in both aquifer and aquitards. A finite-element method with the help of COMSOL Multiphysics will be used to solve the threedimensional model. The grid system is shown in Sect. S2. In this study, four sets of aquitard hydraulic conductivities are employed, i.e., K u = K l = 0.1 K a , K u = K l = 0.02 K a , K u = K l = 0.01 K a and K u = K l = 0.001 K a . A point to note is that the extreme case of K u = K l = 0.1 K a used here is only for the purpose of examining the robustness of comparison, while the real values of K u and K l are usually much lower than 0.1 K a . In other words, the remaining three cases mentioned above are more likely to occur in real applications.
As the first assumption in Sect. 3.3 has been elaborated on in Sect. S2.2, the following discussion will only focus on the second and third assumptions. Figure 3a, b and c represent the snapshots of concentration distributions in the aquifer along the r axis at different times. One may conclude that curves with smaller K u and K l values are closer to the analytical solution. This is because aquitards with smaller K u and K l (when K a remains constant) could make flow closer to the horizontal direction (or parallel with the aquitard-aquifer interface) in the aquifer and closer to the vertical direction (or perpendicular with the aquitard-aquifer interface) in the aquitard, according to the law of refraction (Fetter, 2018). In other words, when the values of K u /K a and K l /K a approach 0, the flow direction becomes horizontal in the aquifer and vertical in the aquitard, and then the numerical model reduces to the analytical model. Therefore, from this figure, one may conclude that the abovementioned second assumption in Sect. 3.3 works well in the aquifer when K u /K a and K l /K a are smaller than 0.01. Figure 4 shows the comparison of the analytical and numerical solutions for aquitards. Figure 4a1-c1 represent the snapshots of concentration distributions obtained from the analytical solution of this study at different times, and Fig. 4a2-c2 represent the snapshots of concentration distributions obtained from the numerical solution. One may find that the contour maps obtained from both solutions are Table 2. Expressions of the coefficients in the solutions expressed in Eq. (26a)-(26g). Figure 4. The vertical profiles (the r −z profiles) of the concentrations. Panels (a1), (b1) and (c1) represent the analytical solutions at t = 250, 300 and 500 d, respectively. Panels (a2), (b2) and (c2) represent the numerical solutions at t = 250, 300 and 500 d, respectively.
y cha, w r wD + 1 4E a almost the same in the aquifer, but very different in the aquitards. Therefore, the abovementioned third assumption in Sect. 3.3 is generally unacceptable in describing solute transport in the aquitard in the SWPP test but works well when the aquifer is of primary concern.

Model applications
As mentioned in Sect.3.3, the new model is a generalization of many previous models, and the conceptual model is closer 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 (α u and α l ), the first-order β ext, D exp(r wD /2)C mD( r wD , t res, D) sβ ext, D + 1 2 r wD −1− sβ ext, D + 1 2 r D | r D →∞ mass transfer coefficient (ω u and ω l ), the retardation factor (R um , R uim , R lm and R lim ), porosity (θ um , θ uim , θ lm and θ lim ), reaction rate (µ um , µ uim , µ lm and µ lim ) and velocity (v um and v lm ). The parameters involved for the aquifer include α r , ω a , R m , R im , θ m , θ im and B. Generally, these parameters could not be measured directly. Otherwise, they have to be obtained by fitting the experimental data using the forward model. Parameter estimation is an inverse problem, and it is generally conducted by an optimization model, such as a genetic algorithm and simulated annealing. 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 the literature as the initial guesses for similar lithology. Table 4 lists some parameter values for sandy and clay aquifers in previous studies. When a result is not sensitive to a particular parameter of concern, the value from previous publications for similar lithology and/or situations could be taken as an estimated value of that parameter if there is no direct measurement of that particular parameter of concern. To prioritize the sensitivity of predictions with respect to the diverse parameters involved in the new model, a global sensitivity analysis is conducted in Sect. 5.2.

A global sensitivity analysis
From the analytical solutions of Eqs. (26)-(28), one may find that BTCs are affected by several parameters, like α u , v um , θ um , ω, α r , θ m and V w . As α l , v lm and θ lm have a similar effect on the results as α u , v um and θ um , they have been excluded in the following analysis. In this section, a global sensitivity analysis is conducted using the model of Morris (1991), which is a one-step-at-a-time method. Morris (1991) employed µ k and σ k to represent the importance of the input parameters for the output concentration, and they have been computed by Morris (1991) and Lin et al. (2019): where M is the total sampling number, assuming that the range of parameter values is divided by M intervals; N is the total parameter number of interest, and it is 7 in this study; k is the kth parameter. In this study, M = 50; where P i is the random value of the ith parameter in the range of P i, 0 , P i, lim ; P i, 0 and P i, lim are the smallest and largest values of P i , as shown in Table S1; is a small increment defined as 1/(M − 1). A larger µ k means a higher sensitive effect of the kth parameter on the output, and a larger σ k represents the kth parameter having a greater interaction effect with others. Figure 5a and b represent the variation in µ k and σ k with time in the wellbore, respectively. The values of µ k are greater for α r and θ m than for the others, as shown in Fig. 5a, indicating that the influence of θ m and α r on the results is more obvious than others. However, the values of σ k are large for α u , θ um , α r , θ m and V w , demonstrating that the interactions of these parameters with others are strong; namely, the influence of them on results can also not be ignored.

Effect of the aquitard
As shown in Sect. 4.2, the new analytical solution is a good approximation for the numerical model in the aquifer when K u /K a and K l /K a are smaller than 0.01. In this section, we try to establish how the aquitards will affect BTCs of the SWPP tests. Since porosity is an important factor of concern, three sets of porosity values are used for the aquitards: θ um = θ lm = 0, 0.1 and 0.25. The other parameters are from the case in Fig. 4. Figure 6 shows the difference between the models with and without aquitards for different flow velocities in the aquitard. The case of θ um = θ lm = 0 represents the model without the aquitard. The difference is not obvious at the beginning of the extraction phase, while such a difference is obvious at the late time. Meanwhile, the smaller aquitard porosity makes the value of BTCs in the aquifer greater at a given time. When the aquitard is ignored, the values of BTCs are the greatest. 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  Brusseau et al. (1991). b Pickens et al. (1981). c Davis et al. (2003). d Javadi et al. (2017). e Liang et al. (2018). f Swami et al. (2016). g Kookana et al. (1992). h Haggerty et al. (1998). i Bouwer and McCarty (1985). j Chun et al. (2009). k Alvarez et al. (1991). References are shown in Sect. S3.  values will be used to analyze its influence: α r = 1.25, 2.50 and 5.00 m. Figure 7 shows BTCs in the well face for different radialdispersivity values. Firstly, the difference is obvious among curves in all phases. Secondly, a larger α r could decrease BTCs at a given time of the injection phase. This could be explained by the boundary condition of Eq. (8). The solute in the mobile domain of the aquifer is transported by both advection and dispersion; thus a larger α r could lower the values of C m in the well face. Thirdly, BTCs increase with increasing α r values in the chaser and rest phases. Fourthly, the peak values of BTCs decrease with increasing α r values.

Data interpretation: field SWPP 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 conducted in the borehole TW3 will be analyzed. The reason for choosing this dataset is that this borehole penetrated several layers, and it had been interpreted by Chen et al. (2017) before (using a model that does not consider the aquitard effect and the wellbore storage). The physical parameters of the SWPP test are r w = 0.1 m, Q inj = Q cha = 7.78 L min −1 , Q res = 0 L min −1 , Q ext = 12 L min −1 , t inj = 180 min, t cha = 26.74 min, t res = 10 080 min and B = 4 m. Other information on the experimental data can be found in Assayag et al. (2009) andYang et al. (2014). Figure 8a shows the fitness of the computed and observed BTCs. The estimated parameters are θ um = 0.05, θ lm = 0.0, θ m = 0.1, θ im = 0.068, α r = 0.5 m, α u = 0.35 m, α l = 0.0 m, R m = R im = R um = R uim = R lm = R lim = 1, µ m = µ im = µ um = µ uim = µ lm = µ lim = 10 −7 s −1 , ω = 0.001 d −1 , h w, inj = h w, cha = 32 m, h w, res = 30 m and h w, ext = 28 m. 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), where the error is defined as where C OBS and C COM are the observed and computed concentrations, respectively, and N is the number of sampling points. How accurate these parameters estimated by best fitting the observed data are in representing the real aquifer will be discussed in the following. The values of the retardation factor and reaction rate demonstrate that the chemical reaction and sorption are weak for the tracer of KBr in the SWPP test. It is not surprising since KBr is commonly treated as a "conservative" tracer. The porosity of the real aquifer ranges from 0.01 to 0.1, according to the well log analysis (Yang et al., 2014), where the estimated values are located. The estimated porosity represents the average values of the aquifer and aquitards. The estimated dispersivity of the aquifer is 0.7134 m by Chen et al. (2017), which is similar to ours. The values of the water level in the test could be observed directly; however, these data are not available, and they have to be estimated in this study. To evaluate the uncertainty in the estimated parameters, the sensitivity of the dispersivity to BTCs is analyzed, as shown in Fig. 8b. One may conclude that the estimated values of this study seem to be representative of reality, since the error is smallest for α r = 0.5 m.

Summary and conclusions
The single-well push-pull (SWPP) test could be applied to estimate the dispersivity, porosity and chemical reaction rates of in situ aquifers. However, previous studies mainly focused on an isolated aquifer, excluding all the possible effects of aquitards bounding the aquifer. In other words, 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 are derived to inspect the ef-fect of overlying and underlying aquitards. Meanwhile, four stages are considered in the new model with wellbore storage: 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.
To derive the analytical solution of the new model, some assumptions are inevitable. For instance, only vertical advection 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 modeling exercises are introduced to test the abovementioned assumptions of the new model. Based on this study, several conclusions could be obtained.
1. A new model of the SWPP test is a generalization of many previous models by considering the aquitard effect, the wellbore storage and the mass transfer rate in both aquifer and aquitards. A sub-model of wellbore storage is developed.
2. The assumption of vertical advection and dispersion in the aquitard and horizontal advection and dispersion in the aquifer are tested by specially designed finiteelement numerical models using COMSOL, and the result shows that this assumption is acceptable when the aquifer is of primary concern, provided that the ratios of the aquitard or 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 or aquifer permeability.
3. The new model is more sensitive to α r and θ m after a global sensitivity analysis, and the values of σ k are large for α u , θ um , α r , θ m and V w , demonstrating that the influence of the aquitard on results could not be ignored.
4. The performance of the new model is better than previous models that exclude the aquitard effect and the wellbore storage in terms of best fitting exercises with field data reported in Chen et al. (2017).
Data availability. All data are available in the Supplement.
Author contributions. QW wrote the paper and developed the mathematic models; JW and WS derived the analytical solutions; HZ revised the paper.