A general model of radial dispersion with wellbore mixing and skin effects

The mechanism of radial dispersion is essential for understanding reactive transport in the subsurface and for estimating the aquifer parameters required in the optimization design of remediation strategies. Many previous studies demonstrated that the injected solute firstly experienced a mixing process in the injection wellbore, then entered a skin zone after leaving the injection wellbore, and finally moved into the aquifer through advective, diffusive, dispersive, and chemical–biological–radiological processes. In this study, a physically based new model and the associated analytical solutions in the Laplace domain are developed by considering the mixing effect, skin effect, scale effect, aquitard effect, and media heterogeneity (in which the solute transport is described in a mobile–immobile framework). This new model is tested against a finite-element numerical model and experimental data. The results demonstrate that the new model performs better than previous models of radial dispersion in interpreting the experimental data. To prioritize the influences of different parameters on the breakthrough curves, a sensitivity analysis is conducted. The results show that the model is sensitive to the mobile porosity and wellbore volume, and the sensitivity coefficient of the wellbore volume increases with the well radius, while it decreases with increasing distance from the wellbore. The new model represents the most recent advancement in radial dispersion study, incorporating many essential processes not considered in previous investigations.

1 Introduction

Radial dispersion refers to a process of reactive transport under the radial flow condition. One unique feature of radial dispersion (as compared to unilateral dispersion, where the flow velocity is unilateral) is that the dispersive transport becomes progressively weaker when the radial distance from the injection/pumping well becomes larger (or the radial flow velocity becomes smaller), and thus the relative importance of molecular diffusion (which is assumed to be constant) versus the dispersion becomes progressively more robust with a more significant radial distance. The radial dispersion problem is both theoretically interesting and practically important in many fields, like chemical engineering (Davis and Davis, 2002), environmental science (Reinhard et al., 1997; Chen et al., 2016), and hydrogeology (Webster et al., 1970). Although numerical modeling is probably inevitable and more powerful than the analytical modeling in describing radial dispersion, especially put forward for heterogeneous aquifers with complex initial and boundary conditions, the numerical errors and computational cost are not always trivial issues and have to be considered by the engineers. As an alternative, many analytical models have been developed for radial dispersion around an injection well under rather simplified conditions. Such analytical models can fulfill a host of tasks, such as (1) prioritizing the importance of different controlling parameters through a sensitivity analysis, (2) benchmarking the numerical solutions to elucidate the possible numerical errors such as numerical dispersion and artificial oscillation, which are notorious for advection-dominated transport problems, and (3) providing a quick screening tool before implementing a full-scale comprehensive study.

Because of the benefits mentioned above, significant efforts have been put forward over many decades in developing advanced analytical models of radial dispersion. Some examples include the works of Hoopes and Harleman (1967), Gelhar and Collins (1971), Tang and Babu (1979), Moench and Ogata (1981), Chen (1985, 1986, 1991), Hsieh (1986), Tang and Peaceman (1987), Yates (1988), Falade and Brigham (1989), Novakowski (1992), Philip (1994), Veling (2001, 2011), Huang and Goltz (2006), Chen et al. (2007, 2011, 2012, 2017), Gao et al. (2009a), Cihan and Tyner (2011), Wang and Zhan (2013), Hsieh and Yeh (2014), Zhou et al. (2017), Wang et al. (2018, 2020), Huang et al. (2019), and Li et al. (2020). A general trend of such developments is to provide more robust models that can better represent physical reality. However, despite the enormous efforts to date, some significant pitfalls still exist and become roadblocks to quick and accurate interpretation of observed data in the experiments. A primary task of this research is to eliminate such pitfalls, which are briefly illustrated in the following.

A well–aquifer system with radial dispersion comprises a wellbore, skin zone, and aquifer formation zone. The skin zone refers to the disturbed region around the well caused by drilling and construction practices or well completion (Yeh and Chang, 2013; Chen et al., 2012; Li et al., 2019, 2020; Huang et al., 2019). It is spatially between the well screen and the aquifer formation zone. Correspondingly, the injected solute may experience three processes from the wellbore to the aquifer formation zone.

Firstly, the injected solute goes through a mixing process with native (or pre-injection) water in the wellbore at the early injection stage, which is called the mixing effect. Probably due to the small radius of the well, the mixing effect has been overlooked by almost all the analytical solutions mentioned above except Novakowski (1992), Wang et al. (2018), Shi et al. (2020), and Wang et al. (2020), e.g., either by assuming that the well radius was infinitesimal or assuming that the solute concentration in the wellbore was the same as the concentration of the injected solution (Hoopes and Harleman, 1967; Veling, 2011; Zhou et al., 2017). Consequently, the solutions developed without considering the wellbore mixing effect may overestimate concentration values in both the wellbore and the aquifer (Novakowski, 1992; Wang et al., 2018; Shi et al., 2020; Wang et al., 2020). The reason is that the solute concentration in the wellbore is initially zero (when the aquifer is free of solute before the injection) and then increases steadily until it is up to the maximum, which is equal to the concentration of the injected solution.

Secondly, the solute enters the skin zone after leaving the wellbore. Compared with the aquifer formation zone of interest, the dimension of the skin zone is much smaller, e.g., ranging from 0.1 m to several meters, and it is ignored or included in the wellbore. In other words, the effect of the skin zone on radial dispersion (named the skin effect) was negligible. However, numerous previous studies demonstrated that the existence of a skin zone might significantly alter the mechanism of groundwater flow and solute transport around a well (Chen et al., 2012; Hsieh and Yeh, 2014; Yeh and Chang, 2013; Li et al., 2019, 2020). This is because the physical properties (such as permeability, porosity, or dispersivity) of the skin zone are often vastly different from their counterparts in the formation zone. Previously, studies on the skin effect mainly concentrated on the groundwater flow process around the well, and they paid less attention to solute transport processes. To date, few studies have considered the skin effect among the abovementioned analytical models on radial dispersion, such as Chen et al. (2012), Hsieh and Yeh (2014), Huang et al. (2019), and Li et al. (2020). Chen et al. (2012) proposed an analytical solution of solute transport with the skin effect to investigate the influences of dispersivity on radial dispersion; soon after, Hsieh and Yeh (2014) extended the model of Chen et al. (2012) by taking into account a third type (Robin) of condition. Huang et al. (2019) demonstrated that the skin effect significantly influences observed breakthrough curves (BTCs) for radially convergent tracer tests. Recently, Li et al. (2020) developed an analytical model for radial reactive transport with the skin effect to investigate the impacts of dispersivity, effective porosity, and mass transfer coefficient in the skin zone on radial dispersion. The abovementioned studies demonstrated that the skin effects are significant for radial dispersion.

Thirdly, the solute moves into the formation zone from the skin zone by advective, diffusive, and dispersive processes. Such processes have been widely described by the traditional advection–dispersion equation (ADE) which is based on Fick's law; however, many recent studies demonstrated that the ADE model mainly worked well for homogeneous (or nearly homogeneous) porous media. As for reactive transport in heterogeneous media, the BTCs may exhibit a host of non-Fickian characteristics, such as early arrival and heavy tailing (Di Dato et al., 2017; Molinari et al., 2015). Alternatively, many non-Fickian transport models have been developed, such as the multirate mass transfer model (MRMT) (Le Borgne and Gouze, 2008; Haggerty et al., 2001; Guo et al., 2020), the mobile–immobile model (MIM) (van Genuchten and Wierenga, 1976; Zhou et al., 2017; Wang et al., 2020), continuous-time random-walk (CTRW) models (Dentz et al., 2015; Hansen et al., 2016), fractional-derivative ADE (fADE) models (Soltanpour Moghadam et al., 2022; Chen et al., 2017), a combination of MRMT and CTRW (Kang et al., 2015), and so on (Zheng et al., 2019; Lu et al., 2018). Although the models of MRMT, CTRW, and fADE perform well in modeling non-Fickian transport, it is not easy to obtain the analytical solutions of these models. Meanwhile, these theories are usually not easy to apply when solving regional-scale transport problems, as pointed out in a recent study (Zheng et al., 2019). MIM is an extension of ADE by considering both flowing and stagnant regions in porous media and mass transfer between them (van Genuchten and Wierenga, 1976; Zhou et al., 2017; Wang et al., 2020). Zhou et al. (2017) and Wang et al. (2020) derived the MIM solutions of radial dispersion. However, the skin effect and the scale effect were ignored in their studies, which will be investigated in this study. Besides the MRMT, MIM, CTRW, and fADE models, another approach to representing the heterogeneity is to use a scale-dependent dispersivity (or dispersion) in the ADE or MIM models (Haddad et al., 2015; Gelhar et al., 1992). Gao et al. (2009a) and Chen et al. (2007) discussed radial dispersion and found that the scale-dependent dispersion effect was not negligible. There is also experimental evidence for the scaling of dispersion, mixing, and reaction (Leitão et al., 1996; Edery et al., 2015).

The differences among the currently available analytical solutions for radial dispersion have been reviewed and summarized in Table 1. As one can see from this table, the mixing effect in the wellbore was ignored in all of the models except for Novakowski (1992), Wang et al. (2018, 2020), and Shi et al. (2020). Only Chen et al. (2012), Hsieh and Yeh (2014), Huang et al. (2019), and Li et al. (2020) took the skin effect into account. The differences among the solutions of Tang and Babu (1979), Moench and Ogata (1981), Hsieh (1986), Tang and Peaceman (1987), Yates (1988), Cihan and Tyner (2011), and Chen et al. (2011) mainly consist of the boundary conditions, source-injection types (instantaneous or continuous), and initial conditions.

Table 1 Summary of the current models for the radial dispersion around the recharge well.

Note: “GE”, “ME”, “SCE”, and “SKE” represent governing equation, mixing effect, scale effect, and skin effect, respectively. “Y” and “N” represent whether the effect is considered or not.

In summary, no existing analytical model has ever considered the mixing, skin, scale, and media heterogeneity effects (which are described using MIM) simultaneously. Although the numerical method is more powerful than the analytical method for problems with complex initial and boundary conditions and heterogeneous aquifers of interest, numerical errors could not be avoided easily for the MIM models of concern here, such as numerical dispersion and numerical oscillation issues (Zheng and Wang, 1999; Wang and Zhan, 2013). Meanwhile, the analytical solutions are usually computationally more efficient than the numerical solutions and can be easily coupled into optimization algorithms for problems related to parameter estimation (Neuman and Mishra, 2012). Therefore, a primary purpose of this study is to develop such an analytical model. Furthermore, the accuracy and robustness of the developed model will be tested against a finite-element numerical simulation and experimental data. Moreover, a sensitivity analysis will be conducted to prioritize the influences of various controlling parameters on the newly developed radial dispersion reactive transport model.

2 Methods

2.1 Mathematical model of radial dispersion

An aquifer is assumed to be confined, homogeneous, horizontally isotropic, have a constant thickness, and be fully penetrated by a well from which the solute is injected. A cylindrical coordinate system is established with the r axis horizontal and the z axis vertically upward. The origin of the coordinate system is located at the intersect of the well center and the middle elevation of the aquifer. A schematic diagram of the problem is available in Fig. S1 in the Supplement.

In this study, we mainly focused on developing analytical solutions of radial dispersion with a Heaviside step source (or step function for short hereinafter), as solutions of a variety of injection scenarios can be easily obtained based on such a step source solution, as shown in Eq. (S2), Eqs. (4a, b), or Eqs. (5a, b). Assuming that t_inj is the duration of the step source, the solute source concentration (C₀) is C_inj(t) when time is smaller than t_inj, while it is C_cha(t) when time is greater than t_inj, in which C_inj(t) and C_cha(t) represent the solute concentrations (M L⁻³) in the wellbore before time t_inj and after time t_inj, respectively. When C_cha(t)=0 and t_inj approach zero but the total injected mass remains finite, the model of the step source reduces to the model of the instantaneous injection. Similarly, the model of the step source becomes the model of the continuous injection source when t_inj becomes infinity.

Similar to Chen et al. (2012) and Hsieh and Yeh (2014), a two-region (skin and formation) model of radial dispersion is employed to describe the skin effect. In the skin zone, the governing equations of radial dispersion are

$$\begin{array}{}\text{(1a)}& {\displaystyle }\begin{array}{rl}{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}{R}_{{\mathrm{m}}_{\mathrm{1}}}{\displaystyle \frac{\partial C_{{\mathrm{m}}_{\mathrm{1}}}}{\partial t}}& ={\displaystyle \frac{{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\mathit{\alpha }}_{\mathrm{1}}\left|v_{{\mathrm{a}}_{\mathrm{1}}}\right|{\displaystyle \frac{\partial C_{{\mathrm{m}}_{\mathrm{1}}}}{\partial r}}\right)\\ & -{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}{v}_{{\mathrm{a}}_{\mathrm{1}}}{\displaystyle \frac{\partial C_{{\mathrm{m}}_{\mathrm{1}}}}{\partial r}}-{\mathit{\omega }}_{\mathrm{1}}\left(C_{{\mathrm{m}}_{\mathrm{1}}}-C_{{\mathrm{im}}_{\mathrm{1}}}\right)\\ & -{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}{\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{1}}}{C}_{{\mathrm{m}}_{\mathrm{1}}},r_{\mathrm{w}}\le r\le r_{\mathrm{s}},\end{array}\text{(1b)}& {\displaystyle }\begin{array}{rl}{\mathit{\theta }}_{{\mathrm{im}}_{\mathrm{1}}}{R}_{{\mathrm{im}}_{\mathrm{1}}}{\displaystyle \frac{\partial C_{{\mathrm{im}}_{\mathrm{1}}}}{\partial t}}=& {\mathit{\omega }}_{\mathrm{1}}\left(C_{{\mathrm{m}}_{\mathrm{1}}}-C_{{\mathrm{im}}_{\mathrm{1}}}\right)-{\mathit{\theta }}_{{\mathrm{im}}_{\mathrm{1}}}{\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{1}}}{C}_{{\mathrm{im}}_{\mathrm{1}}},\\ & r_{\mathrm{w}}\le r\le r_{\mathrm{s}}.\end{array}\end{array}$$

In the formation zone, one has

$$\begin{array}{}\text{(1c)}& \begin{array}{rl}{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{2}}}{R}_{{\mathrm{m}}_{\mathrm{2}}}{\displaystyle \frac{\partial C_{{\mathrm{m}}_{\mathrm{2}}}}{\partial t}}& ={\displaystyle \frac{{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{2}}}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\mathit{\alpha }}_{\mathrm{2}}\left|v_{{\mathrm{a}}_{\mathrm{2}}}\right|{\displaystyle \frac{\partial C_{{\mathrm{m}}_{\mathrm{2}}}}{\partial r}}\right)\\ & -{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{2}}}{v}_{{\mathrm{a}}_{\mathrm{2}}}{\displaystyle \frac{\partial C_{{\mathrm{m}}_{\mathrm{2}}}}{\partial r}}-{\mathit{\omega }}_{\mathrm{2}}\left(C_{{\mathrm{m}}_{\mathrm{2}}}-C_{{\mathrm{im}}_{\mathrm{2}}}\right)\\ & -{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{2}}}{\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{2}}}{C}_{{\mathrm{m}}_{\mathrm{2}}},r>r_{\mathrm{s}},\end{array}\end{array}$$

$$\begin{array}{}\text{(1d)}& \begin{array}{rl}{\mathit{\theta }}_{{\mathrm{im}}_{\mathrm{2}}}{R}_{{\mathrm{im}}_{\mathrm{2}}}{\displaystyle \frac{\partial C_{{\mathrm{im}}_{\mathrm{2}}}}{\partial t}}=& {\mathit{\omega }}_{\mathrm{2}}\left(C_{{\mathrm{m}}_{\mathrm{2}}}-C_{{\mathrm{im}}_{\mathrm{2}}}\right)-{\mathit{\theta }}_{{\mathrm{im}}_{\mathrm{2}}}{C}_{{\mathrm{im}}_{\mathrm{2}}}\\ & r>r_{\mathrm{s}},\end{array}\end{array}$$

where the subscripts “m” and “im” refer to parameters in the mobile and immobile domains, respectively. The subscripts “1” and “2” refer to parameters in the skin and formation regions, respectively. $C_{{\mathrm{m}}_{\mathrm{1}}}$ and $C_{{\mathrm{im}}_{\mathrm{1}}}$ are the mobile and immobile concentrations (M L⁻³) of the skin zone, respectively. $C_{{\mathrm{m}}_{\mathrm{2}}}$ and $C_{{\mathrm{im}}_{\mathrm{2}}}$ are the mobile and immobile concentrations (M L⁻³) of the formation zone, respectively. r is the radial distance (L) from the center of the well. r_w is the well radius. r_s is the radial distance (L) from the center of the well to the outer boundary of the skin zone. v_a is the average radial pore velocity (L T⁻¹) in the aquifer. $v_{{\mathrm{a}}_{\mathrm{1}}}=\frac{u_{{\mathrm{a}}_{\mathrm{1}}}}{{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}}$ and $v_{{\mathrm{a}}_{\mathrm{2}}}=\frac{u_{{\mathrm{a}}_{\mathrm{2}}}}{{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{2}}}}$; $u_{{\mathrm{a}}_{\mathrm{1}}}$ and $u_{{\mathrm{a}}_{\mathrm{2}}}$ represent Darcian velocities (L T⁻¹) in the skin and formation zones, respectively. α₁ and α₂ represent the longitudinal dispersivities (L) in the skin and formation zones, respectively. ${\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{1}}}$, ${\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{1}}}$, ${\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{2}}}$, and ${}_{{\mathrm{im}}_{\mathrm{2}}}$ are reaction rates for the first-order reaction rate, the first-order biodegradation, or the radioactive decay (T⁻¹). ${\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}$, ${\mathit{\theta }}_{{\mathrm{im}}_{\mathrm{1}}}$, ${\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{2}}}$, and ${\mathit{\theta }}_{{\mathrm{im}}_{\mathrm{2}}}$ are porosities. $R_{{\mathrm{m}}_{\mathrm{1}}}$, $R_{{\mathrm{im}}_{\mathrm{1}}}$, $R_{{\mathrm{m}}_{\mathrm{2}}}$, and $R_{{\mathrm{im}}_{\mathrm{2}}}$ are retardation factors (dimensionless). ω₁ and ω₂ represent the first-order mass transfer coefficients (T⁻¹) between the mobile and immobile dissolved phases in the skin and formation zones, respectively. One point to note is that the molecular diffusive effect is assumed to be negligible in the above governing equations.

Assuming that the skin and formation zones are initially free of solute, the initial conditions are

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\left.C_{{\mathrm{m}}_{\mathrm{1}}}(r,t)\right|}_{t=\mathrm{0}}& ={\left.C_{{\mathrm{im}}_{\mathrm{1}}}(r,t)\right|}_{t=\mathrm{0}}={\left.C_{{\mathrm{m}}_{\mathrm{2}}}(r,t)\right|}_{t=\mathrm{0}}\\ & ={\left.C_{{\mathrm{im}}_{\mathrm{2}}}(r,t)\right|}_{t=\mathrm{0}}=r\ge r_{\mathrm{w}}.\end{array}\end{array}$$

The outer boundary condition at an infinite distance is

$$\begin{array}{}\text{(3)}& {\left.C_{{\mathrm{m}}_{\mathrm{2}}}(r,t)\right|}_{r\to \mathrm{\infty }}={\left.C_{{\mathrm{im}}_{\mathrm{2}}}(r,t)\right|}_{r\to \mathrm{\infty }}=\mathrm{0}.\end{array}$$

Two types of models have been widely applied to the boundary condition at the well screen: the mass flux continuity (MFC) model and the resident concentration continuity (RCC) model. The RCC model is

$$\begin{array}{}\text{(4a)}& {\displaystyle }& {\displaystyle }{\left.\left[C_{{\mathrm{m}}_{\mathrm{1}}}(r,t)\right]\right|}_{r=r_{\mathrm{w}}}={\left.\left[C_{\mathrm{inj}}\left(t\right)\right]\right|}_{r=r_{\mathrm{w}}},\mathrm{0}<t\le t_{\mathrm{inj}},\text{(4b)}& {\displaystyle }& {\displaystyle }{\left.\left[C_{{\mathrm{m}}_{\mathrm{1}}}(r,t)\right]\right|}_{r=r_{\mathrm{w}}}={\left.\left[C_{\mathrm{cha}}\left(t\right)\right]\right|}_{r=r_{\mathrm{w}}},t>t_{\mathrm{inj}},\end{array}$$

and the MFC model is

$$\begin{array}{}\text{(5a)}& {\displaystyle }\begin{array}{rl}& {\left.\left[C_{{\mathrm{m}}_{\mathrm{1}}}(r,t)-{\mathit{\alpha }}_{\mathrm{1}}{\displaystyle \frac{\left|v_{{\mathrm{a}}_{\mathrm{1}},\mathrm{inj}}\right|}{v_{{\mathrm{a}}_{\mathrm{1}},\mathrm{inj}}}}{\displaystyle \frac{\partial C_{{\mathrm{m}}_{\mathrm{1}}}(r,t)}{\partial r}}\right]\right|}_{r=r_{\mathrm{w}}}\\ & ={\left.\left[C_{\mathrm{inj}}\left(t\right)\right]\right|}_{r=r_{\mathrm{w}}},\mathrm{0}<t\le t_{\mathrm{inj}},\end{array}\text{(5b)}& {\displaystyle }\begin{array}{rl}& {\left.\left[C_{{\mathrm{m}}_{\mathrm{1}}}(r,t)-{\mathit{\alpha }}_{\mathrm{1}}{\displaystyle \frac{\left|v_{{\mathrm{a}}_{\mathrm{1}},\mathrm{cha}}\right|}{v_{{\mathrm{a}}_{\mathrm{1}},\mathrm{cha}}}}{\displaystyle \frac{\partial C_{{\mathrm{m}}_{\mathrm{1}}}(r,t)}{\partial r}}\right]\right|}_{r=r_{\mathrm{w}}}\\ & ={\left.\left[C_{\mathrm{cha}}\left(t\right)\right]\right|}_{r=r_{\mathrm{w}}},t>t_{\mathrm{inj}},\end{array}\end{array}$$

where $v_{{\mathrm{a}}_{\mathrm{1}},\mathrm{inj}}$ and $v_{{\mathrm{a}}_{\mathrm{1}},\mathrm{cha}}$ refer to velocities in the injection and chasing phases, respectively. It was demonstrated that the mass balance requirement could not be satisfied in the RCC model, while the resident concentration was not continuous in the MFC model (Wang et al., 2018). Many experimental studies demonstrated that the MFC model performed better than the RCC model (Novakowski, 1992). Therefore, the MFC model will be used to describe the boundary condition in the wellbore in this study. Comparing Eqs. (4) and (5), one may find that the main difference between these two models is whether the dispersivity is involved or not. Recently, Wang et al. (2019) pointed out that the conflicts between these two models could be resolved by a scale-dependent dispersivity, which was zero at the well screen and increased with the travel distance of the solute. This is because when the dispersivity is zero in Eqs. (5a) and (5b), the MFC model reduces to the RCC model. The model of the scale-dependent dispersivity will be discussed in Sect. 2.4.

When taking into account the mixing effect in the wellbore, one has

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{V}_{\mathrm{w},\mathrm{inj}}{\displaystyle \frac{\mathrm{d}{C}_{\mathrm{inj}}}{\mathrm{d}t}}=-\mathit{\xi }{v}_{{\mathrm{a}}_{\mathrm{1}},\mathrm{inj}}\left(r_{\mathrm{w}}\right)\left[C_{\mathrm{inj}}\left(t\right)-C_{\mathrm{0}}\right],\mathrm{0}<t\le t_{\mathrm{inj}},\text{(7)}& {\displaystyle }& {\displaystyle }{V}_{\mathrm{w},\mathrm{cha}}{\displaystyle \frac{\mathrm{d}{C}_{\mathrm{cha}}}{\mathrm{d}t}}=-\mathit{\xi }{v}_{{\mathrm{a}}_{\mathrm{1}},\mathrm{cha}}\left(r_{\mathrm{w}}\right)\left[C_{\mathrm{cha}}\left(t\right)\right],t>t_{\mathrm{inj}},\end{array}$$

where V_w,inj is the volume (L³) of water in the wellbore when t≤t_inj, and $V_{\mathrm{w},\mathrm{inj}}=\mathit{\pi }{r}_{\mathrm{w}}^{\mathrm{2}}{h}_{\mathrm{w},\mathrm{inj}}$. h_w,inj is the water level (L) in the wellbore when t≤t_inj. $\mathit{\xi }=\mathrm{2}\mathit{\pi }{r}_{\mathrm{w}}{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}B$, B is the thickness (L) of the aquifer, V_w,cha is the volume (L³) of water in the wellbore when t>t_inj, and $V_{\mathrm{w},\mathrm{cha}}=\mathit{\pi }{r}_{\mathrm{w}}^{\mathrm{2}}{h}_{\mathrm{w},\mathrm{cha}}$. h_w,cha is the water level (L) in the wellbore when t>t_inj, $v_{{\mathrm{a}}_{\mathrm{1}},\mathrm{inj}}\left(r_{\mathrm{w}}\right)$ is the velocity at the well screen in the injection phase, and $v_{{\mathrm{a}}_{\mathrm{1}},\mathrm{inj}}\left(r_{\mathrm{w}}\right)=\frac{Q_{\mathrm{inj}}}{\mathrm{2}\mathit{\pi }B{r}_{\mathrm{w}}{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}}$. $v_{{\mathrm{a}}_{\mathrm{1}},\mathrm{cha}}\left(r_{\mathrm{w}}\right)$ is the velocity at the well screen in the chasing phase, and it equals $\frac{Q_{\mathrm{cha}}}{\mathrm{2}\mathit{\pi }B{r}_{\mathrm{w}}{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}}$. Q_inj and Q_cha are the well flow rates (L³ T⁻¹) in the injection and chasing phases, respectively. The mass balance for the well in Eqs. (6) and (7) is only relevant when velocity exceeds zero because it does not contain terms for possible diffusive losses.

The water level in the wellbore (e.g., h_w,inj and h_w,cha) could be determined by solving the groundwater flow model. In the steady state, one has

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }{Q}_{\mathrm{inj}}=\mathrm{2}\mathit{\pi }rBK{\displaystyle \frac{\mathrm{d}h}{\mathrm{d}r}},\mathrm{0}<t\le t_{\mathrm{inj}},\text{(9)}& {\displaystyle }& {\displaystyle }{Q}_{\mathrm{cha}}=\mathrm{2}\mathit{\pi }rBK{\displaystyle \frac{\mathrm{d}h}{\mathrm{d}r}},t>t_{\mathrm{inj}},\end{array}$$

where K is the hydraulic conductivity (L T⁻¹), and $K=\left\{\begin{array}{cl}{K}_{\mathrm{1}}& \mathrm{when}{r}_{\mathrm{w}}\le r<r_{\mathrm{s}}\\ K_{\mathrm{2}}& \mathrm{when}{r}_{\mathrm{s}}\le r\end{array}\right.$; K₁ and K₂ are the hydraulic conductivities (L T⁻¹) of the skin and formation zones, respectively.

By conducting the integration on Eqs. (8) and (9) from r_w to r_s and from r_s to r_e, respectively, the water level in the wellbore can be obtained as follows:

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{h}_{\mathrm{w},\mathrm{inj}}=h_{\mathrm{0}}+{\displaystyle \frac{Q_{\mathrm{inj}}}{\mathrm{2}\mathit{\pi }B{K}_{\mathrm{1}}}}ln{\displaystyle \frac{r_{\mathrm{s}}}{r_{\mathrm{w}}}}+{\displaystyle \frac{Q_{\mathrm{inj}}}{\mathrm{2}\mathit{\pi }B{K}_{\mathrm{2}}}}ln{\displaystyle \frac{r_{\mathrm{e}}}{r_{\mathrm{s}}}},\mathrm{0}<t\le t_{\mathrm{inj}},\text{(11)}& {\displaystyle }& {\displaystyle }{h}_{\mathrm{w},\mathrm{cha}}=h_{\mathrm{0}}+{\displaystyle \frac{Q_{\mathrm{cha}}}{\mathrm{2}\mathit{\pi }B{K}_{\mathrm{1}}}}ln{\displaystyle \frac{r_{\mathrm{s}}}{r_{\mathrm{w}}}}+{\displaystyle \frac{Q_{\mathrm{cha}}}{\mathrm{2}\mathit{\pi }B{K}_{\mathrm{2}}}}ln{\displaystyle \frac{r_{\mathrm{e}}}{r_{\mathrm{s}}}},t>t_{\mathrm{inj}},\end{array}$$

where r_e is the radial distance (L) from the center of the well to the outer boundary of the formation zone, and h₀ is the hydraulic head (L) at r_e. One could find that a finite radius r_e is needed to keep h_w finite. It seems to contradict the boundary condition of the transport problems, which is at infinity, as shown in Eq. (3). The reason for such a “contradiction” could be explained as follows. In reality, the influence area is limited by the finite injection rate and the finite injection time of the well (from a plane-view perspective), bounded by a circle with a radius r_e where the hydraulic head is almost constant and the flow velocity is almost zero.

At the interface between the skin and formation zones, the concentration and dispersive flux have to be continuous, and one has

$$\begin{array}{}\text{(12)}& {\displaystyle }{C}_{{\mathrm{m}}_{\mathrm{1}}}\left(r_{\mathrm{s}},t\right)=C_{{\mathrm{m}}_{\mathrm{2}}}\left(r_{\mathrm{s}},t\right),t>\mathrm{0},\text{(13)}& {\displaystyle }\begin{array}{rl}& {\left.\left[{\mathit{\alpha }}_{\mathrm{1}}\left|v_{{\mathrm{a}}_{\mathrm{1}}}\right|{\displaystyle \frac{\partial C_{{\mathrm{m}}_{\mathrm{1}}}(r,t)}{\partial r}}\right]\right|}_{r=r_{\mathrm{s}}}={\left.\left[{\mathit{\alpha }}_{\mathrm{2}}\left|v_{{\mathrm{a}}_{\mathrm{2}}}\right|{\displaystyle \frac{\partial C_{{\mathrm{m}}_{\mathrm{2}}}(r,t)}{\partial r}}\right]\right|}_{r=r_{\mathrm{s}}},\\ & t>\mathrm{0}.\end{array}\end{array}$$

Since Darcy fluxes (advective solute fluxes) are continuous, dispersive fluxes must be continuous, which is Eq. (13).

Here, it is worthwhile to comment on the nature of using the MIM approach to describe transport in heterogeneous aquifers. First, it has been commonly observed that the aquifer heterogeneity renders the use of ADE invalid in many cases, as ADE is developed and used primarily for homogeneous aquifers. In particular, ADE fails to explain the early breakthrough and long tailing phenomena that are frequently observed in transport in heterogeneous aquifers, as illustrated in the introduction. Second, a striking feature of a heterogeneous aquifer is that a sequence of mobile and less mobile regions coexists. In contrast, a homogeneous aquifer may be simplified as a single (mobile) region. Ideally, suppose one knows precisely the spatial distribution of those mobile and less mobile regions and their associated flow and transport parameters. In that case, one can use a high-resolution numerical simulator to predict the flow and transport process precisely. Unfortunately, this is not feasible for most practical cases. Therefore, as an alternative, we have adopted the concept of the MIM approach in which two continuums consisting of a mobile domain and an immobile domain coexist over the entire heterogeneous aquifer. Each of these two continuums has uniform flow and transport parameters (such as porosity or retardation factor) for the sake of simplification. Furthermore, mass can transfer between these two continuums in a certain fashion, usually using the first-order rate-limited equation. Third, this alternative approach has successfully explained many phenomena that cannot be explained using ADE, for instance, the early breakthrough and long tailing issues. Later, the two-continuum MIM approach was expanded to multiple-continuum MIM or multirate MIM approaches to better capture the transport features in a heterogeneous aquifer (Vangenuchten and Wierenga, 1976; Elenius and Abriola, 2019; Guo et al., 2019). In summary, the MIM does not incorporate the spatial variation of flow and transport parameters that are mostly unknown. Instead, it is based on an alternative approach, using two or more interrelated continuums, and in each continuum, the flow and transport parameters remain uniform over space. To date, the validation of the MIM model has been tested by numerous experimental studies (Griffioen et al., 1998; Gao et al., 2009b; Elenius and Abriola, 2019).

2.2 Solution of radial dispersion

In this study, dimensionless forms of parameters used in the derivation of analytical solutions are defined as $C_{{\mathrm{m}}_{\mathrm{1}}\mathrm{D}}=\frac{C_{{\mathrm{m}}_{\mathrm{1}}}}{C_{\mathrm{0}}}$, $C_{{\mathrm{im}}_{\mathrm{1}}\mathrm{D}}=\frac{C_{{\mathrm{im}}_{\mathrm{1}}}}{C_{\mathrm{0}}}$, $C_{{\mathrm{m}}_{\mathrm{2}}\mathrm{D}}=\frac{C_{{\mathrm{m}}_{\mathrm{2}}}}{C_{\mathrm{0}}}$, $C_{{\mathrm{im}}_{\mathrm{2}}\mathrm{D}}=\frac{C_{{\mathrm{im}}_{\mathrm{2}}}}{C_{\mathrm{0}}}$, $C_{\mathrm{inj},\mathrm{D}}=\frac{C_{\mathrm{inj}}}{C_{\mathrm{0}}}$, $C_{\mathrm{cha},\mathrm{D}}=\frac{C_{\mathrm{cha}}}{C_{\mathrm{0}}}$, $t_{\mathrm{D}}=\frac{\left|A\right|t}{{\mathit{\alpha }}_{\mathrm{2}}^{\mathrm{2}}{R}_{{\mathrm{m}}_{\mathrm{1}}}}$, $t_{\mathrm{inj},\mathrm{D}}=\frac{\left|A\right|t_{\mathrm{inj}}}{{\mathit{\alpha }}_{\mathrm{2}}^{\mathrm{2}}{R}_{{\mathrm{m}}_{\mathrm{1}}}}$, $r_{\mathrm{D}}=\frac{r}{{\mathit{\alpha }}_{\mathrm{2}}}$, $r_{\mathrm{wD}}=\frac{r_{\mathrm{w}}}{{\mathit{\alpha }}_{\mathrm{2}}}$, $r_{\mathrm{sD}}=\frac{r_{\mathrm{s}}}{{\mathit{\alpha }}_{\mathrm{2}}}$, $r_{\mathrm{0}\mathrm{D}}=\frac{r_{\mathrm{0}}}{{\mathit{\alpha }}_{\mathrm{2}}}$, ${\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{1}}\mathrm{D}}=\frac{{\mathit{\alpha }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{1}}}}{A}$, ${\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{1}}\mathrm{D}}=\frac{{\mathit{\alpha }}_{\mathrm{2}}^{\mathrm{2}}{R}_{{\mathrm{m}}_{\mathrm{1}}}{\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{1}}}}{R_{{\mathrm{im}}_{\mathrm{1}}}A}$, ${\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{2}}\mathrm{D}}=\frac{{\mathit{\alpha }}_{\mathrm{2}}^{\mathrm{2}}{\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{2}}}{R}_{{\mathrm{m}}_{\mathrm{1}}}}{A{R}_{{\mathrm{m}}_{\mathrm{2}}}}$, ${\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{2}}\mathrm{D}}=\frac{{\mathit{\alpha }}_{\mathrm{2}}^{\mathrm{2}}{R}_{{\mathrm{m}}_{\mathrm{1}}}{R}_{{\mathrm{im}}_{\mathrm{2}}}}{R_{{\mathrm{im}}_{\mathrm{2}}}A}$, and $A=\frac{Q}{\mathrm{2}\mathit{\pi }B{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}}$.

The detailed derivation of the analytical solution in the Laplace domain can be seen in Sect. S1 in the Supplement. The analytical solution is

$$\begin{array}{}\text{(14a)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{1}}\mathrm{D}}=& N_{\mathrm{1}}\exp \left({\displaystyle \frac{r_{\mathrm{D}}}{\mathrm{2}\mathit{\lambda }}}\right)A_i\left(y_{\mathrm{1}}\right)+N_{\mathrm{2}}\exp \left({\displaystyle \frac{r_{\mathrm{D}}}{\mathrm{2}\mathit{\lambda }}}\right)B_i\left(y_{\mathrm{1}}\right),\\ & r_{\mathrm{wD}}\le r_{\mathrm{D}}\le r_{\mathrm{sD}},\end{array}\text{(14b)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{C}}_{{\mathrm{im}}_{\mathrm{1}}\mathrm{D}}={\displaystyle \frac{{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{1}}}}{s+{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{1}}}+{\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{1}}\mathrm{D}}}}{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{1}}\mathrm{D}},r_{\mathrm{wD}}\le r_{\mathrm{D}}\le r_{\mathrm{sD}},\end{array}\end{array}$$

$$\begin{array}{}\text{(15a)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{2}}\mathrm{D}}=N_{\mathrm{3}}\exp \left({\displaystyle \frac{r_{\mathrm{D}}}{\mathrm{2}}}\right)A_i\left(y_{\mathrm{2}}\right),r_{\mathrm{D}}>r_{\mathrm{sD}},\text{(15b)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{im}}_{\mathrm{2}}\mathrm{D}}={\displaystyle \frac{{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{2}}}}{s+{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{2}}}+{\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{2}}\mathrm{D}}}}{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{2}}\mathrm{D}},r_{\mathrm{D}}>r_{\mathrm{sD}},\end{array}$$

where A_i(⋅) and B_i(⋅) are the Airy functions of the first and second kinds, respectively. C_inj,D and C_cha,D can be determined by Eqs. (S18) and (S19), which can be seen in Sect. S1. $A_i^{\prime }(\cdot )$ and $B_i^{\prime }(\cdot )$ are the derivatives of the Airy function of the first and second kinds, respectively. $\mathit{\lambda }=\frac{{\mathit{\alpha }}_{\mathrm{1}}}{{\mathit{\alpha }}_{\mathrm{2}}}$, $\mathit{\eta }=\frac{{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}{R}_{{\mathrm{m}}_{\mathrm{1}}}}{{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{2}}}{R}_{{\mathrm{m}}_{\mathrm{2}}}}$, $y_{\mathrm{1}}={\left(\frac{E_{\mathrm{1}}}{\mathit{\lambda }}\right)}^{\mathrm{1}/\mathrm{3}}\left(r_{\mathrm{D}}+\frac{\mathrm{1}}{\mathrm{4}\mathit{\lambda }{E}_{\mathrm{1}}}\right)$, $y_{\mathrm{1}\mathrm{s}}={\left(\frac{E_{\mathrm{1}}}{\mathit{\lambda }}\right)}^{\mathrm{1}/\mathrm{3}}\left(r_{\mathrm{sD}}+\frac{\mathrm{1}}{\mathrm{4}\mathit{\lambda }{E}_{\mathrm{1}}}\right)$, $y_{\mathrm{2}}={\left(E_{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{3}}\left(r_{\mathrm{D}}+\frac{\mathrm{1}}{\mathrm{4}{E}_{\mathrm{2}}}\right)$, $y_{\mathrm{2}\mathrm{s}}={\left(E_{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{3}}\left(r_{\mathrm{sD}}+\frac{\mathrm{1}}{\mathrm{4}{E}_{\mathrm{2}}}\right)$, $E_{\mathrm{1}}=s+{\mathit{\epsilon }}_{{\mathrm{m}}_{\mathrm{1}}}+{\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{1}}\mathrm{D}}-\frac{{\mathit{\epsilon }}_{{\mathrm{m}}_{\mathrm{1}}}{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{1}}}}{s+{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{1}}}+{\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{1}}\mathrm{D}}}$, $E_{\mathrm{2}}=\frac{\mathrm{1}}{\mathit{\eta }}\left(s+{\mathit{\epsilon }}_{{\mathrm{m}}_{\mathrm{2}}}+{\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{2}}\mathrm{D}}-\frac{{\mathit{\epsilon }}_{{\mathrm{m}}_{\mathrm{2}}}{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{2}}}}{s+{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{2}}}+{\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{2}}\mathrm{D}}}\right)$, ${\mathit{\beta }}_{\mathrm{inj}}=\frac{V_{\mathrm{w},\mathrm{inj}}{r}_{\mathrm{wD}}}{\mathit{\xi }{R}_{{\mathrm{m}}_{\mathrm{1}}}{\mathit{\alpha }}_{\mathrm{2}}}$, ${\mathit{\beta }}_{\mathrm{cha}}=\frac{V_{\mathrm{w},\mathrm{cha}}{r}_{\mathrm{wD}}}{\mathit{\xi }{R}_{{\mathrm{m}}_{\mathrm{1}}}{\mathit{\alpha }}_{\mathrm{2}}}$, ${\mathit{\epsilon }}_{{\mathrm{m}}_{\mathrm{1}}}=\frac{{\mathit{\omega }}_{\mathrm{1}}{\mathit{\alpha }}_{\mathrm{2}}^{\mathrm{2}}}{A{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}}$, ${\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{1}}}=\frac{{\mathit{\omega }}_{\mathrm{1}}{\mathit{\alpha }}_{\mathrm{2}}^{\mathrm{2}}{R}_{{\mathrm{m}}_{\mathrm{1}}}}{A{\mathit{\theta }}_{{\mathrm{im}}_{\mathrm{1}}}{R}_{{\mathrm{im}}_{\mathrm{1}}}}$, ${\mathit{\epsilon }}_{{\mathrm{m}}_{\mathrm{2}}}=\frac{{\mathit{\omega }}_{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{2}}^{\mathrm{2}}{R}_{{\mathrm{m}}_{\mathrm{1}}}}{A{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{2}}}{R}_{{\mathrm{m}}_{\mathrm{2}}}}$, and ${\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{2}}}=\frac{{\mathit{\omega }}_{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{2}}^{\mathrm{2}}{R}_{{\mathrm{m}}_{\mathrm{1}}}}{A{\mathit{\theta }}_{{\mathrm{im}}_{\mathrm{2}}}{R}_{{\mathrm{im}}_{\mathrm{2}}}}$. s is the dimensionless Laplace transform parameter in respect to dimensionless time t_D. The expressions for N₁, N₂, and N₃ are listed in Table 2.

Table 2 Expressions of coefficients in solutions of Eqs. (14a)–(15b).

From Eqs. (14) and (15), one may find that it is not easy to invert the Laplace-domain solution to obtain the real-time solution analytically. Alternatively, numerical Laplace transform techniques such as the Fourier series method (Dubner and Abate, 1968), Zakian method (Zakian, 1969), Schapery method (Schapery, 1962), de Hoog method (De Hoog et al., 1982), and Stehfest method (Stehfest, 1970) are called in, where the de Hoog and Stehfest methods perform better for problems related to radial dispersion (Wang and Zhan, 2015). In this study, the MATLAB script of the de Hoog method compiled by Hollenbeck (1998) will be employed to facilitate the computation of the inverse Laplace transform, where the numerical tolerance is set to $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{10}}$.

2.3 Special cases of the new solution

The new solution of this study considers the mixing effect, skin effect, and media heterogeneity (which is described using MIM) simultaneously, and the solute is injected into the well as a step source. This general solution encompasses many previous studies as special cases. For instance, when r_s→∞, the skin effect is excluded. t_inj→∞ implies that the solute is continuously injected into the well. t_inj→0 means that the solute is instantaneously injected into the well. ω=0 implies that the MIM solution reduces to the ADE solution. $V_{\mathrm{w},\mathrm{inj}}=\mathrm{0}$ or r_w=0 shows that the mixing effect is excluded.

Therefore, the new solution reduces to the solutions of Hoopes and Harleman (1967), Gelhar and Collins (1971), Tang and Babu (1979), Moench and Ogata (1981), Hsieh (1986), Tang and Peaceman (1987), and Philip (1994) when r_s→∞, t_inj→∞, ω=0, and $V_{w,inj}=\mathrm{0}$. The solution of Wang et al. (2018) is a special case of this study when r_s→∞, t_inj→∞, and ω=0.

2.4 Extension of the new solution with scale-dependent dispersivity

Due to the heterogeneities of the porous media, the dispersivity was found to be dependent on the travel distance of the solute from the source, and such a phenomenon was first observed in the field-scale experiment (Dagan, 1988; Gelhar et al., 1992; Pickens and Grisak, 1981a). The field-scale effect (i.e., dispersivity growing with distance from the well) is usually considered to be a result of spatial heterogeneity at different scales in the aquifer. Subsequently, the scale-dependent dispersivity phenomenon was also found in controlled laboratory tests due to heterogeneities caused by the bridging effect and microstructures, although the sediments (as the porous media) are well sorted and carefully packed (Silliman and Simpson, 1987; Berkowitz et al., 2000; Wang et al., 2019; Gao et al., 2010). For example, Silliman and Simpson (1987) found that the dispersivity continuously increased with distance, based on the experiments conducted in a $\mathrm{2.4}\times \mathrm{1.07}\times \mathrm{0.10}$ m sandbox. Berkowitz et al. (2000) obtained similar conclusions to Silliman and Simpson (1987) in the laboratory-controlled experiment. Wang et al. (2019) also concluded that the scale-dependent model performed better than the scale-independent model in interpreting observed BTCs of the laboratory-controlled experiment. To date, four types of functions have been widely used to describe scale-dependent dispersivity, including asymptotic, parabolic, exponential, and linear functions, as summarized by Pickens and Grisak (1981b). In this section, the model of the scale-independent dispersivity (e.g., Eqs. 14 and 15 in Sect. 2) will be extended by considering the linear-asymptotic dispersivity model in the formation zone. As for the other types of scale-dependent functions, the analytical solutions could be derived using a similar approach. The formula of the linear distance-dependent dispersivity is

$$\begin{array}{}\text{(16)}& {\mathit{\alpha }}_{\mathrm{2}}\left(r\right)=\left\{\begin{array}{ll}kr,& r_{\mathrm{s}}\le r\le r_{\mathrm{0}}\\ {\mathit{\alpha }}_{\mathrm{0}},& r>r_{\mathrm{0}}\end{array}\right.,\end{array}$$

where r₀ is the distance (L), α₂(r₀)=α₀, k is a constant (dimensionless), and the modified solutions are

$$\begin{array}{}\text{(17a)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{1}}\mathrm{D}}=& T_{\mathrm{1}}\exp \left({\displaystyle \frac{r_{\mathrm{D}}}{\mathrm{2}\mathit{\lambda }}}\right)A_i\left(y_{\mathrm{1}}\right)+T_{\mathrm{2}}\exp \left({\displaystyle \frac{r_{\mathrm{D}}}{\mathrm{2}\mathit{\lambda }}}\right)B_i\left(y_{\mathrm{1}}\right),\\ & r_{\mathrm{wD}}\le r_{\mathrm{D}}\le r_{\mathrm{sD}},\end{array}\text{(17b)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{C}}_{{\mathrm{im}}_{\mathrm{1}}\mathrm{D}}={\displaystyle \frac{{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{1}}}}{s+{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{1}}}+{\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{1}}\mathrm{D}}}}{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{1}}\mathrm{D}},r_{\mathrm{wD}}\le r_{\mathrm{D}}\le r_{\mathrm{sD}},\end{array}\end{array}$$

$$\begin{array}{}\text{(18a)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{2}}\mathrm{D}}=T& {}_{\mathrm{3}}{r}_{\mathrm{D}}^m{K}_m\left({\mathit{\epsilon }}_{\mathrm{1}}{r}_{\mathrm{D}}\right)+T_{\mathrm{4}}{r}_{\mathrm{D}}^m{I}_m\left({\mathit{\epsilon }}_{\mathrm{1}}{r}_{\mathrm{D}}\right),\\ & r_{\mathrm{sD}}\le r_{\mathrm{D}}\le r_{\mathrm{0}\mathrm{D}},\end{array}\text{(18b)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{2}}\mathrm{D}}=& T_{\mathrm{5}}\exp \left({\displaystyle \frac{r_{\mathrm{D}}}{\mathrm{2}}}\right)A_i\left(y_{\mathrm{3}}\right)+T_{\mathrm{6}}\exp \left({\displaystyle \frac{r_{\mathrm{D}}}{\mathrm{2}}}\right)B_i\left(y_{\mathrm{3}}\right),\\ & r_{\mathrm{D}}>r_{\mathrm{0}\mathrm{D}},\end{array}\text{(18c)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{C}}_{{\mathrm{im}}_{\mathrm{2}}\mathrm{D}}={\displaystyle \frac{{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{2}}}}{s+{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{2}}}+{\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{2}}\mathrm{D}}}}{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{2}}\mathrm{D}},r_{\mathrm{D}}>r_{\mathrm{sD}},\end{array}\end{array}$$

where $m=\frac{\mathrm{1}}{\mathrm{2}k}$, K_m(⋅) is the mth-order modified Bessel function of the second kind, and I_m(⋅) is the mth-order modified Bessel function of the first kind. The expressions for T₁, T₂, T₃, T₄, T₅, and T₆ are listed in Table 3. $y_{\mathrm{3}}=({\mathit{\epsilon }}_{\mathrm{1}}{)}^{\mathrm{1}/\mathrm{3}}\left(r_{\mathrm{D}}+\frac{\mathrm{1}}{\mathrm{4}{\mathit{\epsilon }}_{\mathrm{1}}}\right)$, $y_{\mathrm{4}}=({\mathit{\epsilon }}_{\mathrm{1}}{)}^{\mathrm{1}/\mathrm{3}}\left(r_{\mathrm{0}\mathrm{D}}+\frac{\mathrm{1}}{\mathrm{4}{\mathit{\epsilon }}_{\mathrm{1}}}\right)$, C_inj,D, and C_cha,D can be determined by Eqs. (S34)–(S36), and the detailed derivation of Eqs. (17) and (18) can be seen in Sect. S2.

Table 3 Expressions of coefficients in solutions of Eqs. (17a)–(18c).

Substituting Eq. (16) into the dispersivity coefficient (D_α), one has

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{D}_{\mathit{\alpha }}& ={\mathit{\alpha }}_{\mathrm{2}}\left|v_{{\mathrm{a}}_{\mathrm{1}}}\right|+D_{\mathrm{0}}\\ & =\left\{\begin{array}{ll}\frac{krQ}{\mathrm{2}\mathit{\pi }rB{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}}+D_{\mathrm{0}}=\frac{kQ}{\mathrm{2}\mathit{\pi }B{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}}+D_{\mathrm{0}},& r_{\mathrm{s}}\le r\le r_{\mathrm{0}}\\ \frac{{\mathit{\alpha }}_{\mathrm{0}}Q}{\mathrm{2}\mathit{\pi }rB{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{2}}}}+D_{\mathrm{0}},& r>r_{\mathrm{0}}\end{array}\right.,\end{array}\end{array}$$

where D₀ is the molecular diffusion coefficient (L² T⁻¹). A few interesting features are notable here. First, because of the unique feature of a divergent flow field in which the velocity is inversely proportional to the radial distance and the use of a dispersivity function that is proportional to the radial distance when r≤r₀, the dispersion coefficient in Eq. (19) actually becomes constant. However, one must be aware that if other types of dispersivity equations are used (such as exponential and parabolic functions), the dispersion coefficient in Eq. (19) will depend on the radial distance from the well. Second, even when D_α becomes constant for a linear dispersivity function when r≤r₀, the mechanical dispersion is still dominant since the value of D₀ is generally much smaller than the mechanical dispersion term of $\frac{kQ}{\mathrm{2}\mathit{\pi }B{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}}$. For instance, the diffusion coefficient in water ranges from $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{9}}$ to $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹, and it is much smaller in the porous media (Freeze and Cherry, 1979). When k=0.01, Q=0.1 m³ s⁻¹, B=1 m, and ${\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}=\mathrm{0.3}$, one has $\frac{kQ}{\mathrm{2}\mathit{\pi }B{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}}=\mathrm{5.3}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹. Therefore, it is reasonable to ignore the molecular diffusion effect when r≤r₀. The values of $\frac{kQ}{\mathrm{2}\mathit{\pi }B{\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}}$ are dependent on k. The chosen value of k=0.01 is from experimental studies, for instance, k=0.018 in Chen et al. (2007) and k=0.024 and 0.013 in this study, as shown in Table 5.

2.5 Extension of the new solution to a leaky confined aquifer

Regardless of Eqs. (14) and (15) or Eqs. (17) and (18), the aquifer is assumed to be completely isolated from the underlying and overlying aquitards (strictly confined), which might not be true in real applications. As stated before (Zhan et al., 2009a, b), it is nearly impossible to maintain a strictly confined condition in terms of transport. This is because as long as a solute in the aquifer is in contact with the upper or lower aquitard, molecular diffusion will always drive the solute from a high-concentration aquifer into the solute-free aquitard, even if the cross-formation flow in the aquitard does not exist. In fact, such diffusion-driven transport of a solute into the aquitard and the subsequent back-diffusion (from aquitard to aquifer when the aquifer solute concentration drops below the solute concentration in the aquitards) is responsible for many long tails in aquifer BTCs. The importance of aquitards in regulating solute transport has indeed been recognized by a number of investigators, such as Chen (1985, 1986, 1991), Yates (1988), Novakowski (1992), Q. Wang and Zhan (2013), and Zhou et al. (2017).

In this section, the solutions of Eqs. (14) and (15) will be extended considering underlying and overlying aquitards. The detailed derivation of the analytical solution in the Laplace domain can be seen in Sect. S3.

In the aquifer, the solutions are

$$\begin{array}{}\text{(20a)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{1}}\mathrm{D}}=& T_{\mathrm{1}}\exp \left({\displaystyle \frac{r_{\mathrm{D}}}{\mathrm{2}\mathit{\lambda }}}\right)A_i\left({\mathit{\phi }}_{\mathrm{1}}\right)+T_{\mathrm{2}}\exp \left({\displaystyle \frac{r_{\mathrm{D}}}{\mathrm{2}}}\right)B_i\left({\mathit{\phi }}_{\mathrm{1}}\right),\\ & r_{\mathrm{wD}}\le r_{\mathrm{D}}\le r_{\mathrm{sD}},\end{array}\text{(20b)}& {\displaystyle }\begin{array}{rl}& {\stackrel{\mathrm{‾}}{C}}_{{\mathrm{im}}_{\mathrm{1}}\mathrm{D}}={\displaystyle \frac{{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{1}}}}{s+{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{1}}}+{\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{1}}\mathrm{D}}}}{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{1}}\mathrm{D}},r_{\mathrm{wD}}\le r_{\mathrm{D}}\le r_{\mathrm{sD}},\end{array}\end{array}$$

$$\begin{array}{}\text{(21a)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{2}}\mathrm{D}}=T_{\mathrm{3}}\exp \left({\displaystyle \frac{r_{\mathrm{D}}}{\mathrm{2}}}\right)A_i\left({\mathit{\phi }}_{\mathrm{2}}\right),r_{\mathrm{D}}>r_{\mathrm{sD}},\text{(21b)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{im}}_{\mathrm{2}}\mathrm{D}}={\displaystyle \frac{{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{2}}}}{s+{\mathit{\epsilon }}_{{\mathrm{im}}_{\mathrm{2}}}+{\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{2}}\mathrm{D}}}}{\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{2}}\mathrm{D}},r_{\mathrm{D}}>r_{\mathrm{sD}}.\end{array}$$

In the aquitards, the solutions are

$$\begin{array}{}\text{(22a)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{\mathrm{umD}}={\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{1}}\mathrm{D}}\exp \left(a_{\mathrm{2}}{z}_{\mathrm{D}}-a_{\mathrm{2}}\right),r_{\mathrm{wD}}\le r_{\mathrm{D}}\le r_{\mathrm{sD}},\text{(22b)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{\mathrm{umD}}={\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{2}}\mathrm{D}}\exp \left(a_{\mathrm{2}}{z}_{\mathrm{D}}-a_{\mathrm{2}}\right),r_{\mathrm{D}}>r_{\mathrm{sD}},\text{(22c)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{\mathrm{uimD}}={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{uim}}}{s+{\mathit{\epsilon }}_{\mathrm{uim}}+{\mathit{\mu }}_{\mathrm{uimD}}}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{umD}},r_{\mathrm{D}}>r_{\mathrm{wD}},\end{array}$$

$$\begin{array}{}\text{(23a)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{\mathrm{lmD}}={\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{1}}\mathrm{D}}\exp \left(b_{\mathrm{1}}{z}_{\mathrm{D}}+b_{\mathrm{1}}\right),r_{\mathrm{wD}}\le r_{\mathrm{D}}\le r_{\mathrm{sD}},\text{(23b)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{\mathrm{lmD}}={\stackrel{\mathrm{‾}}{C}}_{{\mathrm{m}}_{\mathrm{2}}\mathrm{D}}\exp \left(b_{\mathrm{1}}{z}_{\mathrm{D}}+b_{\mathrm{1}}\right),r_{\mathrm{D}}>r_{\mathrm{sD}},\text{(23c)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{‾}}{C}}_{\mathrm{limD}}={\displaystyle \frac{{\mathit{\epsilon }}_{\lim }}{s+{\mathit{\epsilon }}_{\lim }+{\mathit{\mu }}_{\mathrm{limD}}}}{\stackrel{\mathrm{‾}}{C}}_{\mathrm{lmD}},r_{\mathrm{D}}>r_{\mathrm{wD}},\end{array}$$

where the letters “u” and “l” in the subscripts represent the upper and lower aquitards, respectively, ${\mathit{\phi }}_{\mathrm{w}}={\left(\frac{E_{\mathrm{3}}}{\mathit{\lambda }}\right)}^{\mathrm{1}/\mathrm{3}}\left(r_{\mathrm{wD}}+\frac{\mathrm{1}}{\mathrm{4}\mathit{\lambda }{E}_{\mathrm{3}}}\right)$, ${\mathit{\phi }}_{\mathrm{1}}={\left(\frac{E_{\mathrm{3}}}{\mathit{\lambda }}\right)}^{\mathrm{1}/\mathrm{3}}\left(r_{\mathrm{D}}+\frac{\mathrm{1}}{\mathrm{4}\mathit{\lambda }{E}_{\mathrm{3}}}\right)$, ${\mathit{\phi }}_{\mathrm{2}}=E_{\mathrm{4}}^{\mathrm{1}/\mathrm{3}}\left(r_{\mathrm{D}}+\frac{\mathrm{1}}{\mathrm{4}{E}_{\mathrm{4}}}\right)$, ${\mathit{\phi }}_{\mathrm{1}\mathrm{s}}={\left(\frac{E_{\mathrm{3}}}{\mathit{\lambda }}\right)}^{\mathrm{1}/\mathrm{3}}\left(r_{\mathrm{sD}}+\frac{\mathrm{1}}{\mathrm{4}\mathit{\lambda }{E}_{\mathrm{3}}}\right)$, and ${\mathit{\phi }}_{\mathrm{2}\mathrm{s}}=E_{\mathrm{4}}^{\mathrm{1}/\mathrm{3}}\left(r_{\mathrm{sD}}+\frac{\mathrm{1}}{\mathrm{4}{E}_{\mathrm{4}}}\right)$. The expressions for a₂, b₁, T₁, T₂, and T₃ are listed in Table 4. C_inj,D and C_cha,D can be determined by Eqs. (S71) and (S72), which can be seen in Sect. S3.

Table 4 Expressions of coefficients in solutions of Eqs. (20a)–(23c).

The solutions of Chen (1985, 1986, 1991) and Yates (1988) are special cases of this study when r_s→∞, t_inj→∞, ω=0, and $V_{\mathrm{w},\mathrm{inj}}=\mathrm{0}$. When r_s→∞, t_inj→∞, and $V_{\mathrm{w},\mathrm{inj}}=\mathrm{0}$, the new solution reduces to the solution of Zhou et al. (2017). Novakowski (1992) considered the wellbore mixing effect in an aquifer–aquitard system, while he ignored other factors such as the skin effect, scale-dependent dispersivity, and mass transfer between the mobile and immobile domains in porous media.

3 Results and discussion

3.1 Test of new solutions

To test the new solution of this study, a numerical simulation based on the Galerkin finite-element method is conducted in the COMSOL Multiphysics platform. More details about the numerical simulation setup can be seen in Sect. S4. The parameters used in the numerical simulation are r_w=2.5 cm, r_s=12.5 cm, $Q_{\mathrm{inj}}=Q_{\mathrm{cha}}=\mathrm{100}$ mL s⁻¹, t_inj=300 s, α₁=2.5 cm, α₂=2.5 cm, θ_m=0.30, θ_im=0.01, ω=0.001 d⁻¹, $R_{{\mathrm{m}}_{\mathrm{1}}}=R_{{\mathrm{im}}_{\mathrm{1}}}=R_{{\mathrm{m}}_{\mathrm{2}}}=R_{{\mathrm{im}}_{\mathrm{2}}}=\mathrm{1}$, B=50 cm, ${\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{1}}}={\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{2}}}={\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{1}}}{=}_{{\mathrm{im}}_{\mathrm{2}}}={\mathrm{10}}^{-\mathrm{7}}$ s⁻¹, and $h_{\mathrm{w},\mathrm{inj}}=h_{\mathrm{w},\mathrm{cha}}=B$. These parameters are from the experimental applications of Chao (1999), Chen et al. (2017), and Wang et al. (2018, 2020), in which Wang et al. (2020) summarized the values of reaction rate, retardation factor, dispersivity, porosity, and first-order mass transfer coefficient for sand and clay used in numerous investigations, as shown in Table 4 of Wang et al. (2020). In addition, the values of the retardation factor and reaction rate show that the chemical reaction and sorption are weak for the tracer of potassium bromide (KBr) in the experiment of Chao (1999). This is unsurprising since KBr is commonly treated as a “conservative” tracer.

As it is difficult to describe the wellbore mixing effect in COMSOL Multiphysics, the wellbore concentration is computed by the analytical solutions of Eqs. (14) and (15). Figure 1a and b show the comparison of concentration between the numerical and analytical solutions of this study, and good agreement between these two kinds of solutions is evident for different times and locations. The comparisons between the numerical solution and analytical solutions of Eqs. (20)–(23) are shown in Sect. S4.2, and the agreement is also good between them.

Figure 1 Comparison of the numerical solution by COMSOL Multiphysics and the analytical solution of this study for different times. (a) In the injection phase; (b) in the chasing phase.

3.2 Test of the model using experimental data

To test the influence of the mixing effect, skin effect, scale effect, and heterogeneity of the media on radial dispersion, the experimental data of Chao (1999) are employed. Chao (1999) reported a laboratory experiment of radial dispersion in a sand tank 244 cm in length, 122 cm in width, and 6.35 cm in depth. A well with a radius of 1.0 cm fully penetrated a confined aquifer. Two observation wells were, respectively, located 22.5 cm and 30.4 cm away from the well center. KBr is chosen as a conservative tracer. Before the tracer is introduced into the wellbore, a steady-state flow field is produced by injecting KBr-free water into the aquifer with a constant injection rate of 9.9 mL min⁻¹. The injection time is 5 h (t_inj=300 min) for the tracer while maintaining the same injection rate of 9.9 mL min⁻¹. The experimental data of Chao (1999) were interpreted by Gao et al. (2009a) using the model of Chen et al. (2007), as shown in Fig. 2. “SDM” and “CDM” in the legend of Fig. 2 refer to the scale-dependent dispersivity model and the constant dispersivity model, respectively. Chen et al. (2007) approximated the injection as an instantaneous source (the validity of such a treatment will be addressed later), and the mass M of the instantaneous injection is calculated by

$$\begin{array}{}\text{(24)}& M=C_{\mathrm{0}}{Q}_{\mathrm{inj}}{t}_{\mathrm{inj}}.\end{array}$$

The other parameters of the analytical solution are listed in Table 5. The parameters estimated by Gao et al. (2009a) are also included in Table 5 for comparison. One may find that the goodness of fit between the observed data and the models of Gao et al. (2009a) and Chen et al. (2007) seems good at the observation point close to the well, but they could not capture BTCs at r=30.4 cm. This is probably for the following two reasons. Firstly, the model of Chen et al. (2007) used to best fit the data is an instantaneous slug test model, which is a rather gross approximation of the injection, which lasted about 5 h. A more proper way is to treat the 5 h injection as a step source. Secondly, the solution of Chen et al. (2007) only considered the scale-dependent dispersivity but ignored the mixing effect and the mass transfer between the mobile and immobile domains.

Figure 2 Fitness of observed BTC by the solution of Chen et al. (2007) which considers the scale effect but ignores the mixing and skin effects.

Table 5 Parameter values used in Figs. 2 and 3.

Note: “SDM” represents the scale-dependent dispersivity model. “CDM” represents the constant dispersivity model. “–” in units represents that the variable is dimensionless. “ – ” represents that the variable is not included in the model.

To test the new solutions of this study, we try to best-fit the observed data again using the newly developed model considering the scale-dependent dispersivity, mixing effect, and heterogeneity of the media (described using MIM). As there is no aquitard in the controlled laboratory experiment, the aquitard effect is irrelevant. Meanwhile, as there is no skin, the skin effect is not included either. The best fitness between the analytical solution and the experimental data is an optimization process by minimizing the “error” between them:

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

where C_OBS and C_COM represent the observed and computed concentrations, respectively, Er is the error, and N is the number of sampling points. In this study, the genetic algorithm (GA) is employed to search the optimal parameter values, such as ${\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{2}}}$, α₁, and ω₁ for CDM of Eqs. (14) and (15) and ${\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{2}}}$, α₀, k, and ω₁ for SDM of Eqs. (17) and (18). GA is a stochastic search method based on natural selection and is preferred for optimization. Meanwhile, GA has been packaged into the MATLAB toolbox (Katoch et al., 2020; Whitley, 1994; Deb et al., 2002), and therefore it is efficient, simple to program, and robust. The estimated values of some key parameters are listed in Table 5. The errors between the observed and computed BTCs by different models are listed in Table 6. Figure 3 shows the fitness between the analytical solution and experimental data, with and without scale-dependent dispersivity, respectively. As GA converges after 500 generations (iterations), the fitness is good, as shown in Fig. 3, and the estimated parameters are physically sound.

Figure 3 Fitness of observed BTCs by new solutions of this study, where Eqs. (14) and (15) are without the scale effect and Eqs. (17) and (18) are with the scale effect, respectively.

Table 6 Errors between observed and computed BTCs in Figs. 2 and 3.

Comparing Figs. 2 and 3 shows that the solutions of this study perform better than the model of Chen et al. (2007), since the fitness is good for both observation locations. To better evaluate the overall performance of the models for both locations, we have used Eq. (25) and the coefficient of determination (R²) to compute the errors of best fitness with two BTCs simultaneously in Figs. 2 and 3, and R² is defined by

$$\begin{array}{}\text{(26)}& R^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^N{\left(C_{\mathrm{OBS}}-C_{\mathrm{COM}}\right)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^N{\left(C_{\mathrm{OBS}}-{\stackrel{\mathrm{‾}}{C}}_{\mathrm{OBS}}\right)}^{\mathrm{2}}}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{C}}_{\mathrm{OBS}}$ is the average concentration of observed data. The results are listed in the last two columns of Table 6. This table shows that the new model performs better. For example, when using CDM, the overall errors for the best-fitting two locations are 0.89 (which is the summation of 0.06 and 0.83 in Table 6) for Chen et al. (2007) and 0.39 (which is the summation of 0.34 and 0.05 in Table 6) for this study. When using SDM, however, the overall errors for the best-fitting two locations are 0.78 for Chen et al. (2007) and 0.25 for this study. The overall R² shows the same observation as the overall Er, where the overall R² is closer to 2, implying that the model precision is much higher. Evidently, the model with the scale effect is the best choice for interpreting the experimental data.

We must emphasize that better fitting of one model than the other with the experimental data should not be used as the only evidence of proof of model performance. This is because a model with more fitting parameters usually performs better than the model with fewer fitting parameters. Besides the best-fitting exercises, however, one should pay more attention to see whether the model adequately acknowledges the underlying physiochemical principles governing the transport processes. As far as we can see, the new model proposed in this study has honored the underlying physiochemical principles governing the radial dispersion process properly. In addition, the model performance (as reflected in the best-fitting practice with the experimental data) is also considerably better. Therefore, based on these two considerations, the new model of this study can be regarded as a significant advancement of present knowledge on radial dispersion. Furthermore, the new model is quite general and encompasses almost all the existing models as subsets.

3.3 Sensitivity analysis

As the new model involves several controlling parameters, it is necessary to prioritize the importance of these parameters in their control on the model performance. In this study, a sensitivity analysis involving normalized parameters is conducted as follows (Kabala, 2001; Yang and Yeh, 2009):

$$\begin{array}{}\text{(27)}& {\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. I_j represents any one parameter of interest, like the volume of water in the wellbore (V_w), k, ${\mathit{\theta }}_{\mathrm{m}}={\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}={\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{2}}}$, $\mathit{\omega }={\mathit{\omega }}_{\mathrm{1}}={\mathit{\omega }}_{\mathrm{2}}$, and so on. A larger |SC_i,j| value means a higher sensitivity.

As the expression of the new analytical solution is complex, it is not easy to get the values of SC_i,j directly from Eq. (27). Therefore, a finite-difference scheme is used alternatively to approximate the right-hand-side term (Kabala, 2001; Yang and Yeh, 2009):

$$\begin{array}{}\text{(28)}& {\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 of I_j.

Figure 4 SC_i,j of the parameters r_w, r_s, k, θ_m, ω, R, and μ using the parameters estimated by best-fitting the experimental data. (a) r=22.5 cm; (b) r=30.4 cm.

The main parameters of the new model include the volume of the water in the wellbore (V_w) for the mixing effect, r_s and α_s for the skin zone, θ_m and ω for the MIM model, and k for scale-dependent dispersivity. Figure 4a and b show SC_i,j at r=22.5 and r=30.4 cm, respectively. The parameters used in these two figures are the same as those used in Fig. 3.

Two observations can be found from Fig. 4a and b. Firstly, the results are sensitive to the parameter of θ_m. To test such a finding, we use the model of this study (Eqs. 17 and 18) with the mixing effect to best fit the experimental data of Chao (1999) (shown in Fig. S5 in Sect. S5), and the results show that the influence of the mixing effect could be negligible. Secondly, by comparing Fig. 4a and b, we find that the sensitivity coefficient of V_w, r_s, α_s, and R on BTCs increases with the distance from the wellbore.

Figure 4a and b show that the sensitivity coefficient of V_w on BTCs is minimal, which might contradict the finding reported in some previous studies (Wang et al., 2018). A careful inspection indicates that the well radius and the initial water level in the wellbore are tiny in the experiment of Chao (1999), resulting in a minimal value of V_w. From Eqs. (6) and (7), one can see that V_w could be influenced by the pumping rate, well radius, initial water level (h₀), and hydraulic parameters of the aquifer. In actual field practices, the value of V_w can be significantly larger than what Chao (1999) uses. Therefore, the sensitivity coefficient of V_w on BTCs will be investigated again using the well radius and the initial water level that is more commonly seen in field applications, e.g., r_w=5.0 cm and h₀=31.75 cm, and the other parameters are the same as the ones in Fig. 3.

The sensitivity analysis after such modification shows that the parameter with the highest sensitivity coefficient is still θ_m. However, the parameter with the second-highest sensitivity coefficient becomes the volume of water in the wellbore (Fig. 5). Figures 6 and 7 illustrate SC_i,j of V_w for different r_w and different observation locations and that the sensitivity coefficient of V_w increases with the well radius but decreases with the distance from the wellbore.

Figure 5 SC_i,j of the parameters V_w, r_s, k, θ_m, ω, R, and μ when increasing V_w.

Figure 6 SC_i,j of V_w for different r_w at r=22.5 cm.

Figure 7 SC_i,j of V_w for different observed locations when r_w=5.0 cm.

4 Conclusions

Radial dispersion is an important process in the fields of chemical engineering, environmental science, and hydrogeology. It has been commonly employed to describe the reactive transport in the subsurface or to estimate aquifer transport parameters (dispersivity, porosity, reactive rate, etc.) required in optimization of remediation strategies. However, previous studies did not include all of the mixing effect, skin effect, and mass transfer between the mobile and immobile domains in porous media.

In this study, a new general model is developed considering all the abovementioned factors. The new general model is compared against a finite-element numerical model and existing experimental data. Meanwhile, the new model is also expanded considering the effect of the overlying and underlying aquitards and the scale-dependent dispersivity. The sensitivity analysis is conducted to prioritize influences of various controlling parameters on BTCs. The following conclusions could be summarized.

1. The new general model honors the most relevant processes involved in radial dispersion (wellbore mixing effect, well skin effect, aquitard effect, and mass transfer between the mobile and immobile domains), for which a solution has not yet been presented.

2. The new general model fits the experimental data of Chao (1999) much better than previous models.

3. The results are sensitive to parameters θ_m (mobile porosity) and V_w (the volume of water in the wellbore). When V_w is very small, as in the laboratory experiment of Chao (1999), the sensitivity coefficient approaches 0. However, for typical values of V_w in actual field applications, the sensitivity coefficient of V_w increases significantly, and the value is often ranked as the second highest, after that of θ_m.

4. The sensitivity coefficient of V_w increases with well radius, while it decreases with increasing distance from the wellbore.

Appendix A: Nomenclature

Symbol	Description
Ai(⋅), Bi(⋅)	Airy functions of the first kind and
	second kind, respectively
$A_i^{\prime }(\cdot )$, $B_i^{\prime }(\cdot )$	Derivative of the Airy functions of the
	first kind and second kind, respectively
α0	Longitudinal dispersivity (L) in the
	formation zone at r>r0
α1, α2	Longitudinal dispersivities (L) in the skin
	and formation zones, respectively
B	The thickness (L) of the aquifer
b	Half of the aquifer thickness (L)

$C_{{\mathrm{m}}_{\mathrm{1}}}$, $C_{{\mathrm{im}}_{\mathrm{1}}}$	Resident mobile and immobile
	concentrations (M L−3) of
	the skin zone, respectively
$C_{{\mathrm{m}}_{\mathrm{2}}}$, $C_{{\mathrm{im}}_{\mathrm{2}}}$	Resident mobile and immobile
	concentrations (M L−3) of the formation
	zone, respectively
Cum, Cuim	Resident mobile and immobile
	concentrations (M L−3) of the upper
	aquitard, respectively
Clm, Clim	Resident mobile and immobile
	concentrations (M L−3) of the lower
	aquitard, respectively
Cinj(t), Ccha(t)	Concentrations (M L−3) of the tracer
	in the wellbore in the injection and
	chasing phases, respectively
C0	Concentration (M L−3) of the tracer
	injected into the wellbore
Cw	Concentration (M L−3) of the tracer
	in the wellbore
Du, Dl	Vertical dispersion coefficients
	(L2 T−1) of the upper and
	lower aquitards, respectively
D0	Molecular diffusion coefficient
	(L2 T−1)
h	Hydraulic head (L)
h0	Hydraulic head (L) at re
hw,inj, hw,cha	Water level in the wellbore in the
	injection and chasing phases (L)
k	A constant (dimensionless)
	ranging from 0 to 1
K1, K2	Hydraulic conductivities (L T−1) of
	the skin and formation zones,
	respectively
Kd	Equilibrium distribution coefficient
	(M−1 L3) for the linear sorption
	process
Im(⋅), Km(⋅)	The mth-order modified Bessel
	function of the first and second kinds,
	respectively
Q	Pumping rate (L3 T−1) (negative for
	injection and positive for pumping)
Qinj, Qcha	Well flow rates (L3 T−1) in the
	injection and chasing phases,
	respectively.
r	Radial distance (L) from the center
	of the well
rs	Distance (L) from the center of the
	well to the outer boundary of the
	skin zone
rw	Radius (L) of the well
re	Radial distance (L) from the center of
	the well to the outer boundary
r0	Radial distance (L) for the linear

	distance-dependent dispersivity
$R_{{\mathrm{m}}_{\mathrm{1}}}$, $R_{{\mathrm{im}}_{\mathrm{1}}}$	Retardation factors (dimensionless)
	for the mobile and immobile regions
	of the skin zone
$R_{{\mathrm{m}}_{\mathrm{2}}}$, $R_{{\mathrm{im}}_{\mathrm{2}}}$	Retardation factors (dimensionless)
	for the mobile and immobile regions
	of the formation zone
Rum, Ruim	Retardation factors (dimensionless)
	for the mobile and immobile regions
	of the upper aquitard
Rlm, Rlim	Retardation factors (dimensionless)
	for the mobile and immobile regions
	of the lower aquitard
t	Time (T)
tinj, tcha	Ending times (T) of the injection and
	chasing phases, respectively
$v_{{\mathrm{a}}_{\mathrm{1}}}$, $v_{{\mathrm{a}}_{\mathrm{2}}}$	Average radial pore velocities (L T−1)
	of the skin zone and formation zone,
	respectively
$v_{{\mathrm{a}}_{\mathrm{1}},\mathrm{inj}}$, $v_{{\mathrm{a}}_{\mathrm{1}},\mathrm{cha}}$	Average radial pore velocities (L T−1)
	at the well screen in the injection
	and chasing phases, respectively.
vum, vlm	Vertical velocities (L T−1) of the upper
	and lower aquitards, respectively
αu, αl	Dispersivities (L) of the upper
	and lower aquitards, respectively
${\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{1}}}$, ${\mathit{\mu }}_{{\mathrm{im}}_{\mathrm{1}}}$	Decay constant for radioactive decay
	or reaction rate coefficient (T−1) in the
	mobile and immobile regions of the skin
	zone
${\mathit{\mu }}_{{\mathrm{m}}_{\mathrm{2}}}$, ${}_{{\mathrm{im}}_{\mathrm{2}}}$	Decay constant (T−1) for radioactive
	decay or reaction rate coefficient in the
	mobile and immobile regions of the
	formation zone
μum, μuim	Decay constant (T−1) for radioactive
	decay or reaction rate coefficient in the
	mobile and immobile regions of the
	upper aquitard
μlm, μlim	Decay constant (T−1) for radioactive
	decay or reaction rate coefficient in the
	mobile and immobile regions of the
	lower aquitard
${\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{1}}}$, ${\mathit{\theta }}_{{\mathrm{im}}_{\mathrm{1}}}$	Mobile and immobile porosities
	(dimensionless) in the skin zone
${\mathit{\theta }}_{{\mathrm{m}}_{\mathrm{2}}}$, ${\mathit{\theta }}_{{\mathrm{im}}_{\mathrm{2}}}$	Mobile and immobile porosities
	(dimensionless) in the formation zone
θum, θuim	Mobile and immobile porosities
	(dimensionless) in the upper aquitard
θlm, θlim	Mobile and immobile porosities
	(dimensionless) in the lower aquitard
ρb	Bulk density (M L−3) of the aquifer
	material
ω1, ω2	First-order mass transfer coefficients
	(T−1) in the skin and formation zones,
	respectively

s	Laplace transform variable with
	respect to the time tD
Subscript	Description
D	Dimensionless form
m, im	Mobile and immobile regions,
	respectively
inj, cha	Injection and chasing phases,
	respectively
u, l	Upper and lower aquitard,
	respectively
1, 2	Parameters in the skin and
	formation regions,
	respectively
Abbreviation	Description
ADE	Advection–dispersion equation
BTCs	The observed breakthrough curves
CDM	The constant dispersivity model
CTRW	Continuous-time random-walk models
fADE	Fractional-derivative ADE models
GA	The genetic algorithm
MFC	The mass flux continuity
MIM	Mobile–immobile model
MRMT	The multirate mass transfer model
RCC	The resident concentration continuity
SDM	The scale-dependent dispersivity model