Reactive transport with wellbore storages in a single-well push–pull test

Using the single-well push–pull (SWPP) test to determine the in situ biogeochemical reaction kinetics, a chase phase and a rest phase were recommended to increase the duration of reaction, besides the injection and extraction phases. In this study, we presented multi-species reactive models of the four-phase SWPP test considering the wellbore storages for both groundwater flow and solute transport and a finite aquifer hydraulic diffusivity, which were ignored in previous studies. The models of the wellbore storage for solute transport were proposed based on the mass balance, and the sensitivity analysis and uniqueness analysis were employed to investigate the assumptions used in previous studies on the parameter estimation. The results showed that ignoring it might produce great errors in the SWPP test. In the injection and chase phases, the influence of the wellbore storage increased with the decreasing aquifer hydraulic diffusivity. The peak values of the breakthrough curves (BTCs) increased with the increasing aquifer hydraulic diffusivity in the extraction phase, and the arrival time of the peak value became shorter with a greater aquifer hydraulic diffusivity. Meanwhile, the Robin condition performed well at the rest phase only when the chase concentration was zero and the solute in the injection phase was completely flushed out of the borehole into the aquifer. The Danckwerts condition was better than the Robin condition even when the chase concentration was not zero. The reaction parameters could be determined by directly best fitting the observed data when the nonlinear reactions were described by piece-wise linear functions, while such an approach might not work if one attempted to use nonlinear functions to describe such nonlinear reactions. The field application demonstrated that the new model of this study performed well in interpreting BTCs of a SWPP test.

1 Introduction

A single-well push–pull (SWPP) test is a popular technique to characterize the in situ geological formations and to remedy the polluted aquifer by a series of biogeochemical reactions (Istok, 2012; Phanikumar and McGuire, 2010; Schroth and Istok, 2006). Therefore, the accuracy of the results is not only dependent on the experimental operation, but also on the conceptual model which is expected to properly represent the physical and biogeochemical processes. Unfortunately, most previous studies of the multi-species reactive transport models were based on some assumptions which may not be satisfied in actual applications, although those assumptions usually simplified the mathematical treatment of the problem (Istok, 2012; Wang et al., 2017).

As for the analytical solutions of the SWPP test, they have been widely used for applications, due to the high efficiency and great accuracy of the solutions, like the model of Gelhar and Collins (1971) for a fully penetrating well, the model of Schroth and Istok (2005) for a point source/sink well, and the model of Huang et al. (2010) for a partially penetrating well, assuming that the advection, the dispersion and the first-order reaction were involved in the transport processes. Haggerty et al. (1998) and Snodgrass and Kitanidis (1998) presented a simplified method based on a well-mixed reactor to estimate the first-order and zero-order reaction rate, without involving complex numerical modeling. Schroth and Istok (2006) provided two alternative models: one of them was a plug-flow model and the other was a variably mixed reactor model. Schroth et al. (2000) presented a simplified method for estimating retardation factors, based on the model of Gelhar and Collins (1971). Istok et al. (2001) extended the models of Haggerty et al. (1998) and Snodgrass and Kitanidis (1998) to estimate the Michaelis–Menten kinetic parameters which were used to describe the microbial respiration in the aquifer. Jung and Pruess (2012) presented a closed-form analytical solution for heat transport in a fractured aquifer involving a push-and-pull procedure. However, the above-mentioned analytical or semi-analytical solutions of the SWPP test were based on some over-simplified assumptions. For instance, the hydraulic diffusivity of the aquifer was assumed to be infinite, resulting in a time-independent flow velocity, where the hydraulic diffusivity is the ratio of the radial hydraulic conductivity over the specific storage. The wellbore storage effect on the flow field was assumed to be negligible as well. Therefore, how accurate parameter estimation could be needs to be tested. Recently, Wang et al. (2017) investigated the influences of a finite hydraulic diffusivity on the results and found that it might be significant, since both advective and dispersive transport were related to the flow velocity. One point to note is that the model of Wang et al. (2017) still contains an additional issue that has not been addressed: the wellbore storage influence on solute transport, which will be the focal point of this investigation.

The wellbore storage for solute transport refers to the variation of the solute injected into the wellbore during the processes of the test. A complete SWPP test contains four principle phases: injection of a prepared solution (tracer) into a targeted aquifer; injection of a chaser; rest period; extraction of the mixture solution. The second and third phases are optional, but are recommended to extend the reaction time of the tracer in the aquifer. In the injection phase, the concentration of the solute in the wellbore is smaller than that of the original solution at the early stage, since the original solute could be diluted by the original water in the wellbore, due to the mixing effect. Therefore, excluding the wellbore storage may overestimate the concentration in the wellbore at the early stage of the injection phase before the pre-test water inside the wellbore is completely flushed out of the borehole into the aquifer. In the chaser phase, the concentration of the solute in the wellbore may be greater than the concentration of the chaser, due to the mixing effect. The treatment of excluding the wellbore storage could underestimate the concentration in the wellbore at the early stage of the chase phase, due to the high concentration of solutes in the wellbore at the end of the injection phase. When the chaser phase is absent or the chaser concentration is not zero, the concentration might not be zero in the early stage of the rest phase. As for the chaser concentration, it is usually set to zero. However, under some circumstances, investigators may use a non-zero concentration for the chase phase. For example, Phanikumar and McGuire (2010) used 10 mg L⁻¹ for Cl⁻ and 2 mg L⁻¹ for ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in their chase solutions. Therefore, the concentration at the well screen may not be zero at the early stage of the rest phase when the chase concentration was not zero. All these mixing effects occurring in the wellbore are named the wellbore storage of the solute transport. Obviously, the assumption of ignoring the wellbore storage is not reasonable for the solute transport.

Actually, the above-mentioned assumptions used in the analytical and semi-analytical solutions can be relaxed in the numerical models, such as MODFLOW/MT3DMS (Harbaugh et al., 2000; Zheng and Wang, 1999), FEFLOW (Diersch, 2014), SUTRA (Voss, 1984), and STOMP (Nichols et al., 1997). Huang et al. (2010), Sun (2016), Haggerty et al. (1998), and Schroth and Istok (2006), respectively, employed such four software packages to carry out numerical simulations of SWPP tests, mainly involving advection, dispersion and first-order reaction. Unfortunately, the traditional three-dimensional models in the Cartesian coordinate system may create some errors in describing the wellbore storage of solute transport in the wellbore-confined aquifer, which is explained in the Supplement.

This study addresses multi-species reactive transport associated with SWPP tests with a better conceptual model that acknowledges the realistic circumstances that have been either overlooked or overly simplified in previous investigations. Firstly, we will employ a more realistic finite hydraulic diffusivity instead of an infinite hydraulic diffusivity to describe the flow field. Secondly, we will propose a better way to handle the boundary condition of transport at the wellbore by considering the wellbore storage effect for both groundwater and solute transport during the SWPP tests. Thirdly, the new model is tested using a field test dataset reported in McGuire et al. (2002). Fourthly, the sensitivity analysis and uniqueness analysis will be employed to investigate the assumptions used in previous studies on the parameter estimation.

2 Problem statement of the SWPP test

A cylindrical coordinate system is adopted with the r axis horizontal and the z axis vertically upward, as shown in Fig. 1. The origin is at the center of the well and is located in the plane of symmetry of the aquifer. The well fully penetrates a confined aquifer with a constant thickness. The aquifer is homogeneous, and the influence of the regional flow is ignored.

Figure 1 The schematic diagram of the SWPP test at the beginning of the rest phase when the chase concentration is not 0.

2.1 Revisit of the previous model

The general form of the governing equation for a multi-species reactive SWPP test is

$$\begin{array}{ll}{\displaystyle \frac{\partial C_i}{\partial t}}+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{b}}}{\mathit{\theta }}}{\displaystyle \frac{\partial S_i}{\partial t}}& {\displaystyle }=-\sum _{j=\mathrm{1}}^{N-\mathrm{1}}\left[H\left(t-t_j^{\ast }\right)-H\left(t-t_{j+\mathrm{1}}^{\ast }\right)\right]\\ \text{(1)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_j{C}_i^{n_j}\pm F_j,t>\mathrm{0},\end{array}$$

where C_i is the aqueous-phase concentration of the ith reactive solute, S_i is the solid-phase concentration of the ith reactive solute, t is the time in the SWPP test, ρ_b is the bulk density, θ is the porosity, H is the Heaviside step function, λ_j and n_j are the constant and orders, N is the number of the segment, $t_j^{\ast }$ and $t_{j+\mathrm{1}}^{\ast }$ are the times at two ends of segment j, and F_j is Monod/Michaelis–Menten kinetics. For the purpose of simplicity, we only present the reactive processes of the chemicals as described by Eq. (1), while the expressions of the transport (e.g., dispersion, diffusion, and advection) could be seen in Phanikumar and McGuire (2010), who used it to describe biogeochemical reactive transport of an arbitrary number of species including Monod/Michaelis–Menten kinetics, and the sorption models could be isotherm (Freundlich, Langmuir and linear sorption), one-site kinetic and two-site kinetic. As for the studies of Gelhar and Collins (1971), Schroth and Istok (2005), Huang et al. (2010), Haggerty et al. (1998), Snodgrass and Kitanidis (1998), Schroth and Istok (2006), Schroth et al. (2000), Istok et al. (2001), Jung and Pruess (2012), Wang et al. (2017), and so on, the governing equation is a special case of Eq. (1).

As mentioned in the Introduction, several assumptions may be debatable in previous studies and could be the source of errors for the actual applications. Firstly, the transport model is composed of a set of advection–dispersion equations (ADEs) built on the basis of flow velocity which is assumed to be time-independent (Chen et al., 2017; Gelhar and Collins, 1971; Huang et al., 2010; Phanikumar and McGuire, 2010):

$$\begin{array}{}\text{(2)}& v_{\mathrm{r}}={\displaystyle \frac{Q}{\mathrm{2}\mathit{\pi }rB\mathit{\theta }}},r\ge r_{\mathrm{w}},\end{array}$$

where r_w is the well radius; r is the radius distance from the center of the well; B is the aquifer thickness; Q is the flow rate of the well; $v_{\mathrm{r}}=u_{\mathrm{r}}/\mathit{\theta }$ is the average radial pore velocity and u_r is the radial Darcian velocities. Equation (2) implies that the hydraulic diffusivity of the aquifer is infinite; thus, the flow velocity is independent of time. Meanwhile, the wellbore storage is negligible or the well radius r_w is assumed to be infinitesimal in formulating Eq. (2).

The second assumption of the model is the boundary condition of the well screen in the rest phase of the SWPP test, in which a Robin condition (or a third-type condition) is employed to describe the aqueous solute transport (Chen et al., 2017; Phanikumar and McGuire, 2010; Wang et al., 2017):

$$\begin{array}{}\text{(3)}& {\left.\left(v_{\mathrm{r}}C-D_{\mathrm{r}}{\displaystyle \frac{\partial C}{\partial r}}\right)\right|}_{r\to r_{\mathrm{w}}}=\mathrm{0},t_{\mathrm{inj}}+t_{\mathrm{cha}}<t\le t_{\mathrm{inj}}+t_{\mathrm{cha}}+t_{\mathrm{res}},\end{array}$$

where t_inj, t_cha, t_res, and t_ext represent the durations of the injection, chase, rest, and extraction phases, respectively; C is the resident concentration of the aqueous phase to represent C_i in Eq. (1); D_r is the dispersion coefficient, which is

$$\begin{array}{}\text{(4)}& D_{\mathrm{r}}={\mathit{\alpha }}_{\mathrm{r}}{v}_{\mathrm{r}}+D_{\mathrm{0}},\end{array}$$

in which α_r is the radial dispersivity; D₀ is the effective diffusion coefficient in the aquifer.

Thirdly, a constant solute concentration in the wellbore is applied in the injection and chase phases without considering the solute diluted effect in the wellbore (Chen et al., 2017; Gelhar and Collins, 1971; Istok, 2012; Phanikumar and McGuire, 2010; Wang et al., 2017):

$$\begin{array}{}\text{(5a)}& {\displaystyle }& {\displaystyle }{\left.\left(v_{\mathrm{r}}C-D_{\mathrm{r}}{\displaystyle \frac{\partial C}{\partial r}}\right)\right|}_{r\to r_{\mathrm{w}}}=v_{\mathrm{r}}{C}_{\mathrm{0}}^{\mathrm{inj}},\mathrm{0}<t\le t_{\mathrm{inj}},\text{(5b)}& {\displaystyle }& {\displaystyle }\text{or}{\left.C\right|}_{r\to r_{\mathrm{w}}}=C_{\mathrm{0}}^{\mathrm{inj}},\mathrm{0}<t\le t_{\mathrm{inj}},\end{array}$$

$$\begin{array}{}\text{(6a)}& {\displaystyle }& {\displaystyle }{\left.\left(v_{\mathrm{r}}C-D_{\mathrm{r}}{\displaystyle \frac{\partial C}{\partial r}}\right)\right|}_{r\to r_{\mathrm{w}}}=v_{\mathrm{r}}{C}_{\mathrm{0}}^{\mathrm{cha}},t_{\mathrm{inj}}<t\le t_{\mathrm{inj}}+t_{\mathrm{cha}},\text{(6b)}& {\displaystyle }& {\displaystyle }\text{or}{\left.C\right|}_{r\to r_{\mathrm{w}}}=C_{\mathrm{0}}^{\mathrm{cha}},t_{\mathrm{inj}}<t\le t_{\mathrm{inj}}+t_{\mathrm{cha}},\end{array}$$

where $C_{\mathrm{0}}^{\mathrm{inj}}$ and $C_{\mathrm{0}}^{\mathrm{cha}}$ represent the solute concentrations injected into the wellbore during the injection and chase phases, respectively. A detailed discussion about the above-mentioned assumptions can be seen in Phanikumar and McGuire (2010) or Wang et al. (2017).

Fourthly, the solute transport caused by dispersion and advection was assumed to be negligible in estimating the reaction rates. For instance, one of the simplest models of such reactions may be the first-order reaction

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial C}{\partial t}}=-\mathit{\lambda }C,\end{array}$$

where λ is the first-order reaction rate constant. Besides the first-order reaction, Eq. (7) could be used to describe the first-order biodegradation and radioactive decay. Haggerty et al. (1998) presented a simplified method to estimate λ for the SWPP test:

$$\begin{array}{}\text{(8)}& \ln \left({\displaystyle \frac{C_{\mathrm{rec}}\left(t^{\ast }\right)}{C_{\mathrm{tra}}\left(t^{\ast }\right)}}\right)=\ln \left({\displaystyle \frac{\mathrm{1}-\exp \left(-\mathit{\lambda }{t}_{\mathrm{inj}}\right)}{\mathit{\lambda }{t}_{\mathrm{inj}}}}\right),\end{array}$$

where t^∗ is the time since the end of injection; C_rec(t^∗) is the reactant concentration; C_tra(t^∗) is the concentration of a conservative tracer. To obtain the value of λ, the reactant and the conservative tracer should be fully mixed and injected into the aquifer simultaneously to conduct the SWPP test. After measuring the data of C_rec(t^∗) and C_tra(t^∗) in the extraction phase, one could fit the data of $\ln \left(\frac{C_{\mathrm{rec}}\left(t^{\ast }\right)}{C_{\mathrm{tra}}\left(t^{\ast }\right)}\right)\sim t^{\ast }$ using a linear function, and the slope of t^∗ is the estimation of λ. Snodgrass and Kitanidis (1998) derived a similar model for estimating λ:

$$\begin{array}{}\text{(9)}& \ln \left({\displaystyle \frac{C_{\mathrm{rea}}\left(t^{\ast }\right)}{C_{\mathrm{tra}}\left(t^{\ast }\right)}}\right)=\ln \left({\displaystyle \frac{C_{\mathrm{rec}}^{\mathrm{0}}}{C_{\mathrm{tra}}^{\mathrm{0}}}}\right)-\mathit{\lambda }{t}^{\ast }.\end{array}$$

Comparing Eq. (8) with Eq. (9), one could find that the difference is the first terms on the right-hand sides of equations, while λ is the slope for both Eqs. (8) and (9). Although the accuracy of both models has been tested by a number of investigators, previous studies on reactive transport were based on an assumption that the aquifer hydraulic diffusivity was infinite (e.g., Eq. 1 of Reinhard et al., 1997, and Eq. 2 of Haggerty et al., 1998).

Actually, the assumptions of Eqs. (2)–(9) are debatable for the actual applications, and may cause errors in modeling the solute transport in the SWPP test. The second and third assumptions relate to the wellbore storage of the solute transport in the SWPP test. In the following section, the new models will be proposed to investigate the potential errors when these assumptions are involved.

2.2 A revised model with a finite hydraulic diffusivity

As for the first assumption in Sect. 2.1, Wang et al. (2017) demonstrated that it might result in non-negligible errors in parameter estimation, particularly for the estimation of dispersivity. A minor point to note is that the model of Wang et al. (2017) mainly focused on conservative solute transport, rather than reactive transport. Nevertheless, the pore velocity of transient flow is calculated by Darcy's law:

$$\begin{array}{}\text{(10)}& v_{\mathrm{r}}={\displaystyle \frac{K_{\mathrm{r}}}{\mathit{\theta }}}{\displaystyle \frac{\partial s}{\partial r}},\end{array}$$

where K_r is the radial hydraulic conductivity; s is drawdown which could be obtained by solving the following mass balance equation with the proper initial and boundary conditions:

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle \frac{\partial v_{\mathrm{r}}}{\partial r}}+{\displaystyle \frac{v_{\mathrm{r}}}{r}}={\displaystyle \frac{S_{\mathrm{s}}}{\mathit{\theta }}}{\displaystyle \frac{\partial s\left(r,t\right)}{\partial r}},r\ge r_{\mathrm{w}},\text{(12)}& {\displaystyle }& {\displaystyle }{\left.s\left(r,t\right)\right|}_{t=\mathrm{0}}=\mathrm{0},\text{(13)}& {\displaystyle }& {\displaystyle }{\left.v_{\mathrm{r}}\right|}_{r\to \mathrm{\infty }}=\mathrm{0},\text{(14)}& {\displaystyle }& {\displaystyle }{\left.\left(\mathrm{2}\mathit{\pi }B{v}_{\mathrm{r}}\right)\right|}_{r\to r_{\mathrm{w}}}-{\displaystyle \frac{\mathit{\pi }{r}_{\mathrm{w}}^{\mathrm{2}}}{\mathit{\theta }}}{\displaystyle \frac{\mathrm{d}{s}_{\mathrm{w}}\left(t\right)}{\mathrm{d}t}}=Q,\end{array}$$

where S_s is the specific storage of the aquifer; s_w is the drawdown inside the wellbore.

As for the second assumption in the rest phase, as shown in Eq. (3), it implies that the concentration of the solute is zero in the wellbore. This assumption works when the chase concentration is zero and the prepared solution is completely pushed out of the borehole into the aquifer at the end of the chase phase. However, the chase concentration might be non-zero, as demonstrated in Phanikumar and McGuire (2010) and McGuire et al. (2002). Consequently, the concentration in the early stage of the rest phase, which is close to the concentration at the end of the chase phase, is not zero. This is because the water level in the wellbore is greater than the hydraulic head in the surrounding aquifer due to the wellbore storage, resulting in a positive flux from the wellbore into the aquifer. Correspondingly, when the chase concentration is not zero or the prepared solution in the injection phase is not completely pushed out of the wellbore, the concentration in the wellbore may not be zero in the early stage of the rest phase. In this study, we employed the Danckwerts condition for transport at the well screen in the rest period (Danckwerts, 1953):

$$\begin{array}{}\text{(15)}& {\left.{\displaystyle \frac{\partial C}{\partial r}}\right|}_{r\to r_{\mathrm{w}}}=\mathrm{0},t_{\mathrm{inj}}+t_{\mathrm{cha}}<t\le t_{\mathrm{inj}}+t_{\mathrm{cha}}+t_{\mathrm{res}}.\end{array}$$

Actually, Eq. (15) acknowledges the continuity of concentration and continuity of mass flux simultaneously across the well screen, namely ${\left.C\right|}_{r\to r_{\mathrm{w}}^{-}}={\left.C\right|}_{r\to r_{\mathrm{w}}^{+}}$ and ${\left.\left(v_{\mathrm{r}}C\right)\right|}_{r\to r_{\mathrm{w}}^{-}}={\left.\left(v_{\mathrm{r}}C-D_{\mathrm{r}}\frac{\partial C}{\partial r}\right)\right|}_{r\to r_{\mathrm{w}}^{+}}$, where the − and + signs in the subscript of r_w represent approaching of the well screen from inside the well and outside the well, respectively.

The third assumption mentioned in Sect. 2.1 seems not reasonable at the early stage of the injection and chase phases, because the concentration of the injected solute will be affected by the finite volume of water in the wellbore. Take the chase phase as an example: it is impossible to immediately reduce the solute concentration inside the wellbore from a certain level during the tracer injection phase to zero when switching to the chase phase, even when the solute concentration in the chase phase is zero. This is because the wellbore with a finite radius contains a certain finite mass of solute at the moment of switching from injection of a tracer to injection of a chaser. Therefore, it will take some time to completely flush out the residual tracer inside the wellbore after the start of the chase phase, and a larger wellbore will take a longer time to flush out the residual tracer inside the wellbore. This means that the concentration at the wellbore–aquifer interface will not drop to zero immediately after the start of the chase phase. Instead, it will take a finite period of time to gradually approach zero during the chase phase. Similarly, the boundary condition of the well screen in the injection phase might not be appropriate in previous studies if the wellbore storage effect is of concern. Therefore, the value of a solute concentration inside the wellbore should be smaller than or equal to $C_{\mathrm{0}}^{\mathrm{inj}}$ in the injection phase and greater than or equal to $C_{\mathrm{0}}^{\mathrm{cha}}$ in the chase phase.

Here, we will develop a new approach to take care of the concentration in the wellbore in the injection and chase phases based on the mass balance principle, i.e.,

$$\begin{array}{ll}{\displaystyle }\mathrm{\Delta }m& {\displaystyle }=C_{\mathrm{0}}^{\mathrm{inj}}Q\mathrm{\Delta }t=C_{\mathrm{w}}^{t+\mathrm{\Delta }t}\left(V^t+Q\mathrm{\Delta }t\right)-C_{\mathrm{w}}^t{V}^t,\\ \text{(16)}& {\displaystyle }& {\displaystyle }\mathrm{0}<t\le t_{\mathrm{inj}},{\displaystyle }\mathrm{\Delta }m& {\displaystyle }=C_{\mathrm{0}}^{\mathrm{cha}}Q\mathrm{\Delta }t=C_{\mathrm{w}}^{t+\mathrm{\Delta }t}\left(V^t+Q\mathrm{\Delta }t\right)-C_{\mathrm{w}}^t{V}^t,\\ \text{(17)}& {\displaystyle }& {\displaystyle }{t}_{\mathrm{inj}}<t\le t_{\mathrm{inj}}+t_{\mathrm{cha}},\end{array}$$

where Δm represents the mass entering into the well during time interval Δt; $C_{\mathrm{w}}^t$ and $C_{\mathrm{w}}^{t+\mathrm{\Delta }t}$ represent the solute concentrations in the wellbore at times t and t+Δt, respectively; V^t represents the volume of water in the wellbore at time t. The initial values of $C_{\mathrm{w}}^t$ and V^t at the injection phase are

$$\begin{array}{}\text{(18)}& {\displaystyle }& {\displaystyle }{\left.C_{\mathrm{w}}^t\right|}_{t=\mathrm{0}}=\mathrm{0},\text{(19)}& {\displaystyle }& {\displaystyle }{\left.V^t\right|}_{t=\mathrm{0}}=\mathit{\pi }{r}_{\mathrm{w}}^{\mathrm{2}}\left({\left.H_{\mathrm{w}}\right|}_{t=\mathrm{0}}\right),\end{array}$$

where H_w represents the water depth of the wellbore.

In the chase phase, one has

$$\begin{array}{}\text{(20)}& {\displaystyle }& {\displaystyle }{\left.C_{\mathrm{w}}^t\right|}_{t=t_{\mathrm{inj}}^{-}}={\left.C_{\mathrm{w}}^t\right|}_{t=t_{\mathrm{inj}}^{+}},\text{(21)}& {\displaystyle }& {\displaystyle }{\left.V^t\right|}_{t=t_{\mathrm{inj}}^{+}}=\mathit{\pi }{r}_{\mathrm{w}}^{\mathrm{2}}\left({\left.H_{\mathrm{w}}\right|}_{t=t_{\mathrm{inj}}^{+}}\right),\end{array}$$

where the − and + signs in the subscripts of Eqs. (20)–(21) hereinafter represent approaching of the limit from the left- and right-hand sides of t_inj, respectively.

2.3 Capability of the new SWPP model of this study

Different from the model of Wang et al. (2017), the multi-species reactive transport models are used to describe the nonlinear biogeochemical reactive processes considering wellbore effects not only for groundwater flow, but also for solute concentrations. The new model of this study is an extension of Phanikumar and McGuire (2010) that ignored the wellbore storage for both groundwater flow and solute transport, and assumed that the aquifer hydraulic diffusivity was infinite. The Danckwerts condition rather than the Robin condition is applied at the well screen in the rest phase of this study. Therefore, the new model is more powerful in describing an arbitrary number of species and user-defined reaction rate expressions, including Monod/Michaelis–Menten kinetics.

3 Numerical solution of the SWPP test

In this study, we will use a finite-difference method to solve the model of the SWPP test, where the finite-difference scheme of the groundwater flow is the same as Wang et al. (2017), and the scheme of the transport governing equation (ADE) is similar to the model of Phanikumar and McGuire (2010). However, the flow velocity used in the advective term of ADE is computed by solving the model of groundwater flow rather than directly using Eq. (2), which was employed by Phanikumar and McGuire (2010).

To minimize numerical errors and to increase computational efficiency, we employ a non-uniform grid system for simulations (Wang et al., 2014), which is

$$\begin{array}{}\text{(22)}& r_i={\displaystyle \frac{r_{i-\mathrm{1}/\mathrm{2}}+r_{i+\mathrm{1}/\mathrm{2}}}{\mathrm{2}}},i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\cdots },N_{\mathrm{r}},\end{array}$$

where N_r represents the number of nodes in discretization of the spatial domain [r_w, r_e]; r_w and r_e, respectively, represent the distances of inner and outer boundary nodes; r_i is the radial distance of a node; $r_{i+\mathrm{1}/\mathrm{2}}$ is calculated as follows:

$$\begin{array}{ll}{\displaystyle }{\log }_{\mathrm{10}}\left(r_{i+\mathrm{1}/\mathrm{2}}\right)& {\displaystyle }={\log }_{\mathrm{10}}\left(r_{\mathrm{w}}\right)+i\left[{\displaystyle \frac{{\log }_{\mathrm{10}}\left(r_{\mathrm{e}}\right)-{\log }_{\mathrm{10}}\left(r_{\mathrm{w}}\right)}{N}}\right],\\ \text{(23)}& {\displaystyle }& {\displaystyle }i=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\cdots },N_{\mathrm{r}}.\end{array}$$

The value of $r_{i-\mathrm{1}/\mathrm{2}}$ can be calculated using the similar way. Equations (22)–(23) represent a space domain discretized logarithmically, and the spatial steps are smaller near the wellbore and become progressively greater away from the wellbore.

Similarly, we logarithmically discretize the temporal domain:

$$\begin{array}{}\text{(24)}& t_i={\displaystyle \frac{t_{i-\mathrm{1}/\mathrm{2}}+t_{i+\mathrm{1}/\mathrm{2}}}{\mathrm{2}}},i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\cdots },M,\end{array}$$

where M represents the number of nodes in discretization of the temporal domain; t_i is the time of node i; $t_{i+\mathrm{1}/\mathrm{2}}$ is calculated as follows in the injection phase:

$$\begin{array}{ll}{\displaystyle }{\log }_{\mathrm{10}}\left(t_{i+\mathrm{1}/\mathrm{2}}\right)& {\displaystyle }={\log }_{\mathrm{10}}\left(t_{\mathrm{0}}\right)+i\left[{\displaystyle \frac{{\log }_{\mathrm{10}}\left(t_{\mathrm{inj}}\right)-{\log }_{\mathrm{10}}\left(t_{\mathrm{0}}\right)}{M}}\right],\\ \text{(25)}& {\displaystyle }& {\displaystyle }i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\cdots },M,\end{array}$$

where t₀ is a very small positive value representing the first time step, such as $t_{\mathrm{0}}=\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{7}}$ h.

As for the chase, one has

$$\begin{array}{ll}{\displaystyle }{t}_{i+\mathrm{1}/\mathrm{2}}& {\displaystyle }={\mathrm{10}}^{{\log }_{\mathrm{10}}\left(t_{\mathrm{0}}\right)+i\left[\frac{{\log }_{\mathrm{10}}\left(t_{\mathrm{cha}}\right)-{\log }_{\mathrm{10}}\left(t_{\mathrm{0}}\right)}{M}\right]}+t_{\mathrm{inj}},\\ \text{(26)}& {\displaystyle }& {\displaystyle }i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\cdots },M.\end{array}$$

Similarly, in the rest phase, one has

$$\begin{array}{ll}{\displaystyle }{t}_{i+\mathrm{1}/\mathrm{2}}& {\displaystyle }={\mathrm{10}}^{{\log }_{\mathrm{10}}\left(t_{\mathrm{0}}\right)+i\left[\frac{{\log }_{\mathrm{10}}\left(t_{\mathrm{res}}\right)-{\log }_{\mathrm{10}}\left(t_{\mathrm{0}}\right)}{M}\right]}+t_{\mathrm{inj}}+t_{\mathrm{cha}},\\ \text{(27)}& {\displaystyle }& {\displaystyle }i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\cdots },M.\end{array}$$

In the extraction phase, one has

$$\begin{array}{ll}{\displaystyle }{t}_{i+\mathrm{1}/\mathrm{2}}& {\displaystyle }={\mathrm{10}}^{{\log }_{\mathrm{10}}\left(t_{\mathrm{0}}\right)+i\left[\frac{{\log }_{\mathrm{10}}\left(t_{\mathrm{ext}}\right)-{\log }_{\mathrm{10}}\left(t_{\mathrm{0}}\right)}{M}\right]}+t_{\mathrm{inj}}+t_{\mathrm{cha}}+t_{\mathrm{res}},\\ \text{(28)}& {\displaystyle }& {\displaystyle }i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\cdots },M.\end{array}$$

Before using the new model of this study, it is necessary to evaluate the numerical errors (like artificial oscillation and numerical dispersion) of the solution. Unfortunately, the benchmark analytical solutions of the SWPP test with a finite hydraulic diffusivity are not available to date. Alternatively, the accuracy of the finite-difference solution could be tested by comparison with the numerical solution of Wang et al. (2017), which was proven to be accurate and robust. Figure 2 shows the comparison of BTCs between the solution of Wang et al. (2017) and of this study, where the parameters used are similar to Fig. 3 of Phanikumar and McGuire (2010): B=8 m, r_w=0.052 m, α_r=1 m, θ=0.38 m, D₀=0 m² h⁻¹, t_inj=94.32 h, t_cha=0 h, t_res=0 h, t_ext=405.6 h, injection flow rate Q_inj=0.1 m³ h⁻¹, and extraction flow rate $Q_{\mathrm{ext}}=-\mathrm{0.11}$ m³ h⁻¹. It shows a small oscillation in the numerical solutions, which might be caused by the numerical errors.

Figure 2 Comparison of BTCs between the solutions of Wang et al. (2017) and of this study, where $C_{\mathrm{0}}^{\mathrm{inj}}$ represents the concentration of the prepared solute in the injection phase.

By comparing the solution of this study with Wang et al. (2017), one may conclude that the solution of this study appears to be accurate and reliable since the mean square error between two solutions is smaller than 0.05 for all cases in Fig. 2. In the wellbore (r=r_w), the concentration is equal to $C_{\mathrm{0}}^{\mathrm{inj}}$, as shown in Fig. 2. This is due to the boundary condition of the wellbore, e.g.,

$$\begin{array}{}\text{(29)}& {\left.C\right|}_{r\to r_{\mathrm{w}}}=C_{\mathrm{0}}^{\mathrm{inj}},\mathrm{0}<t\le t_{\mathrm{inj}}.\end{array}$$

In the aquifer, the values of BTCs increase with the decreasing distance from the wellbore.

It is also necessary to test the accuracy of the new models against the numerical software packages. Since the code of the original MODFLOW/MT3DMS package is open source and could be downloaded freely from the website of the United States Geological Survey, it is preferred by many modelers and is selected as the basis of comparison in this study. Unfortunately, such an open-source MODFLOW/MT3DMS package may create some errors in describing the solute transport in the wellbore-confined aquifer. The errors come from an assumption that the water volume in the wellbore is computed by a product of the wellbore cross section and the aquifer thickness, which is incorrect. The actual water volume in the wellbore should be computed by a product of the wellbore cross section and the water level in the wellbore (see the Supplement for a detailed explanation). Figure 3 shows the comparisons of BTCs between the open-source MODFLOW/MT3DMS package and the new model of this study. The water level of the wellbore is assumed to be equal to the aquifer thickness in the new models for the purpose of comparison, although it may not be true, and the other parameters used are the same as the ones in Fig. 2. Therefore, the agreement between the two models demonstrates the accuracy of the new model. Figure 3 shows that the concentration in the wellbore is not unit in the injection phase, and this is because the new model considers the wellbore storage for both groundwater flow and solute transport. It is worthwhile pointing out that an advanced version of MODFLOW/MT3DMS, namely MODFLOW-SURFACT, includes a fracture-well package (FWL4 and FWL5) to overcome the problems in the original open-source MODFLOW well package. The FWL4 and FWL5 packages calculate the water volume using simulated heads, not aquifer thicknesses (see the MODFLOW-SURFACT manual, Vol I, Sect. 3.2, Eq. 24 for details). FEFLOW also has a similar package, referred to as a discrete feature to simulate a pumping/extraction well, if one chooses to do so. Additionally, with a FEFLOW model, the model mesh can be highly discretized to accurately represent well dimensions using a subset of elements (in centers). The modeler can assign a porosity of the unit for those elements representing the wells, rather than assuming the same porosity of the surrounding materials. In the future, we will conduct a comprehensive comparative investigation of the method proposed in this study and those of MODFLOW-SURFACT and FEFLOW for understanding the effects of well mixing and wellbore storage for both flow and transport processes involving an aquifer–well system.

Figure 3 Comparison of BTCs between the solutions of MODFLOW/MT3DMS and of this study.

4 Discussions: effect of wellbore storage on the SWPP test under a transient flow field

Revisiting the assumptions used in previous studies as mentioned in Sects. 2.1 and 2.2, one may find that the flow field and the wellbore storage are key factors for the SWPP test. This is not surprising, since the flow velocity is included not only in the advective term, but also in the dispersive term. The wellbore storage which is dependent on the volume of pre-test water in the wellbore may influence the concentration of the solute injected into the wellbore. As the influence of the hydraulic diffusivity solute transport in the SWPP test has been investigated in Wang et al. (2017), in this section, we mainly investigate the influence the wellbore storage on the reactive transport in the SWPP test in the transient flow field.

The variation of the transient flow field is mainly controlled by the hydraulic diffusivity of the aquifer and the wellbore storage. In the following discussion, we choose three representative types of porous media to test the influence of the hydraulic diffusivity on the results of the SWPP test, including fine sand, medium sand, and coarse sand. According to Domenico and Schwartz (1990) and Batu (1998), one could obtain the values of the hydraulic diffusivity for the above-mentioned three types of media: 4.17×10 m² h⁻¹ (with $K_{\mathrm{r}}=\mathrm{4.17}\times {\mathrm{10}}^{-\mathrm{3}}$ m h⁻¹ and $S_{\mathrm{s}}=\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}$ m⁻¹) for the fine sand, 4.17×10² m² h⁻¹ (with $K_{\mathrm{r}}=\mathrm{4.17}\times {\mathrm{10}}^{-\mathrm{2}}$ m h⁻¹ and $S_{\mathrm{s}}=\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}$ m⁻¹) for the medium sand, and 4.17×10⁴ m² h⁻¹ (with $K_{\mathrm{r}}=\mathrm{4.17}\times {\mathrm{10}}^{-\mathrm{1}}$ m h⁻¹ and $S_{\mathrm{s}}=\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{5}}$ m⁻¹) for the coarse sand. Generally, the hydraulic diffusivity of the aquifer correlates with the grain size of the media, and the value is smaller for the smaller grain size, e.g., fine sand.

The parameters related to the solute transport mainly come from the studies of Phanikumar and McGuire (2010), who interpreted the field experimental data of the SWPP test conducted by McGuire et al. (2002). Except for parameters specifically mentioned otherwise, the default values used in the following section are $C_{\mathrm{0}}^{\mathrm{inj}}=\mathrm{100}$ mg L⁻¹, $C_{\mathrm{0}}^{\mathrm{cha}}=\mathrm{10}$ mg L⁻¹, $B=\mathrm{0}.$1 m, r_w=0.0125 m, α_r=0.01 m, θ=0.33, D₀=0 m² h⁻¹, t_inj=0.6 h, t_cha=0.067 h, t_res=0.0333 h, t_ext=3.6 h, Q_inj=0.0333 m³ h⁻¹, Q_cha=0.0255 m³ h⁻¹, and Q_ext=0.0333 m³ h⁻¹, which can be found in Fig. 5 of Phanikumar and McGuire (2010).

4.1 The rest phase

Figure 4a and b show the comparison of BTCs between the Robin and Danckwerts conditions at the wellbore for different porous media, where $C_{\mathrm{0}}^{\mathrm{cha}}=\mathrm{10}$ mg L⁻¹ in Fig. 4a and $C_{\mathrm{0}}^{\mathrm{cha}}=\mathrm{0.0}$ mg L⁻¹ in Fig. 4b. For the purpose of comparison, the boundary conditions at the wellbore in the injection and chase phases are still described by Eqs. (5)–(6).

Figure 4a shows that the difference of BTCs between two boundary conditions is significant at the early stage of the extraction phase when $C_{\mathrm{0}}^{\mathrm{cha}}=\mathrm{10}$ mg L⁻¹, and BTCs of the Danckwerts condition are above BTCs of the Robin condition. With time going, such a difference becomes negligible. As for the curves of the Robin condition, the solute concentration in the wellbore is 0 in the chase phase; correspondingly, the concentration starts from 0 at the early stage of the extraction phase. Actually, the solute concentration in the wellbore may be non-zero in the rest phase due to the wellbore storage and finite hydraulic diffusivity when $C_{\mathrm{0}}^{\mathrm{cha}}=\mathrm{10}$ mg L⁻¹. Another interesting observation is that the properties of the porous media could also influence the difference of BTCs between two boundary conditions. Obviously, a smaller hydraulic diffusivity would result in a larger difference between them; e.g., such a difference is greater for the fine sand aquifer.

Figure 4 Comparison of BTCs in the wellbore between the Robin and Danckwerts conditions: (a) $C_{\mathrm{0}}^{\mathrm{cha}}=\mathrm{10.0}$ mg L⁻¹; (b) $C_{\mathrm{0}}^{\mathrm{cha}}=\mathrm{0}$.

Figure 4b shows the comparison of BTCs for different boundary conditions in the wellbore when $C_{\mathrm{0}}^{\mathrm{cha}}=\mathrm{0.0}$, and one could find that the difference of BTCs between the Robin and Danckwerts conditions is negligible, which implies that the Robin condition performs well when $C_{\mathrm{0}}^{\mathrm{cha}}=\mathrm{0.0}$, while this is not for the case when $C_{\mathrm{0}}^{\mathrm{cha}}\ne \mathrm{0.0}$.

4.2 The injection and chase phases

Figure 5 shows the comparison of BTCs in the wellbore for different boundary conditions and different porous media. The parameters used in this case are the same as the ones in Sect. 4.1. The initial head is 1 m. The boundary condition of the wellbore in the rest phase is described by the Danckwerts condition.

Two interesting observations can be seen from this figure. Firstly, the difference of BTCs between the two boundary conditions at the wellbore is obvious, and such a difference is larger for the medium sand than for the coarse sand, implying that it increases with the decreasing hydraulic diffusivity. Secondly, the values of BTCs obtained from Eqs. (16) to (17) are greater at the early stage of the extraction phase, while the peak values of BTC are smaller. In other words, the model of Eqs. (5)–(6) may underestimate the concentration in the early stage of the extraction phase while overestimating the peak values of BTCs.

These observations can be explained as follows. The model of Eqs. (5)–(6) assumes that the volume of water in the wellbore is negligible, and the concentration in the wellbore is close to 10.0 mg L⁻¹ in the rest phase, due to $C_{\mathrm{0}}^{\mathrm{cha}}=\mathrm{10.0}$ mg L⁻¹. As for the model of Eqs. (16)–(17), the volume of water in the wellbore is non-negligible and could dilute the concentration in the injection phase; i.e., the solute concentration in the wellbore could not immediately rise to $C_{\mathrm{0}}^{\mathrm{inj}}$ at the early stage of the injection phase, thus resulting in smaller peak values of BTCs. Similarly, the concentration in the wellbore could not immediately reduce to $C_{\mathrm{0}}^{\mathrm{cha}}$ at the early stage of the chase phase, which makes the concentration larger at the early stage of the extraction phase based on the model of Eqs. (16)–(17).

5 Uniqueness of estimated parameters

Physical and chemical parameters are important in predicting the contaminant transport in the aquifer, and the values of these parameters are generally estimated by best fitting the observed BTCs in the SWPP test using a simplified model, ignoring a number of relevant factors such as the influences of the flow field and the wellbore storage. The discussions in Sect. 4 demonstrate that the negligence of such factors in reactive transport might cause errors and invalidate the whole parameter estimation exercise. Besides porosity, dispersivity, and reaction rates, the new model of this study appears to be useful for estimating the values of hydraulic conductivity and specific storage by best fitting the observed BTCs in the SWPP test. For instance, the values could be determined by minimizing the sum of absolute differences between the observed and calculated BTCs in the wellbore:

$$\begin{array}{ll}{\displaystyle }F& {\displaystyle }=\sum _{i=\mathrm{1}}^O\left|{\left.C_{\mathrm{CAL}}\left(K_{\mathrm{r}},S_{\mathrm{s}},{\mathit{\alpha }}_{\mathrm{r}},\mathit{\theta },{\mathit{\lambda }}_j,t,r_{\mathrm{w}}\right)\right|}_{t=t_i}\right.\\ \text{(30)}& {\displaystyle }& {\displaystyle }\left.-{\left.C_{\mathrm{OBS}}\left(t,r_{\mathrm{w}}\right)\right|}_{t=t_i}\right|,\end{array}$$

where ${\left.C_{\mathrm{CAL}}\left(K_{\mathrm{r}},S_{\mathrm{s}},{\mathit{\alpha }}_{\mathrm{r}},{\mathit{\lambda }}_j,t,r_{\mathrm{w}}\right)\right|}_{t=t_i}$ and ${\left.C_{\mathrm{OBS}}\left(t,r_{\mathrm{w}}\right)\right|}_{t=t_i}$ represent the concentrations calculated by the new model of the SWPP test and the observed concentrations at t=t_i, respectively; o is the number of observed data.

Figure 5 The BTCs in the wellbore for the different boundary conditions at the wellbore in the injection and chase phases.

Although the number of observation points is usually much greater than the number of parameters needed to be estimated, one may still wonder whether Eq. (30) is practically reliable for estimation of multiple parameters simultaneously. To answer this question, two approaches are employed in the following: sensitivity analysis and uniqueness analysis. The sensitivity analysis is used to check whether the solution is sensitive to the parameters or not, while the uniqueness analysis is to check whether the multiple input parameter values could map to the same output results.

5.1 Sensitivity analysis

McCuen (1985) proposed a sensitivity model of a dependent variable, which was normalized as (Kabala, 2001; Yang and Yeh, 2009)

$$\begin{array}{}\text{(31)}& {\mathrm{SC}}_{i,j}=I_j{\displaystyle \frac{\partial C_i}{\partial I_j}},\end{array}$$

where SC_i,j is the sensitivity coefficient of the jth parameter I_j at the ith time; C_i is the concentration at the ith time. In this study, the differentiation of Eq. (31) will be approximated by a finite-difference scheme:

$$\begin{array}{}\text{(32)}& {\mathrm{SC}}_{i,j}=I_j{\displaystyle \frac{C_i\left(I_j+\mathrm{\Delta }{I}_j\right)-C_i\left(I_j\right)}{\mathrm{\Delta }{I}_j}},\end{array}$$

where ΔI_j is a small increment.

From the mathematical models of the groundwater flow, one may find that both hydraulic conductivity and specific storage could affect the flow field. Since greater hydraulic conductivity or smaller specific storage could shorten the time in approaching the steady state, we will employ the hydraulic diffusivity for the sensitivity analysis, which is the ratio of the two parameters. Figure 6 shows the sensitivity of the hydraulic diffusivity on BTCs, and one may find that it is not sensitive to the hydraulic diffusivity when the values of hydraulic diffusivity are sufficiently large. This might be because the time in approaching the steady state is very short when the hydraulic diffusivity values are sufficiently large (for instance, greater than 4.17×10² m² h⁻¹), and the influence of the transient flow could be ignored. Therefore, the steady-state assumption could be used to approximate the flow field in the SWPP test when the hydraulic diffusivity is greater than 4.17×10² m² h⁻¹. Otherwise, the steady-state assumption is not recommended. Figure 7 shows that the BTCs in the wellbore are sensitive to both dispersivity and porosity.

Figure 6 Sensitivity analysis of the hydraulic diffusivity on BTCs in the extraction phase.

5.2 Uniqueness analysis of physical parameters

Besides the sensitivity analysis, the uniqueness analysis is also important for the parameter estimation, which is used to check whether there exist two or more sets of parameters for the same BTCs. Similar to the treatment in previous studies, we firstly use the transient model of this study to reproduce BTCs based on a set of given input parameters, and then estimate the values of parameters by best fitting such BTCs. If the values of the input parameters are different from the estimated parameter when the fitness is very good, one could conclude that the solution is not unique and the parameters estimated from Eq. (30) may not be reliable.

There are four physical parameters in the new model of this study, i.e., hydraulic conductivity, specific storage, dispersivity, and porosity, and one chemical parameter (reaction rate). Wang et al. (2017) investigated the uniqueness of solutions for the flow field, and the results showed that BTCs of the SWPP test were not unique for the flow-related parameters. For instance, BTCs with a steady-state flow field were almost the same as BTCs with a transient flow field, as shown in Figs. 10 and 11 of Wang et al. (2017). It implies that one may not inversely determine the hydraulic parameters of a flow field only by best fitting observed BTCs in the wellbore, and additional aquifer tests are required to supplement the SWPP test to determine the flow-related parameters. However, Wang et al. (2017) did not investigate the uniqueness of porosity and dispersion when the hydraulic parameters were given, which will be discussed in this study.

Figure 7 Sensitivity analysis of dispersivity and porosity on BTCs in the extraction phase.

Figure 8 shows comparison of BTCs for different dispersivities and porosities but for the same hydraulic parameters, and one could see that the curves of α_r=0.01 m and θ=0.33 are almost the same as the curves of α_r=0.006 m and θ=0.9. Therefore, Eq. (30) may not be used to determine α_r and θ simultaneously. Fortunately, the porosity could be measured in the laboratory from core samples or determined by the SWPP test with drift flow (Hall et al., 1991; Paradis et al., 2018). When the values of K_r, S_s, and θ are given, the dispersivity could be determined uniquely by Eq. (30).

In summary, it seems impossible to determine all parameters (K_r, S_sα_r, and θ) simultaneously by only best fitting the observed BTCs in the wellbore of the SWPP test using Eq. (30). Therefore, before determining the parameters related to the solute transport (α_r and θ), the hydraulic parameters (K_r and S_s) needed to be estimated by supplementary aquifer tests, or by best fitting the pressure data measured during the SWPP test, e.g., the pumping phase. The value of α_r could be determined by Eq. (30) when the porosity is given.

Table 1 Reaction parameters estimated by linear functions.

5.3 Chemical parameter estimation

The models estimating the reaction rate are based on several assumptions in previous studies, e.g., Eqs. (8)–(9) as demonstrated in Sect. 2.1. To test the applicability of those equations, we will use the model of this study to reproduce the data of $\ln \left(C_{\mathrm{rec}}/C_{\mathrm{tra}}\right)\sim t^{\ast }$ based on a set of given parameters, and then using Eqs. (8)–(9) (which is based on an infinite hydraulic diffusivity presumption) to estimate λ (denoted as $\stackrel{\mathrm{̃}}{\mathit{\lambda }}$) by best fitting $\ln \left(C_{\mathrm{rec}}/C_{\mathrm{tra}}\right)\sim t^{\ast }$. Two species involved in this case are Cl⁻¹ and ${\mathrm{SO}}_{\mathrm{4}}^{-\mathrm{2}}$, in which C_tra and C_rec represent the concentrations of Cl⁻¹ and ${\mathrm{SO}}_{\mathrm{4}}^{-\mathrm{2}}$, respectively. Figure 9 shows the fitness of the simulated $\ln \left(C_{\mathrm{rec}}/C_{\mathrm{tra}}\right)\sim t^{\ast }$ in the wellbore using a linear function, with the detailed information shown in Table 1. Two sets of λ are employed in the discussions for the reactant, e.g., λ=0.1 and 0.2 h⁻¹. One may conclude that the simplified models of Eqs. (8)–(9) with an infinite hydraulic diffusivity perform well in the estimation of λ for reactive transport under the finite hydraulic diffusivity condition.

Figure 8 Comparison of BTCs for different dispersivities and porosities but for the same hydraulic parameters.

This simplified model of Eq. (9) has been widely used to estimate λ, due to the advantages that λ could be determined directly by best fitting the observed $\ln \left(C_{\mathrm{rec}}/C_{\mathrm{tra}}\right)\sim t^{\ast }$, without knowledge of the aquifer properties, such as porosity, dispersivity, and hydraulic diffusivity. However, this model is proposed based on the first-order reaction assumption, which is a linear function as shown in Eq. (7). Whether this model works for nonlinear reactions or not is still unknown and will be investigated in the following section.

Assuming that the extraction time since the rest phase ended could be divided into N−1 segments, Phanikumar and McGuire (2010) employed the Heaviside unit step function to describe a type of nonlinear biogeochemical reaction:

$$\begin{array}{}\text{(33)}& {\displaystyle \frac{\partial C}{\partial t}}=-\sum _{j=\mathrm{1}}^{N-\mathrm{1}}\left[H\left(t-t_j^{\ast }\right)-H\left(t-t_{j+\mathrm{1}}^{\ast }\right)\right]{\mathit{\lambda }}_j{C}_i^{n_j},\end{array}$$

where λ_j is the reaction constant in the temporal segment j, and the Heaviside step function H(⋅) is

$$\begin{array}{}\text{(34)}& H\left(t-t_j^{\ast }\right)-H\left(t-t_{j+\mathrm{1}}^{\ast }\right)=\left\{\begin{array}{l}\mathrm{0}\text{if}t<t_j^{\ast }\\ \mathrm{1}\text{if}{t}_j^{\ast }<t<t_{j+\mathrm{1}}^{\ast }\\ \mathrm{0}\text{if}{t}_{j+\mathrm{1}}^{\ast }<t\end{array}\right..\end{array}$$

Equation (33) is a series of piece-wise linear (n_j=1) or nonlinear (n_j≠1) functions, which are an extension of Eq. (7).

Figure 9 Fitness of $\ln \left[C_{\mathrm{rec}}/C_{\mathrm{tra}}\right]\sim t^{\ast }$ produced by the numerical solution of this study with the first-order reaction in the different porous media.

Figure 10 Computed $\ln \left[C_{\mathrm{rec}}/C_{\mathrm{tra}}\right]\sim t^{\ast }$ by the model of this study using a piece-wise linear function to describe the nonlinear chemical reactions.

To test the influence of the hydraulic diffusivity on the accuracy of this model in estimating λ_j for the nonlinear reactions, the model of this study is used to reproduce the data of $\ln \left(C_{\mathrm{rec}}/C_{\mathrm{tra}}\right)\sim t^{\ast }$ with a set of specific λ_j, n_j and $t_j^{\ast }$ for three types of porous media. Figures 10 and 11 represent the computed $\ln \left(C_{\mathrm{rec}}/C_{\mathrm{tra}}\right)\sim t^{\ast }$ based on the model of the chemical reactions described by the piece-wise linear and nonlinear functions, respectively. The values of λ_j and $t_j^{\ast }$ of Fig. 10 are obtained by best fitting the observation data using a piece-wise linear function (e.g., n_j=1) proposed by Phanikumar and McGuire (2010). The circle represents the experiment data observed by McGuire et al. (2002). The parameters related to the chemical reactions in Fig. 11 are from Phanikumar and McGuire (2010) by best fitting the observation data using a nonlinear function: λ_j=0.25, n_j=0.25, $N=j=\mathrm{1}$. Comparing Figs. 10 and 11, one may find that the influence of the hydraulic diffusivity on the computed $\ln \left(C_{\mathrm{rec}}/C_{\mathrm{tra}}\right)\sim t^{\ast }$ is negligible for the chemical reaction described by the piece-wise linear function, which is similar to the first-order reaction as shown in Fig. 9. However, the influence of the hydraulic diffusivity on the relationship of $\ln \left(C_{\mathrm{rec}}/C_{\mathrm{tra}}\right)\sim t^{\ast }$ cannot be ignored if one attempts to use nonlinear functions to describe such a chemical reaction. The difference between the curves of different porous media is obvious in Fig. 11. The agreement between the observed and computed data is satisfactory for the medium and coarse sands, but not for the fine sand in Fig. 11. This is because the hydraulic diffusivity values of the medium and coarse sands are larger than that of the fine sand, and thus are close to the assumption of an infinite hydraulic diffusivity used in Phanikumar and McGuire (2010).

Therefore, one may conclude that λ_j, n_j and $t_j^{\ast }$ could be determined by directly best fitting the observed $\ln \left(C_{\mathrm{rec}}/C_{\mathrm{tra}}\right)\sim t^{\ast }$ when the nonlinear reactions are described by the piece-wise linear functions, in a similar way to estimating the linear reaction rate by Eq. (7). However, such an approach may not work if one attempts to use nonlinear functions to describe such reactions.

Figure 11 Computed $\ln \left[C_{\mathrm{rec}}/C_{\mathrm{tra}}\right]\sim t^{\ast }$ by the model of this study using a nonlinear function to describe the nonlinear chemical reactions.

Figure 12 BTCs for the different porous media with a piece-wise linear function to describe chemical reactions: (a) Cl¹⁻ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in the aquifer at $r=r_{\mathrm{w}}+\mathrm{0.15}$ m, (b) Cl¹⁻ in the wellbore.

6 Field applications

To test the model of this study, the field data of a SWPP test conducted in a single well by McGuire et al. (2002) will be employed. In this test, the prepared solution contains Na₂SO₄ (as a reactant) and NaCl (as a conservative tracer). The reactant and the tracer were well mixed and then injected into a targeted aquifer.

6.1 Revisit of the previous model

Phanikumar and McGuire (2010) interpreted such data using a model containing several assumptions mentioned in Sect. 2.1. The parameters used in their model were B=0.1 m, r_w=0.0125 m, α_r=0.001 m, θ=0.33 m, D₀=0 m² h⁻¹, t_inj=0.6 h, t_cha=0.067 h, t_res=0.0333 h, t_ext=3.6 h, Q_inj=0.0333 m³ h⁻¹, Q_cha=0.0255 m³ h⁻¹, and $Q_{\mathrm{ext}}=-\mathrm{0.011}$ m³ h⁻¹. The concentrations of NaCl were $C_{\mathrm{0}}^{\mathrm{inj}}=\mathrm{100}$ mg L⁻¹ in the injection phase and $C_{\mathrm{0}}^{\mathrm{cha}}=\mathrm{10}$ mg L⁻¹ in the chase phase. As for the reactant of Na₂SO₄, the concentrations were $C_{\mathrm{0}}^{\mathrm{inj}}=\mathrm{20}$ and $C_{\mathrm{0}}^{\mathrm{cha}}=\mathrm{2}$ mg L⁻¹.

To demonstrate the importance of the wellbore storage of the solute transport, which was ignored in Wang et al. (2017), the observed and computed BTCs are compared based on the estimated parameters in Phanikumar and McGuire (2010), as shown in Fig. 12a and b. The computed BTCs in Fig. 12a and b are located at $r=r_{\mathrm{w}}+\mathrm{0.15}$ m and r=r_w, respectively. The legend of “PPTEST” represents the solution of Phanikumar and McGuire (2010), and the others are produced by the new model, ignoring the wellbore storage effect on the solute transport.

The results showed that the fitness between the observed BTCs in the wellbore and computed BTCs by “PPTEST” was very good, as shown in Fig. 12a of this study or Fig. 5 of Phanikumar and McGuire (2010). However, by carefully checking the report of Phanikumar (2010), we found that the computed BTCs were at a radial distance of 0.15 m from the wellbore, rather than at the wellbore itself in Phanikumar (2010). They did not provide a convincing argument why to choose BTCs in the aquifer to represent BTCs in the wellbore, and thus the use of “0.15 m” in their analysis appears to be an artifact, rather than being physically based. Figure 12b shows the comparison of the computed and observed BTCs in the wellbore for different hydraulic diffusivities. Obviously, the new model ignoring the wellbore storage of the solute transport could not be used to interpret experimental data, since the computed BTCs are zero at the early stage of the extraction phase.

Figure 13 Spatial distribution of the flow velocity in the extraction phase.

Figure 14 Fitness of the field SWPP test data by the new model of this study.

From Fig. 12a and b, several interesting observations could be made. Firstly, the difference of BTCs among different porous media is obvious. BTCs of the coarse sand aquifer are close to the solution of “PPTEST”, as shown in Fig. 12a. This is because the hydraulic diffusivity of the coarse sand aquifer is the largest, which is close to the assumption used in “PPTEST” that hydraulic diffusivity is infinity. Secondly, the wellbore concentration is 10 mg L⁻¹ at the early stage of the extraction phase for Cl⁻. This is mainly due to the chosen boundary condition at the well screen, which has been discussed in detail in Sect. 4.1. Thirdly, the peak values of BTCs increase with the decreasing hydraulic diffusivity, and the arrival times of peak values increase with the decreasing hydraulic diffusivity. Such an observation is also found in Fig. 4a and b. Fourthly, the configuration of BTCs in the aquifer (at $r=r_{\mathrm{w}}+\mathrm{0.15}$ m) computed by the model of this study shows that the concentration firstly decreases with time and then increases with time, as shown in Fig. 12a. This observation could be explained by the corresponding flow field, as shown in Fig. 13. Looking at the flow velocity in the aquifer at $r=r_{\mathrm{w}}+\mathrm{0.15}$ m, one may find that the flow direction is still outward from the wellbore in the early stage of the extraction phase, due to the finite hydraulic diffusivity. The outward flow will persist for a finite period of time, depending on the value of the hydraulic diffusivity, and then reverse its direction to flow towards the wellbore for the rest of the extraction phase. This feature is very different from the results with an infinite hydraulic diffusivity assumption, in which the flow direction is always towards the wellbore for the entire extraction phase.

6.2 Fitness of this study

We try to use the new model to interpret BTCs of the SWPP test, considering a finite hydraulic diffusivity, a finite wellbore storage, and new boundary conditions of the wellbore at the injection, chase and rest phases, assuming the initial head of the flow field is 1 m. In a trial-and-error process of best fitting the observed BTC data, we only estimate parameters of K_r, S_s and α_r, while the other parameters are the same as those used to produce Fig. 5 of Phanikumar and McGuire (2010). Figure 14 demonstrates the fitness of the observed BTC data in the wellbore when K_r=1.0 m h⁻¹, $S_{\mathrm{s}}=\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{5}}$ m⁻¹ and α_r=0.015 m. Since the hydraulic diffusivity of this case is greater than the hydraulic diffusivity of medium sand (4.17×10² m² h⁻¹), the influence of the flow field could be negligible. In this study, we mainly estimated the value of dispersivity, where the porosity is fixed and comes from the reference of Phanikumar and McGuire (2010). Therefore, the dispersivity is uniquely determined.

7 Summary and conclusions

A complete SWPP test includes injection, chase, rest and extraction phases, where the second and third phases are not necessary but are recommended to increase the duration of reaction. Due to the complex mechanics of biogeochemical reactions, aquifer properties, and so on, previous mathematical or numerical models contain some assumptions which may oversimplify the actual physics; for instance, the hydraulic diffusivity of the aquifer is infinite. The Robin or the third-type boundary condition was often used in previous studies at the well screen in the injection, chase and rest phases by ignoring the mixing effect of the volume of water in the wellbore (namely, wellbore storage). In this study, we presented a multi-species reactive SWPP model considering the wellbore storage for both groundwater flow and solute transport, and a finite aquifer hydraulic diffusivity. The models of wellbore storage for both solute transports are derived based on the mass balance. The Danckwerts boundary condition instead of the Robin condition is employed for solute transport across the well screen in the rest phase. The robustness of the new model is tested by the field data. Meanwhile, the sensitivity analysis and uniqueness analysis of BTCs in wellbore are conducted. The following conclusions can be drawn from this study.

1. The influence of wellbore storage for the solute transport increases with the decreasing hydraulic diffusivity in the injection and chase phases, and the model of Eqs. (16)–(17) underestimates the concentration in the early stage of the injection phase while overestimating the peak values of BTCs.

2. The values of λ_j, n_j and $t_j^{\ast }$ could be determined by directly best fitting the observed $\ln \left(C_{\mathrm{rec}}/C_{\mathrm{tra}}\right)\sim t^{\ast }$ when the nonlinear reactions are described by the piece-wise linear functions, while such an approach may not work if one attempts to use nonlinear functions to describe such nonlinear reactions.

3. The Robin condition used to describe the wellbore flux in the rest phase works well only when the chase concentration is zero and the prepared solution in the injection phase is completely pushed out of the borehole into the aquifer, while the Danckwerts boundary condition performs better even when the chase concentration is not zero.

4. In the extraction phase, the peak values of BTCs increase with the decreasing hydraulic diffusivity, and the arrival time of the peak value becomes shorter when the hydraulic diffusivity is smaller.

5. It seems impossible to determine all parameters simultaneously by only best fitting the observed BTCs in the wellbore of the SWPP test using Eq. (30). The hydraulic parameters needed to be estimated by supplementary aquifer tests before determining the parameters related to the solute transport. The value of α_r could be determined by Eq. (30) when the porosity is given.

Appendix A: Nomenclature

B	Aquifer thickness (L)
Ci	Aqueous phase concentration of the ith reactive solute (ML−3)
C	Resident concentration of the aqueous phase to represent Ci in Eq. (1) (ML−3)
$C_{\mathrm{0}}^{\mathrm{inj}},C_{\mathrm{0}}^{\mathrm{cha}}$	Solute concentrations injected into the wellbore during the injection and chase phases (ML−3), respectively
$C_{\mathrm{w}}^t,C_{\mathrm{w}}^{t+\mathrm{\Delta }t}$	Solute concentrations in the wellbore at time t and t+Δt (ML−3), respectively
$C_{\mathrm{rec}}\left(t^{\ast }\right),C_{\mathrm{tra}}\left(t^{\ast }\right)$	Reactant concentration and the concentration of a conservative tracer (ML−3), respectively
Dr	Dispersion coefficient (L2 T−1)
et∗	Time since the end of injection (T)
Fj	Monod/Michaelis–Menten kinetics function (dimensionless)
Kr	Radial hydraulic conductivity (LT−1)
M	Number of nodes in discretization of the temporal domain (dimensionless)
Δ m	Mass entering into the well during time interval Δt (M)
N	Number of the segment for chemical reactions (dimensionless)
Nr	Number of nodes in discretization of the spatial domain (dimensionless)
Q	Flow rate of the well (L3 T−1)
r	Radius distance from the center of the well (L)
ri	Radial distance of node (L)
rw	Well radius (L)
re	Distance of the outer boundary of the aquifer (L)
s	Drawdown (L)
sw	Drawdown inside the wellbore (L)
Si	Solid phase concentration of the ith reactive solute (ML−3)
Ss	Specific storage of aquifer (L−1)
t	Time in the SWPP test (T)
ti	Time of node i (T)
tinj, tcha, tres, text	Durations (T) of the injection, chase, rest, and extraction phases, respectively
t∗	Time since the end of injection (T)
$t_j^{\ast }$	Times at two ends of segment j (T)
ur	Radial Darcian velocities (LT−1)
$v_{\mathrm{r}}=u_{\mathrm{r}}/\mathit{\theta }$	Average radial pore velocity (LT−1)
Vt	Volume of water in the wellbore at the time t (L3)
ρb	Bulk density of the aquifer material (ML−3)
θ	Porosity (dimensionless)
H	Heaviside step function (dimensionless)
nj, λj	Orders (dimensionless) and constant (dimensionless) in the temporal segment j, respectively
αr	Radial dispersivity (L)
λ	First-order reaction rate constant (dimensionless)
ADE	Advection dispersion equation
BTC	Breakthrough curve
PPTEST	Solution of Phanikumar and McGuire (2010)
SWPP	Single well push–pull