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).

1 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 single-well 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 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 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 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 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 (1985, 1987), Moench (1989), Chen et al. (2007, 2012, 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.

Figure 1 The schematic diagram of the SWPP test.

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 mobile–immobile 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 advection–dispersion 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.

2 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.

2.1 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

$$\begin{array}{}\text{(1a)}& {\displaystyle }\begin{array}{rl}& {\mathit{\theta }}_{\mathrm{m}}{R}_{\mathrm{m}}{\displaystyle \frac{\partial C_{\mathrm{m}}}{\partial t}}={\displaystyle \frac{{\mathit{\theta }}_{\mathrm{m}}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{D}_r{\displaystyle \frac{\partial C_{\mathrm{m}}}{\partial r}}\right)-{\mathit{\theta }}_{\mathrm{m}}{v}_{\mathrm{a}}{\displaystyle \frac{\partial C_{\mathrm{m}}}{\partial r}}\\ & -{\mathit{\omega }}_{\mathrm{a}}\left(C_{\mathrm{m}}-C_{\mathrm{im}}\right)-{\mathit{\theta }}_{\mathrm{m}}{\mathit{\mu }}_{\mathrm{m}}{C}_{\mathrm{m}}\\ & -{\left.\left({\displaystyle \frac{{\mathit{\theta }}_{\mathrm{um}}{v}_{\mathrm{um}}}{\mathrm{2}B}}{C}_{\mathrm{um}}-{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{um}}{D}_{\mathrm{u}}}{\mathrm{2}B}}{\displaystyle \frac{\partial C_{\mathrm{um}}}{\partial z}}\right)\right|}_{z=B}\\ & +{\left.\left(-{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{lm}}{v}_{\mathrm{lm}}}{\mathrm{2}B}}{C}_{\mathrm{lm}}-{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{lm}}{D}_{\mathrm{l}}}{\mathrm{2}B}}{\displaystyle \frac{\partial C_{\mathrm{lm}}}{\partial z}}\right)\right|}_{z=-B},r\ge r_{\mathrm{w}},\end{array}\text{(1b)}& {\displaystyle }{\mathit{\theta }}_{\mathrm{im}}{R}_{\mathrm{im}}{\displaystyle \frac{\partial C_{\mathrm{im}}}{\partial t}}={\mathit{\omega }}_{\mathrm{a}}\left(C_{\mathrm{m}}-C_{\mathrm{im}}\right)-{\mathit{\theta }}_{\mathrm{im}}{\mathit{\mu }}_{\mathrm{im}}{C}_{\mathrm{im}},r\ge r_{\mathrm{w}},\end{array}$$

$$\begin{array}{}\text{(2a)}& {\displaystyle }\begin{array}{rl}& {\mathit{\theta }}_{\mathrm{um}}{R}_{\mathrm{um}}{\displaystyle \frac{\partial C_{\mathrm{um}}}{\partial t}}={\mathit{\theta }}_{\mathrm{um}}{D}_{\mathrm{u}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{C}_{\mathrm{um}}}{\partial z^{\mathrm{2}}}}-{\mathit{\theta }}_{\mathrm{um}}{v}_{\mathrm{um}}{\displaystyle \frac{\partial C_{\mathrm{um}}}{\partial z}}\\ & -{\mathit{\omega }}_{\mathrm{u}}\left(C_{\mathrm{um}}-C_{\mathrm{uim}}\right)-{\mathit{\theta }}_{\mathrm{um}}{\mathit{\mu }}_{\mathrm{um}}{C}_{\mathrm{um}},z\ge B,\end{array}\text{(2b)}& {\displaystyle }\begin{array}{rl}& {\mathit{\theta }}_{\mathrm{uim}}{R}_{\mathrm{uim}}{\displaystyle \frac{\partial C_{\mathrm{uim}}}{\partial t}}={\mathit{\omega }}_{\mathrm{u}}\left(C_{\mathrm{um}}-C_{\mathrm{uim}}\right)-{\mathit{\theta }}_{\mathrm{uim}}{\mathit{\mu }}_{\mathrm{uim}}{C}_{\mathrm{uim}},\\ & z\ge B,\end{array}\end{array}$$

$$\begin{array}{}\text{(3a)}& {\displaystyle }\begin{array}{rl}& {\mathit{\theta }}_{\mathrm{lm}}{R}_{\mathrm{lm}}{\displaystyle \frac{\partial C_{\mathrm{lm}}}{\partial t}}={\mathit{\theta }}_{\mathrm{lm}}{D}_{\mathrm{l}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{C}_{\mathrm{lm}}}{\partial z^{\mathrm{2}}}}-{\mathit{\theta }}_{\mathrm{lm}}{v}_{\mathrm{lm}}{\displaystyle \frac{\partial C_{\mathrm{lm}}}{\partial z}}\\ & -{\mathit{\omega }}_{\mathrm{l}}\left(C_{\mathrm{lm}}-C_{\lim }\right)-{\mathit{\theta }}_{\mathrm{lm}}{\mathit{\mu }}_{\mathrm{lm}}{C}_{\mathrm{lm}},z\le -B,\end{array}\text{(3b)}& {\displaystyle }\begin{array}{rl}& {\mathit{\theta }}_{\lim }{R}_{\lim }{\displaystyle \frac{\partial C_{\lim }}{\partial t}}={\mathit{\omega }}_{\mathrm{l}}\left(C_{\mathrm{lm}}-C_{\lim }\right)-{\mathit{\theta }}_{\lim }{\mathit{\mu }}_{\lim }{C}_{\lim },\\ & z\le -B,\end{array}\end{array}$$

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⁻³) of the aquifer; C_um, C_uim, C_lm and C_lim are concentrations (ML⁻³) 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²T⁻¹); D_u and D_l are vertical dispersion coefficients (L²T⁻¹) of the upper aquitard and lower aquitard, respectively; v_a represents the average velocity (LT⁻¹) in the aquifer and $v_{\mathrm{a}}=\frac{u_{\mathrm{a}}}{{\mathit{\theta }}_{\mathrm{m}}}$; u_a is Darcian velocity (LT⁻¹); v_um and v_lm are vertical velocities (LT⁻¹) in the aquitards; μ_m, μ_im, μ_um, μ_uim, μ_lm and μ_lim are reaction rates; θ_m, θ_im, θ_um, θ_uim, θ_lm and θ_lim are the porosities (dimensionless); $R_{\mathrm{m}}=\mathrm{1}+\frac{{\mathit{\rho }}_{\mathrm{b}}{K}_{\mathrm{d}}}{{\mathit{\theta }}_{\mathrm{m}}}$, $R_{\mathrm{im}}=\mathrm{1}+\frac{{\mathit{\rho }}_{\mathrm{b}}{K}_{\mathrm{d}}}{{\mathit{\theta }}_{\mathrm{im}}}$, $R_{\mathrm{um}}=\mathrm{1}+\frac{{\mathit{\rho }}_{\mathrm{b}}{K}_{\mathrm{d}}}{{\mathit{\theta }}_{\mathrm{um}}}$, $R_{\mathrm{uim}}=\mathrm{1}+\frac{{\mathit{\rho }}_{\mathrm{b}}{K}_{\mathrm{d}}}{{\mathit{\theta }}_{\mathrm{uim}}}$, $R_{\mathrm{lm}}=\mathrm{1}+\frac{{\mathit{\rho }}_{\mathrm{b}}{K}_{\mathrm{d}}}{{\mathit{\theta }}_{\mathrm{lm}}}$ and $R_{\lim }=\mathrm{1}+\frac{{\mathit{\rho }}_{\mathrm{b}}{K}_{\mathrm{d}}}{{\mathit{\theta }}_{\lim }}$ are the retardation factors (dimensionless); K_d is the equilibrium distribution coefficient (M⁻¹L³); ρ_b is the bulk density (ML⁻³); ω_a, ω_u and ω_l are the first-order mass transfer coefficients (T⁻¹).

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

$$\begin{array}{}\text{(4a)}& {\displaystyle }{D}_r={\mathit{\alpha }}_r\left|v_r\right|+D_r^{*},\text{(4b)}& {\displaystyle }{D}_{\mathrm{u}}={\mathit{\alpha }}_{\mathrm{u}}\left|v_{\mathrm{u}}\right|+D_{\mathrm{u}}^{*},\text{(4c)}& {\displaystyle }{D}_{\mathrm{l}}={\mathit{\alpha }}_{\mathrm{l}}\left|v_{\mathrm{l}}\right|+D_{\mathrm{l}}^{*},\end{array}$$

where α_r, α_u and α_l are dispersivities (L) of the aquifer, upper aquitard and lower aquitard, respectively; $D_r^{*}$, $D_{\mathrm{u}}^{*}$ and $D_{\mathrm{l}}^{*}$ are the diffusion coefficients (L²T⁻¹).

Initial conditions are

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\left.C_{\mathrm{m}}\left(r,t\right)\right|}_{t=\mathrm{0}}={\left.C_{\mathrm{im}}\left(r,t\right)\right|}_{t=\mathrm{0}}={\left.C_{\mathrm{um}}\left(r,z,t\right)\right|}_{t=\mathrm{0}}\\ & ={\left.C_{\mathrm{uim}}\left(r,z,t\right)\right|}_{t=\mathrm{0}}={\left.C_{\mathrm{lm}}\left(r,z,t\right)\right|}_{t=\mathrm{0}}\\ & ={\left.C_{\lim }\left(r,z,t\right)\right|}_{t=\mathrm{0}}=\mathrm{0},r\ge r_{\mathrm{w}}.\end{array}\end{array}$$

The boundary conditions at infinity are

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\left.C_{\mathrm{m}}\left(r,t\right)\right|}_{r\to \mathrm{\infty }}={\left.C_{\mathrm{im}}\left(r,t\right)\right|}_{r\to \mathrm{\infty }}={\left.C_{\mathrm{um}}\left(r,z,t\right)\right|}_{z\to \mathrm{\infty }}\\ & ={\left.C_{\mathrm{uim}}\left(r,z,t\right)\right|}_{z\to \mathrm{\infty }}={\left.C_{\mathrm{lm}}\left(r,z,t\right)\right|}_{z\to -\mathrm{\infty }}\\ & ={\left.C_{\lim }\left(r,z,t\right)\right|}_{z\to -\mathrm{\infty }},=\mathrm{0},r\ge r_{\mathrm{w}}.\end{array}\end{array}$$

Due to the concentration continuity at the aquifer–aquitard interface, one has

$$\begin{array}{}\text{(7a)}& {\displaystyle }{C}_{\mathrm{m}}\left(r,t\right)=C_{\mathrm{um}}\left(r,z=B,t\right),\text{(7b)}& {\displaystyle }{C}_{\mathrm{m}}\left(r,t\right)=C_{\mathrm{lm}}\left(r,z=-B,t\right).\end{array}$$

The flux concentration continuity (FCC) is applied on the surface of wellbore, and one has

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}& {\left.\left[v_{\mathrm{a}}{C}_{\mathrm{m}}\left(r,t\right)-{\mathit{\alpha }}_r\left|v_{\mathrm{a}}\right|{\displaystyle \frac{\partial C_{\mathrm{m}}\left(r,t\right)}{\partial r}}\right]\right|}_{r=r_{\mathrm{w}}}\\ & ={\left.\left[v_{\mathrm{a}}{C}_{\mathrm{inj},\mathrm{m}}\left(t\right)\right]\right|}_{r=r_{\mathrm{w}}},\mathrm{0}<t\le t_{\mathrm{inj}},\end{array}\text{(9)}& {\displaystyle }\begin{array}{rl}& {\left.\left[v_{\mathrm{a}}{C}_{\mathrm{m}}\left(r,t\right)-{\mathit{\alpha }}_r\left|v_{\mathrm{a}}\right|{\displaystyle \frac{\partial C_{\mathrm{m}}\left(r,t\right)}{\partial r}}\right]\right|}_{r=r_{\mathrm{w}}}\\ & ={\left.\left[v_{\mathrm{a}}{C}_{\mathrm{cha},\mathrm{m}}\left(t\right)\right]\right|}_{r=r_{\mathrm{w}}},t_{\mathrm{inj}}<t\le t_{\mathrm{cha}},\end{array}\text{(10)}& {\displaystyle }{\left.\left[C_{\mathrm{m}}\left(r,t\right)\right]\right|}_{r=r_{\mathrm{w}}}=C_{\mathrm{res},\mathrm{m}}\left(r_{\mathrm{w}},t\right),t_{\mathrm{cha}}<t\le t_{\mathrm{res}},\text{(11)}& {\displaystyle }\begin{array}{rl}& {\left.\left[v_{\mathrm{a}}{C}_{\mathrm{m}}\left(r,t\right)-{\mathit{\alpha }}_r\left|v_{\mathrm{a}}\right|{\displaystyle \frac{\partial C_{\mathrm{m}}\left(r,t\right)}{\partial r}}\right]\right|}_{r=r_{\mathrm{w}}}\\ & ={\left.\left[v_{\mathrm{a}}{C}_{\mathrm{ext},\mathrm{m}}\left(t\right)\right]\right|}_{r=r_{\mathrm{w}}},t_{\mathrm{res}}<t\le t_{\mathrm{ext}},\end{array}\end{array}$$

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⁻³) 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):

$$\begin{array}{}\text{(12a)}& {\displaystyle }\begin{array}{rl}& V_{\mathrm{w},\mathrm{inj}}{\displaystyle \frac{\mathrm{d}{C}_{\mathrm{inj},\mathrm{m}}}{\mathrm{d}t}}=-\mathit{\xi }{v}_{\mathrm{a}}\left(r_{\mathrm{w}}\right)\left[C_{\mathrm{inj},\mathrm{m}}\left(t\right)-C_{\mathrm{0}}\right],\\ & \mathrm{0}<t\le t_{\mathrm{inj}},\end{array}\text{(12b)}& {\displaystyle }{\left.C_{\mathrm{inj},\mathrm{m}}\left(t\right)\right|}_{t=\mathrm{0}}=\mathrm{0},\mathrm{0}<t\le t_{\mathrm{inj}},\text{(12c)}& {\displaystyle }{V}_{\mathrm{w},\mathrm{inj}}=\mathit{\pi }{r}_{\mathrm{w}}^{\mathrm{2}}{h}_{\mathrm{w},\mathrm{inj}},\text{(12d)}& {\displaystyle }\mathit{\xi }=\mathrm{2}\mathit{\pi }{r}_{\mathrm{w}}{\mathit{\theta }}_{\mathrm{m}}\mathrm{2}B,\end{array}$$

where h_w, inj is the wellbore water depth (L) in the injection phase, C₀ is concentration (ML⁻³) of prepared tracer.

As for the chaser phase, the models describing the concentration variation in the wellbore could be obtained using mass balance principle:

$$\begin{array}{}\text{(13a)}& {\displaystyle }\begin{array}{rl}& V_{\mathrm{w},\mathrm{cha}}{\displaystyle \frac{\mathrm{d}{C}_{\mathrm{cha},\mathrm{m}}}{\mathrm{d}t}}=-\mathit{\xi }{v}_{\mathrm{a}}\left(r_{\mathrm{w}}\right)\left[C_{\mathrm{cha},\mathrm{m}}\left(t\right)\right],\\ & t_{\mathrm{inj}}<t\le t_{\mathrm{cha}},\end{array}\text{(13b)}& {\displaystyle }{\left.C_{\mathrm{cha},\mathrm{m}}\left(t\right)\right|}_{t=t_{\mathrm{inj}}}={\left.C_{\mathrm{inj},\mathrm{m}}\left(t\right)\right|}_{t=t_{\mathrm{inj}}},t_{\mathrm{inj}}<t\le t_{\mathrm{cha}},\text{(13c)}& {\displaystyle }{V}_{\mathrm{w},\mathrm{cha}}=\mathit{\pi }{r}_{\mathrm{w}}^{\mathrm{2}}{h}_{\mathrm{w},\mathrm{cha}},\end{array}$$

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)

$$\begin{array}{}\text{(14a)}& {\displaystyle }\begin{array}{rl}& V_{\mathrm{w},\mathrm{ext}}{\left.{\displaystyle \frac{\mathrm{d}{C}_{\mathrm{ext},\mathrm{m}}}{\mathrm{d}t}}\right|}_{r=r_{\mathrm{w}}}=-\mathit{\xi }{\mathit{\alpha }}_r{v}_{\mathrm{a}}\left(r_{\mathrm{w}}\right){\left.{\displaystyle \frac{\mathrm{d}{C}_{\mathrm{ext},\mathrm{m}}}{\mathrm{d}t}}\right|}_{r=r_{\mathrm{w}}},\\ & t_{\mathrm{res}}<\mathrm{t}\le t_{\mathrm{ext}},\end{array}\text{(14b)}& {\displaystyle }{\left.C_{\mathrm{ext},\mathrm{m}}\left(t\right)\right|}_{t=t_{\mathrm{res}}}={\left.C_{\mathrm{res},\mathrm{m}}\left(t\right)\right|}_{t=t_{\mathrm{res}}},t_{\mathrm{res}}<t\le t_{\mathrm{ext}},\text{(14c)}& {\displaystyle }{V}_{\mathrm{w},\mathrm{ext}}=\mathit{\pi }{r}_{\mathrm{w}}^{\mathrm{2}}{h}_{\mathrm{w},\mathrm{ext}},\end{array}$$

where h_w, ext is the wellbore water depth (L) in the extraction phase.

2.2 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

$$\begin{array}{}\text{(15)}& v_{\mathrm{a}}\left(r_{\mathrm{w}}\right)={\displaystyle \frac{Q}{\mathrm{4}\mathit{\pi }{r}_{\mathrm{w}}B{\mathit{\theta }}_{\mathrm{m}}}},r\ge r_{\mathrm{w}},\end{array}$$

where Q is the pumping rate (L³T⁻¹), 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_{\mathrm{a}}/S_{\mathrm{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_{\mathrm{c}}\approx \frac{l^{\mathrm{2}}}{\mathrm{2}D}$. 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⁻¹, $S_{\mathrm{a}}={\mathrm{10}}^{-\mathrm{5}}$ m⁻¹ and l=10 m (which are representative of an aquifer consisting of medium sands), one has $t_{\mathrm{c}}\approx \frac{l^{\mathrm{2}}}{\mathrm{2}D}=\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{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 quasi-steady state when time is greater than t_c, since two curves of $t=\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{3}}$ d and $t=\mathrm{10.0}\times {\mathrm{10}}^{-\mathrm{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.

Figure 2 Comparison of BTCs at the well screen computed by the solution of this study and Chen et al. (2017).

Figure 3 Comparison of the concentration distribution between the analytical and numerical solutions along the r axis at z=0 m. “ANA” and “NUM” represent the analytical and numerical solutions, respectively. (a) At the end of the injection phase: t=250 d. (b) At the end of the chasing phase: t=300 d. (c) In the extraction phase: t=500 d.

The water levels in the wellbore in Eqs. (12)–(14) could be calculated by the models of Moench (1985):

$$\begin{array}{}\text{(16)}& h_{\mathrm{w}}={lim}_{t\to \mathrm{\infty }}\left\{L^{-\mathrm{1}}\left[{\stackrel{\mathrm{‾}}{h}}_{\mathrm{w}}\left(p\right)\right]\right\},\end{array}$$

where p is Laplace transform variable; L⁻¹ represents the inverse Laplace transform; the over bar represents the Laplace-domain variable, and

$$\begin{array}{}\text{(17)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{h}}_{\mathrm{w}}\left(p\right)=h_{\mathrm{0}}\\ & -{\displaystyle \frac{Q}{\mathrm{8}\mathit{\pi }\text{KB}}}{\displaystyle \frac{\mathrm{2}\left[K_{\mathrm{0}}\left(x\right)+x{S}_{\mathrm{w}}{K}_{\mathrm{1}}\left(x\right)\right]}{p\left\{p{W}_{\mathrm{D}}\left[K_{\mathrm{0}}\left(x\right)+x{S}_{\mathrm{w}}{K}_{\mathrm{1}}\left(x\right)\right]+x{K}_{\mathrm{1}}\left(x\right)\right\}}},\end{array}\text{(18)}& {\displaystyle }{W}_{\mathrm{D}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}{\text{BS}}_{\mathrm{a}}}},\text{(19)}& {\displaystyle }x={\displaystyle \frac{\left(p+\stackrel{\mathrm{‾}}{q}\right)}{\mathrm{2}}},\text{(20)}& {\displaystyle }\stackrel{\mathrm{‾}}{q}={\left(\mathit{\gamma }{}^{\prime }\right)}^{\mathrm{2}}m{}^{\prime }\coth \left(m{}^{\prime }\right)+{\left(\mathit{\gamma }{}^{\prime \prime }\right)}^{\mathrm{2}}m{}^{\prime \prime }\coth \left(m{}^{\prime \prime }\right),\text{(21)}& {\displaystyle }m{}^{\prime }={\displaystyle \frac{{\left(\frac{S_{\mathrm{u}}{B}_{\mathrm{u}}p}{S_{\mathrm{a}}B}\right)}^{\mathrm{1}/\mathrm{2}}}{\mathit{\gamma }{}^{\prime }}},\text{(22)}& {\displaystyle }m{}^{\prime \prime }={\displaystyle \frac{{\left(\frac{S_{\mathrm{l}}{B}_{\mathrm{l}}p}{S_{\mathrm{a}}B}\right)}^{\mathrm{1}/\mathrm{2}}}{\mathit{\gamma }{}^{\prime \prime }}},\end{array}$$

$$\begin{array}{}\text{(23)}& {\displaystyle }\mathit{\gamma }{}^{\prime }=r_{\mathrm{w}}{\left({\displaystyle \frac{K_{\mathrm{u}}}{\mathrm{2}{K}_{\mathrm{a}}B{B}_{\mathrm{u}}}}\right)}^{\mathrm{1}/\mathrm{2}},\text{(24)}& {\displaystyle }\mathit{\gamma }{}^{\prime \prime }=r_{\mathrm{w}}{\left({\displaystyle \frac{K_{\mathrm{l}}}{\mathrm{2}{K}_{\mathrm{a}}B{B}_{\mathrm{l}}}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where K_u and K_l are hydraulic conductivities (LT⁻¹); S_u and S_l are specific storages (L⁻¹); S_w is the wellbore skin factor (dimensionless); B_u and B_l are thicknesses (L); and K₀(⋅) and K₁(⋅) 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: $C_{\mathrm{mD}}=\frac{C_{\mathrm{m}}}{C_{\mathrm{0}}}$, $C_{\mathrm{imD}}=\frac{C_{\mathrm{im}}}{C_{\mathrm{0}}}$, $C_{\mathrm{inj},\mathrm{mD}}=\frac{C_{\mathrm{inj},\mathrm{m}}}{C_{\mathrm{0}}}$, $C_{\mathrm{inj},\mathrm{imD}}=\frac{C_{\mathrm{inj},\mathrm{im}}}{C_{\mathrm{0}}}$, $C_{\mathrm{cha},\mathrm{mD}}=\frac{C_{\mathrm{cha},\mathrm{m}}}{C_{\mathrm{0}}}$, $C_{\mathrm{cha},\mathrm{imD}}=\frac{C_{\mathrm{cha},\mathrm{im}}}{C_{\mathrm{0}}}$, $C_{\mathrm{res},\mathrm{mD}}=\frac{C_{\mathrm{res},\mathrm{m}}}{C_{\mathrm{0}}}$, $C_{\mathrm{res},\mathrm{imD}}=\frac{C_{\mathrm{res},\mathrm{im}}}{C_{\mathrm{0}}}$, $C_{\mathrm{ext},\mathrm{mD}}=\frac{C_{\mathrm{ext},\mathrm{m}}}{C_{\mathrm{0}}}$, $C_{\mathrm{ext},\mathrm{imD}}=\frac{C_{\mathrm{ext},\mathrm{im}}}{C_{\mathrm{0}}}$, $C_{\mathrm{umD}}=\frac{C_{\mathrm{um}}}{C_{\mathrm{0}}}$, $C_{\mathrm{uimD}}=\frac{C_{\mathrm{uim}}}{C_{\mathrm{0}}}$, $C_{\mathrm{lmD}}=\frac{C_{\mathrm{lm}}}{C_{\mathrm{0}}}$, $C_{\mathrm{limD}}=\frac{C_{\lim }}{C_{\mathrm{0}}}$, $t_{\mathrm{D}}=\frac{\left|A\right|}{{\mathit{\alpha }}_r^{\mathrm{2}}{R}_{\mathrm{m}}}t$, $r_{\mathrm{D}}=\frac{r}{{\mathit{\alpha }}_r}$, $r_{\mathrm{wD}}=\frac{r_{\mathrm{w}}}{{\mathit{\alpha }}_r}$, $z_{\mathrm{D}}=\frac{z}{B}$, ${\mathit{\mu }}_{\mathrm{mD}}=\frac{{\mathit{\alpha }}_r^{\mathrm{2}}{\mathit{\mu }}_{\mathrm{m}}}{A}$, ${\mathit{\mu }}_{\mathrm{imD}}=\frac{{\mathit{\alpha }}_r^{\mathrm{2}}{R}_{\mathrm{m}}{\mathit{\mu }}_{\mathrm{im}}}{R_{\mathrm{im}}A}$, ${\mathit{\mu }}_{\mathrm{umD}}=\frac{{\mathit{\alpha }}_r^{\mathrm{2}}{\mathit{\mu }}_{\mathrm{um}}}{A}$, ${\mathit{\mu }}_{\mathrm{uimD}}=\frac{{\mathit{\alpha }}_r^{\mathrm{2}}{R}_{\mathrm{m}}{\mathit{\mu }}_{\mathrm{uim}}}{R_{\mathrm{im}}A}$, ${\mathit{\mu }}_{\mathrm{lmD}}=\frac{{\mathit{\alpha }}_r^{\mathrm{2}}{\mathit{\mu }}_{\mathrm{lm}}}{A}$, ${\mathit{\mu }}_{\mathrm{limD}}=\frac{{\mathit{\alpha }}_r^{\mathrm{2}}{R}_{\mathrm{m}}{\mathit{\mu }}_{\lim }}{R_{\mathrm{im}}A}$ and $A=\frac{Q}{\mathrm{4}\mathit{\pi }B{\mathit{\theta }}_{\mathrm{m}}}$. The detailed derivation of the new solution is listed in Sect. S1 in the Supplement.

3.1 Solutions in the Laplace domain

As for the injection phase of the SWPP test, the solutions in the Laplace domain are

$$\begin{array}{}\text{(25a)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}}\left(r_{\mathrm{D}},s\right)={\mathit{\varphi }}_{\mathrm{1}}\exp \left({\displaystyle \frac{y_{\mathrm{inj}}}{\mathrm{2}}}\right)A_i\left(E^{\mathrm{1}/\mathrm{3}}{y}_{\mathrm{inj}}\right),\\ & r_{\mathrm{D}}\ge r_{\mathrm{wD}},\end{array}\text{(25b)}& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{\mathrm{imD}}={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{im}}}{\left(s+{\mathit{\mu }}_{\mathrm{imD}}+{\mathit{\epsilon }}_{\mathrm{im}}\right)}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}},r_{\mathrm{D}}\ge r_{\mathrm{wD}},\text{(25c)}& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{\mathrm{umD}}={\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}}\exp (a_{\mathrm{2}}{z}_{\mathrm{D}}-a_{\mathrm{2}}),z_{\mathrm{D}}\ge \mathrm{1},\text{(25d)}& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{\mathrm{uimD}}={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{uim}}}{s+{\mathit{\epsilon }}_{\mathrm{uim}}+{\mathit{\mu }}_{\mathrm{uimD}}}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{umD}},z_{\mathrm{D}}\ge \mathrm{1},\text{(25e)}& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{\mathrm{lmD}}={\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}}\exp (b_{\mathrm{1}}{z}_{\mathrm{D}}+b_{\mathrm{1}}),z_{\mathrm{D}}\le -\mathrm{1},\text{(25f)}& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{\mathrm{limD}}={\displaystyle \frac{{\mathit{\epsilon }}_{\lim }}{s+{\mathit{\epsilon }}_{\lim }+{\mathit{\mu }}_{\mathrm{limD}}}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{lmD}},z_{\mathrm{D}}\le -\mathrm{1},\end{array}$$

where s represents the Laplace transform parameter for t_D (which is proportional to p); A_i(⋅) is the Airy function; $A_i{}^{\prime }(\cdot )$ is the derivative of the Airy function; and the expressions for a₂, b₁, E, y_inj, y_inj, w, ε_m, ε_im, ε_um, ε_uim, ε_lm, ε_lim, β_inj and ϕ₁ are listed in Table 1.

In the chaser phase, the solutions of the SWPP test in the Laplace domain are

$$\begin{array}{}\text{(26a)}& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}}=\mathrm{\Psi }\left(r_{\mathrm{D}}\right)+{\mathit{\delta }}_{\mathrm{1}}+{\mathit{\delta }}_{\mathrm{2}}{r}_{\mathrm{D}},r_{\mathrm{D}}\ge r_{\mathrm{wD}},\text{(26b)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{C}}_{\mathrm{imD}}={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{im}}}{\left(s+{\mathit{\mu }}_{\mathrm{imD}}+{\mathit{\epsilon }}_{\mathrm{im}}\right)}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}}+{\displaystyle \frac{C_{\mathrm{imD}}\left(r_{\mathrm{D}},t_{\mathrm{inj},\mathrm{D}}\right)}{\left(s+{\mathit{\mu }}_{\mathrm{imD}}+{\mathit{\epsilon }}_{\mathrm{im}}\right)}},\\ & r_{\mathrm{D}}\ge r_{\mathrm{wD}},\end{array}\text{(26c)}& {\displaystyle }\begin{array}{rl}& \mathrm{\Psi }\left(r_{\mathrm{D}},E_{\mathrm{a}};\mathit{\eta }\right)={\int }_{r_{\mathrm{wD}}}^{\mathrm{\infty }}g\left(r_{\mathrm{D}},E_{\mathrm{a}};\mathit{\eta }\right)\mathit{\phi }\left(\mathit{\eta }\right)\mathrm{d}\mathit{\eta },\\ & r_{\mathrm{D}}\ge r_{\mathrm{wD}},\end{array}\text{(26d)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{umD}}& ={\int }_{\mathrm{1}}^{\mathrm{\infty }}{g}_{\mathrm{u}}\left(z_{\mathrm{D}},E_{\mathrm{u}};{\mathit{\eta }}_{\mathrm{u}}\right)f_{\mathrm{u}}\left({\mathit{\eta }}_{\mathrm{u}}\right)\mathrm{d}{\mathit{\eta }}_{\mathrm{u}}\\ & +{\displaystyle \frac{z_{\mathrm{D}}-z_{e\mathrm{D}}}{\mathrm{1}-z_{e\mathrm{D}}}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}}\left(r_{\mathrm{D}},s\right),z_{\mathrm{D}}\ge \mathrm{1},\end{array}\text{(26e)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{uimD}}& ={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{uim}}}{s+{\mathit{\epsilon }}_{\mathrm{uim}}+{\mathit{\mu }}_{\mathrm{uimD}}}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{umD}}\\ & +{\displaystyle \frac{C_{\mathrm{uimD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{inj},\mathrm{D}}\right)}{s+{\mathit{\epsilon }}_{\mathrm{uim}}+{\mathit{\mu }}_{\mathrm{uimD}}}},z_{\mathrm{D}}\ge \mathrm{1},\end{array}\text{(26f)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{lmD}}& ={\int }_{-\mathrm{1}}^{-\mathrm{\infty }}{g}_{\mathrm{l}}\left(z_{\mathrm{D}},E_{\mathrm{l}};{\mathit{\eta }}_{\mathrm{l}}\right)f_{\mathrm{l}}\left({\mathit{\eta }}_{\mathrm{l}}\right)\mathrm{d}{\mathit{\eta }}_{\mathrm{l}}\\ & +{\displaystyle \frac{z_{e\mathrm{D}}+z_{\mathrm{D}}}{z_{e\mathrm{D}}-\mathrm{1}}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},s\right),z_{\mathrm{D}}\le -\mathrm{1},\end{array}\text{(26g)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{limD}}& ={\displaystyle \frac{{\mathit{\epsilon }}_{\lim }}{s+{\mathit{\epsilon }}_{\lim }+{\mathit{\mu }}_{\mathrm{limD}}}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{lmD}}\\ & +{\displaystyle \frac{C_{\mathrm{limD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{inj},\mathrm{D}}\right)}{s+{\mathit{\epsilon }}_{\lim }+{\mathit{\mu }}_{\mathrm{limD}}}},z_{\mathrm{D}}\le -\mathrm{1},\end{array}\end{array}$$

where η varies between r_wD and ∞ – e.g., $r_{\mathrm{wD}}\le \mathit{\eta }\le \mathrm{\infty }$; η_u varies between 1 and ∞; η_l varies between −1 and −∞; $C_{\mathrm{mD}}\left(r_{\mathrm{D}},t_{\mathrm{inj},\mathrm{D}}\right)$ and $C_{\mathrm{imD}}\left(r_{\mathrm{D}},t_{\mathrm{inj},\mathrm{D}}\right)$ are the concentrations (ML⁻³) 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_{\mathrm{umD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{inj},\mathrm{D}}\right)$ and $C_{\mathrm{uimD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)$ represent the concentrations (ML⁻³) 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_{\mathrm{lmD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{inj},\mathrm{D}}\right)$ and $C_{\mathrm{limD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{inj},\mathrm{D}}\right)$ are the concentrations (ML⁻³) 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\left(r_{\mathrm{D}},E_{\mathrm{a}};\mathit{\eta }\right)$,  $g_{\mathrm{u}}\left(z_{\mathrm{D}},E_{\mathrm{u}};{\mathit{\eta }}_{\mathrm{u}}\right)$ and $g_{\mathrm{l}}\left(z_{\mathrm{D}},E_{\mathrm{l}};{\mathit{\eta }}_{\mathrm{l}}\right)$ are Green's functions; and the expressions for $g\left(r_{\mathrm{D}},E_{\mathrm{a}};\mathit{\eta }\right)$, $g_{\mathrm{u}}\left(z_{\mathrm{D}},E_{\mathrm{u}};{\mathit{\eta }}_{\mathrm{u}}\right)$, $g_{\mathrm{l}}\left(z_{\mathrm{D}},E_{\mathrm{l}};{\mathit{\eta }}_{\mathrm{l}}\right)$, δ₁, δ₂, s₁, s₂, E_a, E_u, E_l, y_cha, y_cha, w, F, φ(r_D), f_u(η_u), X, M₁, M₂, M₃, M₄, N₁, N₂, N₃, N₄, T₁, T₂, T₃, T₄ and β_cha, D are listed in Table 2.

For the rest phase, the solutions of the SWPP test in the Laplace domain are

$$\begin{array}{}\text{(27a)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}}={\displaystyle \frac{C_{\mathrm{mD}}(r_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}})+\frac{{\mathit{\epsilon }}_{\mathrm{m}}{C}_{\mathrm{imD}}\left(r_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)}{\left(s+{\mathit{\mu }}_{\mathrm{imD}}+{\mathit{\epsilon }}_{\mathrm{im}}\right)}}{\left(s+{\mathit{\epsilon }}_{\mathrm{m}}+{\mathit{\mu }}_{\mathrm{mD}}-\frac{{\mathit{\epsilon }}_{\mathrm{m}}{\mathit{\epsilon }}_{\mathrm{im}}}{s+{\mathit{\mu }}_{\mathrm{imD}}+{\mathit{\epsilon }}_{\mathrm{im}}}\right)}},\\ & r_{\mathrm{D}}\ge r_{\mathrm{wD}},\end{array}\text{(27b)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{C}}_{\mathrm{imD}}={\displaystyle \frac{C_{\mathrm{imD}}\left(r_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)}{\left(s+{\mathit{\mu }}_{\mathrm{imD}}+{\mathit{\epsilon }}_{\mathrm{im}}\right)}}+{\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{im}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}}}{\left(s+{\mathit{\mu }}_{\mathrm{imD}}+{\mathit{\epsilon }}_{\mathrm{im}}\right)}},\\ & r_{\mathrm{D}}\ge r_{\mathrm{wD}},\end{array}\text{(27c)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{C}}_{\mathrm{umD}}={\displaystyle \frac{C_{\mathrm{umD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)+\frac{{\mathit{\epsilon }}_{\mathrm{um}}{C}_{\mathrm{uimD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)}{s+{\mathit{\epsilon }}_{\mathrm{uim}}+{\mathit{\mu }}_{\mathrm{uimD}}}}{\left(s+{\mathit{\epsilon }}_{\mathrm{um}}+{\mathit{\mu }}_{\mathrm{umD}}-\frac{{\mathit{\epsilon }}_{\mathrm{um}}{\mathit{\epsilon }}_{\mathrm{uim}}}{s+{\mathit{\epsilon }}_{\mathrm{uim}}+{\mathit{\mu }}_{\mathrm{uimD}}}\right)}},\\ & z_{\mathrm{D}}\ge \mathrm{1},\end{array}\text{(27d)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{uimD}}& ={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{uim}}}{s+{\mathit{\epsilon }}_{\mathrm{uim}}+{\mathit{\mu }}_{\mathrm{umD}}}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{umD}}\\ & +{\displaystyle \frac{C_{\mathrm{uimD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)}{s+{\mathit{\epsilon }}_{\mathrm{uim}}+{\mathit{\mu }}_{\mathrm{umD}}}},z_{\mathrm{D}}\ge \mathrm{1},\end{array}\text{(27e)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{C}}_{\mathrm{lmD}}={\displaystyle \frac{C_{\mathrm{lmD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)+\frac{{\mathit{\epsilon }}_{\mathrm{lm}}{C}_{\mathrm{limD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)}{s+{\mathit{\epsilon }}_{\lim }+{\mathit{\mu }}_{\mathrm{limD}}}}{\left(s+{\mathit{\epsilon }}_{\mathrm{lm}}+{\mathit{\mu }}_{\mathrm{lmD}}-\frac{{\mathit{\epsilon }}_{\mathrm{lm}}{\mathit{\epsilon }}_{\lim }}{s+{\mathit{\epsilon }}_{\lim }+{\mathit{\mu }}_{\mathrm{limD}}}\right)}},\\ & z_{\mathrm{D}}\le -\mathrm{1},\end{array}\text{(27f)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{limD}}& ={\displaystyle \frac{{\mathit{\epsilon }}_{\lim }}{s+{\mathit{\epsilon }}_{\lim }+{\mathit{\mu }}_{\mathrm{lmD}}}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{lmD}}\\ & +{\displaystyle \frac{C_{\mathrm{limD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)}{s+{\mathit{\epsilon }}_{\lim }+{\mathit{\mu }}_{\mathrm{lmD}}}},z_{\mathrm{D}}\le -\mathrm{1},\end{array}\end{array}$$

where $C_{\mathrm{mD}}\left(r_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)$ and $C_{\mathrm{imD}}\left(r_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)$ are the concentrations (ML⁻³) of the aquifer at the end of the chaser phase, which could be calculated by Eq. (26a) and (26b) after applying the inverse Laplace transform, $C_{\mathrm{umD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)$ and $C_{\mathrm{uimD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)$ are the concentrations (ML⁻³) of the upper aquitard at the end of the chaser phase, which could be computed by Eq. (26d) and (26e) after applying the inverse Laplace transform, and $C_{\mathrm{lmD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)$ and $C_{\mathrm{limD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{cha},\mathrm{D}}\right)$ are the concentrations (ML⁻³) of the lower aquitard at the end of the chaser phase, which could be calculated by Eq. (26f) and (26g) after applying the inverse Laplace transform.

As for the extraction phase of the SWPP test, the solutions in the Laplace domain are

$$\begin{array}{}\text{(28a)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}}\left(r_{\mathrm{D}},s\right)& =\exp \left(-r_{\mathrm{D}}/\mathrm{2}\right)\left[U\left(r_{\mathrm{D}},\mathit{\zeta };\mathit{\epsilon }\right)\right.\\ & \left.+{\mathit{\sigma }}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{2}}{r}_{\mathrm{D}}\right],r_{\mathrm{D}}\ge r_{\mathrm{wD}},\end{array}\text{(28b)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{imD}}& ={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{im}}}{\left(s+{\mathit{\mu }}_{\mathrm{imD}}+{\mathit{\epsilon }}_{\mathrm{im}}\right)}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}}\\ & +{\displaystyle \frac{C_{\mathrm{imD}}\left(r_{\mathrm{D}},\text{res}\right)}{s+{\mathit{\mu }}_{\mathrm{imD}}+{\mathit{\epsilon }}_{\mathrm{im}}}},r_{\mathrm{D}}\ge r_{\mathrm{wD}},\end{array}\text{(28c)}& {\displaystyle }U\left(r_{\mathrm{D}},\mathit{\zeta };\mathit{\epsilon }\right)={\int }_{r_{\mathrm{wD}}}^{\mathrm{\infty }}g\left(r_{\mathrm{D}},\mathit{\zeta };\mathit{\epsilon }\right)f\left(\mathit{\epsilon }\right)d\mathit{\epsilon },\text{(28d)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{umD}}& ={\int }_{\mathrm{1}}^{\mathrm{\infty }}{g}_{\mathrm{u}}\left(z_{\mathrm{D}},E_{\mathrm{u}};b_{\mathrm{u}}\right)f_{\mathrm{u}}\left(b_{\mathrm{u}}\right)\mathrm{d}{b}_{\mathrm{u}}\\ & +{\displaystyle \frac{z_{\mathrm{D}}-z_{e\mathrm{D}}}{\mathrm{1}-z_{e\mathrm{D}}}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}}\left(r_{\mathrm{D}},s\right),z_{\mathrm{D}}\ge \mathrm{1},\end{array}\text{(28e)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{uimD}}& ={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{uim}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{umD}}}{s+{\mathit{\epsilon }}_{\mathrm{uim}}+{\mathit{\mu }}_{\mathrm{uimD}}}}\\ & +{\displaystyle \frac{C_{\mathrm{uimD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{res},\mathrm{D}}\right)}{s+{\mathit{\epsilon }}_{\mathrm{uim}}+{\mathit{\mu }}_{\mathrm{uimD}}}},z_{\mathrm{D}}\ge \mathrm{1},\end{array}\text{(28f)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{lmD}}& ={\int }_{-\mathrm{1}}^{-\mathrm{\infty }}{g}_{\mathrm{l}}\left(z_{\mathrm{D}},E_{\mathrm{l}};b_{\mathrm{l}}\right)f_{\mathrm{l}}\left(b_{\mathrm{l}}\right)\mathrm{d}{b}_{\mathrm{l}}\\ & +{\displaystyle \frac{z_{\mathrm{D}}+z_{e\mathrm{D}}}{z_{e\mathrm{D}}-\mathrm{1}}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{mD}}\left(r_{\mathrm{D}},s\right),z_{\mathrm{D}}\le -\mathrm{1},\end{array}\text{(28g)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{limD}}& ={\displaystyle \frac{{\mathit{\epsilon }}_{\lim }{\stackrel{\mathrm{‾}}{C}}_{\mathrm{lmD}}}{s+{\mathit{\epsilon }}_{\lim }+{\mathit{\mu }}_{\mathrm{limD}}}}\\ & +{\displaystyle \frac{C_{\mathrm{limD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{res},\mathrm{D}}\right)}{s+{\mathit{\epsilon }}_{\lim }+{\mathit{\mu }}_{\mathrm{limD}}}},z_{\mathrm{D}}\le -\mathrm{1},\end{array}\end{array}$$

where $C_{\mathrm{mD}}\left(r_{\mathrm{D}},t_{\mathrm{res},\mathrm{D}}\right)$ and $C_{\mathrm{imD}}\left(r_{\mathrm{D}},t_{\mathrm{res},\mathrm{D}}\right)$ are the concentrations (ML⁻³) 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_{\mathrm{umD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{res},\mathrm{D}}\right)$ and $C_{\mathrm{uimD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{res},\mathrm{D}}\right)$ are the concentrations (ML⁻³) 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_{\mathrm{lmD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{res},\mathrm{D}}\right)$ and $C_{\mathrm{limD}}\left(r_{\mathrm{D}},z_{\mathrm{D}},t_{\mathrm{res},\mathrm{D}}\right)$ are the concentrations (ML⁻³) 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_{\mathrm{wD}}\le \mathit{\epsilon }\le \mathrm{\infty })$; $g\left(r_{\mathrm{D}},\mathit{\zeta };\mathit{\epsilon }\right)$, $g_{\mathrm{u}}\left(z_{\mathrm{D}},E_{\mathrm{u}};b_{\mathrm{u}}\right)$ and $g_{\mathrm{l}}\left(z_{\mathrm{D}},E_{\mathrm{l}};b_{\mathrm{l}}\right)$ are Green's functions; and the expressions for $g\left(r_{\mathrm{D}},\mathit{\zeta };\mathit{\epsilon }\right)$, $g_{\mathrm{u}}\left(z_{\mathrm{D}},E_{\mathrm{u}};b_{\mathrm{u}}\right)$, $g_{\mathrm{l}}\left(z_{\mathrm{D}},E_{\mathrm{l}};b_{\mathrm{l}}\right)$, σ₁, σ₂, Λ, ζ, f(ε), f_u(b_u), f_l(b_l), H₁∼H₄, I₁∼I₄, m₁∼m₂, n₁∼n₂, P₁∼P₄, W, y_ext, y_ext, w and β_ext, D are listed in Table 3.

Table 1 Expressions of the coefficients in the solutions expressed in Eq. (25a)–(25f).

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.

Table 2 Expressions of the coefficients in the solutions expressed in Eq. (26a)–(26g).

Table 3 Expressions of the coefficients in the solutions expressed in Eqs. (28a)–(28g).

3.2 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 real-time 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).

3.3 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 ${\mathit{\omega }}_{\mathrm{a}}={\mathit{\omega }}_{\mathrm{u}}={\mathit{\omega }}_{\mathrm{l}}=D_{\mathrm{u}}=D_{\mathrm{l}}=v_{\mathrm{um}}=v_{\mathrm{lm}}=V_{\mathrm{w},\mathrm{inj}}=V_{\mathrm{w},\mathrm{cha}}=V_{\mathrm{w},\mathrm{ext}}=t_{\mathrm{cha}}=t_{\mathrm{res}}=\mathrm{0}$, to the solution of Chen et al. (2017) when ${\mathit{\omega }}_{\mathrm{u}}={\mathit{\omega }}_{\mathrm{l}}=D_{\mathrm{u}}=D_{\mathrm{l}}=v_{\mathrm{um}}=v_{\mathrm{lm}}=V_{\mathrm{w},\mathrm{inj}}=V_{\mathrm{w},\mathrm{cha}}=V_{\mathrm{w},\mathrm{ext}}=\mathrm{0}$ and that of Wang et al. (2018) when ${\mathit{\omega }}_{\mathrm{a}}={\mathit{\omega }}_{\mathrm{u}}={\mathit{\omega }}_{\mathrm{l}}=D_{\mathrm{u}}=D_{\mathrm{l}}=v_{\mathrm{um}}=v_{\mathrm{lm}}=t_{\mathrm{cha}}=t_{\mathrm{res}}=\mathrm{0}$. $t_{\mathrm{cha}}=t_{\mathrm{res}}=\mathrm{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.

4 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.

4.1 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_{\mathrm{m}}=R_{\mathrm{im}}=R_{\mathrm{um}}=R_{\mathrm{uim}}=R_{\mathrm{lm}}=R_{\lim }=$1, ${\mathit{\theta }}_{\mathrm{um}}={\mathit{\theta }}_{\mathrm{lm}}=\mathrm{0.1}$, ${\mathit{\alpha }}_r={\mathit{\alpha }}_{\mathrm{u}}={\mathit{\alpha }}_{\mathrm{l}}=\mathrm{0.1}$ m, ${\mathit{\mu }}_{\mathrm{m}}={\mathit{\mu }}_{\mathrm{im}}={\mathit{\mu }}_{\mathrm{um}}={\mathit{\mu }}_{\mathrm{uim}}={\mathit{\mu }}_{\mathrm{lm}}={\mathit{\mu }}_{\lim }={\mathrm{10}}^{-\mathrm{6}}$ d⁻¹, r_w=0.2 m, Q_inj=2.5 m³ d⁻¹, Q_cha=2.5 m³ d⁻¹, Q_res=0 m³ d⁻¹, $Q_{\mathrm{ext}}=-\mathrm{2.5}$ m³ d⁻¹, t_inj=100 d, t_cha=50 d, t_res=40 d, B=5 m, θ_m=0.3, θ_im=0.15, ${\mathit{\theta }}_{\mathrm{uim}}={\mathit{\theta }}_{\lim }=\mathrm{0.1}$ and ω=0.001 d⁻¹. The parameters $h_{\mathrm{w},\mathrm{inj}}=h_{\mathrm{w},\mathrm{cha}}=h_{\mathrm{w},\mathrm{res}}=h_{\mathrm{w},\mathrm{ext}}=$0 represent V_w, inj =0, $V_{\mathrm{w},\mathrm{cha}}=\mathrm{0}$ and $V_{\mathrm{w},\mathrm{ext}}=\mathrm{0}$ and imply that wellbore storage is neglected. The values of $v_{\mathrm{um}}=v_{\mathrm{lm}}=\mathrm{0}$ m d⁻¹ mean that aquitards are neglected. As shown in Fig. 2, both solutions agree well for the mobile and immobile domains.

4.2 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 three-dimensional model. The grid system is shown in Sect. S2.

In this study, four sets of aquitard hydraulic conductivities are employed, i.e., $K_{\mathrm{u}}=K_{\mathrm{l}}=\mathrm{0.1}{K}_{\mathrm{a}}$, $K_{\mathrm{u}}=K_{\mathrm{l}}=\mathrm{0.02}{K}_{\mathrm{a}}$, $K_{\mathrm{u}}=K_{\mathrm{l}}=\mathrm{0.01}{K}_{\mathrm{a}}$ and $K_{\mathrm{u}}=K_{\mathrm{l}}=\mathrm{0.001}{K}_{\mathrm{a}}$. A point to note is that the extreme case of $K_{\mathrm{u}}=K_{\mathrm{l}}=\mathrm{0.1}{K}_{\mathrm{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.

The initial drawdown and the initial concentration are 0 for aquifer and aquitards. The hydraulic parameters are K_a=0.1 m d⁻¹ and $S_{\mathrm{a}}=S_{\mathrm{u}}=S_{\mathrm{l}}={\mathrm{10}}^{-\mathrm{4}}$ m⁻¹, and the other parameters are $R_{\mathrm{m}}=R_{\mathrm{im}}=R_{\mathrm{um}}=R_{\mathrm{uim}}=R_{\mathrm{lm}}=R_{\lim }=\mathrm{1}$, ${\mathit{\theta }}_{\mathrm{um}}={\mathit{\theta }}_{\mathrm{lm}}=\mathrm{0.1}$, α_r=2.5 m, ${\mathit{\alpha }}_{\mathrm{u}}={\mathit{\alpha }}_{\mathrm{l}}=\mathrm{0.5}$ m, ${\mathit{\mu }}_{\mathrm{m}}={\mathit{\mu }}_{\mathrm{im}}={\mathit{\mu }}_{\mathrm{um}}={\mathit{\mu }}_{\mathrm{uim}}={\mathit{\mu }}_{\mathrm{lm}}={\mathit{\mu }}_{\lim }={\mathrm{10}}^{-\mathrm{7}}$ s⁻¹, r_w=0.5 m, $Q_{\mathrm{inj}}=Q_{\mathrm{cha}}=\mathrm{50}$ m³ d⁻¹, Q_res=0 m³ d⁻¹, $Q_{\mathrm{ext}}=-\mathrm{50}$ m³ d⁻¹, t_inj=250 d, t_cha=50 d, t_res=50 d, B=10 m, θ_m=0.25, θ_im=0.05 and ω=0.01 d⁻¹. The comparison of concentration between the analytical and numerical solutions is shown in Figs. 3 and 4.

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 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.

5 Discussions

5.1 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 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.

Table 4 A partial list of parameters from the literature.

^a 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). ⁱ Bouwer and McCarty (1985). ^j Chun et al. (2009). ^k Alvarez et al. (1991). References are shown in Sect. S3.

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):

$$\begin{array}{}\text{(29a)}& {\displaystyle }{\mathit{\mu }}_k={\sum }_{l=\mathrm{1}}^M\left(\left|{\mathrm{EE}}_k^l\right|/M\right),k=\mathrm{1},\mathrm{2},\mathrm{\dots }N,\text{(29b)}& {\displaystyle }{\mathit{\sigma }}_k=\sqrt{{\displaystyle \frac{\mathrm{1}}{M}}{\sum }_{l=\mathrm{1}}^M{\left({\mathrm{EE}}_k^l-{\mathit{\mu }}_k\right)}^{\mathrm{2}}},k=\mathrm{1},\mathrm{2},\mathrm{\dots }N,\end{array}$$

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;

$${\mathrm{EE}}_k^l={\displaystyle \frac{{\displaystyle \begin{array}{c}{C}_{\mathrm{mD}}\left(P_{\mathrm{1}},P_{\mathrm{2}},\mathrm{\cdots },P_k+l\mathrm{\Delta },\mathrm{\cdots },P_N\right)\\ -C_{\mathrm{mD}}\left(P_{\mathrm{1}},P_{\mathrm{2}},\mathrm{\cdots },P_k,\mathrm{\cdots },P_N\right)\end{array}}}{l\mathrm{\Delta }}},$$

where P_i is the random value of the ith parameter in the range of $\left(P_{i,\mathrm{0}},P_{i,\lim }\right)$; 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 $\mathrm{1}/(M-\mathrm{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.

Figure 5 Sensitivity analysis. (a) Variation in μ_k with time in the wellbore. (b) Variation in σ_k with time in the wellbore.

5.3 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: ${\mathit{\theta }}_{\mathrm{um}}={\mathit{\theta }}_{\mathrm{lm}}=\mathrm{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 ${\mathit{\theta }}_{\mathrm{um}}={\mathit{\theta }}_{\mathrm{lm}}=\mathrm{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.

Figure 6 Comparison of BTCs between the model with and without aquitards for different porosities.

5.4 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 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 radial-dispersivity 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.

Figure 7 BTCs in the wellbore for different α_r.

6 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_{\mathrm{inj}}=Q_{\mathrm{cha}}=\mathrm{7.78}$ L min⁻¹, Q_res=0 L min⁻¹, Q_ext=12 L min⁻¹, 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) and Yang 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_{\mathrm{m}}=R_{\mathrm{im}}=R_{\mathrm{um}}=R_{\mathrm{uim}}=R_{\mathrm{lm}}=R_{\lim }=\mathrm{1}$, ${\mathit{\mu }}_{\mathrm{m}}={\mathit{\mu }}_{\mathrm{im}}={\mathit{\mu }}_{\mathrm{um}}={\mathit{\mu }}_{\mathrm{uim}}={\mathit{\mu }}_{\mathrm{lm}}={\mathit{\mu }}_{\lim }={\mathrm{10}}^{-\mathrm{7}}$ s⁻¹, ω=0.001 d⁻¹, $h_{\mathrm{w},\mathrm{inj}}=h_{\mathrm{w},\mathrm{cha}}=\mathrm{32}$ m, $h_{\mathrm{w},\mathrm{res}}=\mathrm{30}$ m and $h_{\mathrm{w},\mathrm{ext}}=\mathrm{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

$$\begin{array}{}\text{(30)}& \text{Error}={\sum }_{i=\mathrm{1}}^N{\left(C_{\text{OBS}}-C_{\text{COM}}\right)}^{\mathrm{2}},\end{array}$$

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.

Figure 8 Fitness of observed BTC. (a) Fitness of the observed data by different models. (b) Influence of the dispersivity of the aquifer on BTCs.

7 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 effect 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 finite-element 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).