Interactive comment on “ Reactive Transport with Wellbore Storages in a Single-Well Push-Pull Test ” by Quanrong Wang and Hongbin

Response to Referee #1: General Comments: 1. Based on Wang et al. (2017) and Phanikumar and McGuire (2010), this study proposes an extended model to interpret single-well push-pull (SWPP) tests in homogeneous reservoirs considering transient flow, well-bore storage and reactive transport. The improvements from Wang et al. (2017) concern the change in boundary conditions to better consider the effect of well-bore storage and the extension to reactive transport involving complex chemical reactions based on Phanikumar and McGuire (2010). To show the interest of this extended model, the authors revisited the interpretation of a published SWPP test data set by McGuire et al. (2002) and compared with the model of Phanikumar and McGuire (2010) (available on-line) used to interpret this data set. The authors highlighted that


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 −1 for Cl − and 2 mg L −1 for SO 2− 4 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.

Revisit of the previous model
The general form of the governing equation for a multispecies reactive SWPP test is 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 and t * j +1 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): 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 r = u r /θ 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): 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 in which α r is the radial dispersivity; D 0 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): or or where C inj 0 and C cha 0 represent the solute concentrations injected into the wellbore during the injection and chase phases, respectively.A detailed discussion about the abovementioned 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 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) 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, andEq. 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.

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: 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 condi-tions: 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 nonzero, 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): Actually, Eq. ( 15) acknowledges the continuity of concentration and continuity of mass flux simultaneously across the well screen, namely , 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 wellboreaquifer 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 inj 0 in the injection phase and greater than or equal to C cha 0 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., where m represents the mass entering into the well during time interval t; C t w and C t+ t w 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 t w and V t at the injection phase are where H w represents the water depth of the wellbore.
In the chase phase, one has where the − and + signs in the subscripts of Eqs. ( 20)-( 21) hereinafter represent approaching of the limit from the leftand right-hand sides of t inj , respectively.

Capability of the new SWPP model of this study
Different from the model of Wang et al. (2017), the multispecies 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.

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 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+1/2 is calculated as follows: The value of r i−1/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: where M represents the number of nodes in discretization of the temporal domain; t i is the time of node i; t i+1/2 is calculated as follows in the injection phase: where t 0 is a very small positive value representing the first time step, such as t 0 = 1.0 × 10 −7 h.As for the chase, one has Similarly, in the rest phase, one has In the extraction phase, one has 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  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 inj 0 , as shown in Fig. 2.This is due to the boundary condition of the wellbore, e.g., 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 MOD-FLOW/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).FE-FLOW 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.
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 2 h −1 (with K r = 4.17 × 10 −3 m h −1 and S s = 1.0 × 10 −4 m −1 ) for the fine sand, 4.17 × 10 2 m 2 h −1 (with K r = 4.17 × 10 −2 m h −1 and S s = 1.0 × 10 −4 m −1 ) for the medium sand, and 4.17 × 10 4 m 2 h −1 (with K r = 4.17 × 10 −1 m h −1 and S s = 1.0 × 10 −5 m −1 ) 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 inj 0 = 100 mg L −1 , C cha 0 = 10 mg L −1 , B = 0.1 m, r w = 0.0125 m, α r = 0.01 m, θ = 0.33, D 0 = 0 m 2 h −1 , 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 3 h −1 , Q cha = 0.0255 m 3 h −1 , and Q ext = 0.0333 m 3 h −1 , which can be found in Fig. 5 of Phanikumar and McGuire (2010).

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 cha 0 = 10 mg L −1 in Fig. 4a and C cha 0 = 0.0 mg L −1 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 cha 0 = 10 mg L −1 , 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 cha 0 = 10 mg L −1 .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 4b shows the comparison of BTCs for different boundary conditions in the wellbore when C cha 0 = 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 cha 0 = 0.0, while this is not for the case when C cha 0 = 0.0.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.

The injection and chase phases
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 −1 in the rest phase, due to C cha 0 = 10.0 mg L −1 .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 inj 0 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 cha 0 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).

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: where C CAL K r , S s , α r , λ j , t, r w t=t i and C OBS (t, r w )| 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.
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.(1985) proposed a sensitivity model of a dependent variable, which was normalized as (Kabala, 2001;Yang and Yeh, 2009)

McCuen
where SC i,j is the sensitivity coefficient of the j th 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: 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 2 m 2 h −1 ), 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 2 m 2 h −1 .Otherwise, the steady-state assumption is not recommended.Figure 7 shows that the BTCs in the wellbore are sensitive to both dispersivity and porosity.

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 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 param-Q.Wang and H. Zhan: Reactive transport with wellbore storages eters (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.

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 (C rec /C tra ) ∼ t * 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 λ) by best fitting ln (C rec /C tra ) ∼ t * .Two species involved in this case are Cl −1 and SO −2 4 , in which C tra and C rec represent the concentrations of Cl −1 and SO −2 4 , respectively.Figure 9 shows the fitness of the simulated ln (C rec /C tra ) ∼ t * 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 −1 .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.
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 (C rec /C tra ) ∼ t * , 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: where λ j is the reaction constant in the temporal segment j , and the Heaviside step function H (•) is 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).
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 (C rec /C tra ) ∼ t * with a set of specific λ j , n j and t * j for three types of porous media.Figures 10 and 11 represent the computed ln (C rec /C tra ) ∼ t * 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 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 = 1.Comparing Figs. 10 and 11, one may find that the influence of the hydraulic diffusivity on the computed ln (C rec /C tra ) ∼ t * 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 (C rec /C tra ) ∼ t * 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 could be determined by directly best fitting the observed ln (C rec /C tra ) ∼ t * 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.

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 2 SO 4 (as a reactant) and NaCl (as a conservative tracer).The reactant and the tracer were well mixed and then injected into a targeted aquifer.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.

Revisit of the previous model
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 −1 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 w + 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 w + 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.

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 −1 , S s = 1.0 × 10 −5 m −1 and α r =0.015 m.Since the hydraulic diffusivity of this case is greater than the hydraulic diffusivity of medium sand (4.17 × 10 2 m 2 h −1 ), 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.

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 could be determined by directly best fitting the observed ln (C rec /C tra ) ∼ t * 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.

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

Figure 2 .
Figure 2. Comparison of BTCs between the solutions of Wang et al. (2017) and of this study, where C inj 0 represents the concentration of the prepared solute in the injection phase.

Figure 3 .
Figure 3.Comparison of BTCs between the solutions of MOD-FLOW/MT3DMS and of this study.

Figure 5
Figure5shows 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.

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

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

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

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

KFigure 9 .
Figure 9. Fitness of ln C rec /C tra ∼ t * produced by the numerical solution of this study with the first-order reaction in the different porous media.

Figure 10 .
Figure 10.Computed ln C rec /C tra ∼ t * by the model of this study using a piece-wise linear function to describe the nonlinear chemical reactions.

Figure 11 .
Figure 11.Computed ln C rec /C tra ∼ t * by the model of this study using a nonlinear function to describe the nonlinear chemical reactions.

Figure 12 .
Figure 12.BTCs for the different porous media with a piece-wise linear function to describe chemical reactions: (a) Cl 1− and SO 2−4 Phanikumar and McGuire (2010) interpreted such data using a model containing several assumptions mentioned in Sect.2.1.The parameters used in their model wereB = 0.1 m, r w = 0.0125 m, α r = 0.001 m, θ = 0.33 m, D 0 = 0 m 2 h −1 , 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 3 h −1 , Q cha = 0.0255 m 3 h −1 , and Q ext = −0.011m 3 h −1 .The concentrations of NaCl were C inj 0 = 100 mg L −1 in the injection phase and C cha 0 = 10 mg L −1 in the chase phase.As for the reactant of Na 2 SO 4 , the concentrations were C inj 0 = 20 and C cha 0 = 2 mg L −1 .To demonstrate the importance of the wellbore storage of the solute transport, which was ignored inWang et al. (2017),

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

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

Table 1 .
Reaction parameters estimated by linear functions.