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 ﬁrstly experienced a mixing process in the injection wellbore, then entered a skin zone after leaving the injection wellbore, and ﬁnally 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 ﬁnite-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 inﬂuences 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 coefﬁcient 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.

ternative, 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.
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., 2019Li et al., , 2020Huang 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 ini-tially 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., 2019Li et al., , 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 randomwalk (CTRW) models Hansen et al., 2016), fractional-derivative ADE (fADE) models (Soltan-pour Moghadam et al., 2022;Chen et al., 2017), a combination of MRMT and CTRW , 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 scaledependent dispersivity (or dispersion) in the ADE or MIM models (Haddad et al., 2015;Gelhar et al., 1992).  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. (2018Wang et al. ( , 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 andTyner (2011), andChen et al. (2011) mainly consist of the boundary conditions, source-injection types (instantaneous or continuous), and initial conditions.
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 prior-itize the influences of various controlling parameters on the newly developed radial dispersion reactive transport model.

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 0 ) 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 −3 ) 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 In the formation zone, one has  Hsieh (1986), Tang and Peaceman (1987), Yates (1988), Cihan and Tyner (2011), Chen et al. (2011) Chen (1985, 1991 Leaky 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. 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 m 1 and C im 1 are the mobile and immobile concentrations (M L −3 ) of the skin zone, respectively. C m 2 and C im 2 are the mobile and immobile concentrations (M L −3 ) 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 −1 ) in the aquifer. v a 1 = u a 1 θ m 1 and v a 2 = u a 2 θ m 2 ; u a 1 and u a 2 represent Darcian velocities (L T −1 ) in the skin and formation zones, respectively. α 1 and α 2 represent the longitudinal dispersivities (L) in the skin and formation zones, respectively. µ m 1 , µ im 1 , µ m 2 , and im 2 are reaction rates for the first-order reaction rate, the first-order biodegradation, or the radioactive decay (T −1 ). θ m 1 , θ im 1 , θ m 2 , and θ im 2 are porosities. R m 1 , R im 1 , R m 2 , and R im 2 are retardation factors (dimensionless). ω 1 and ω 2 represent the first-order mass transfer coefficients (T −1 ) 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 The outer boundary condition at an infinite distance is 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 and the MFC model is where v a 1 ,inj and v a 1 ,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 where V w,inj is the volume (L 3 ) of water in the wellbore when t ≤ t inj , and V w,inj = π r 2 w h w,inj . h w,inj is the water level (L) in the wellbore when t ≤ t inj . ξ = 2π r w θ m 1 B, B is the thickness (L) of the aquifer, V w,cha is the volume (L 3 ) of water in the wellbore when t > t inj , and V w,cha = π r 2 w h w,cha . h w,cha is the water level (L) in the wellbore when t > t inj , v a 1 ,inj (r w ) is the velocity at the well screen in the injection phase, and v a 1 ,inj (r w ) = Q inj 2πBr w θ m 1 . v a 1 ,cha (r w ) is the velocity at the well screen in the chasing phase, and it equals Q cha 2πBr w θ m 1 . Q inj and Q cha are the well flow rates (L 3 T −1 ) 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 where K is the hydraulic conductivity (L T −1 ), and K = K 1 when r w ≤ r < r s K 2 when r s ≤ r ; K 1 and K 2 are the hydraulic conductivities (L T −1 ) 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: where r e is the radial distance (L) from the center of the well to the outer boundary of the formation zone, and h 0 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 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).

Solution of radial dispersion
In this study, dimensionless forms of parameters used in the derivation of analytical solutions are defined as C m 1 D = The detailed derivation of the analytical solution in the Laplace domain can be seen in Sect. S1 in the Supplement. The analytical solution is 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 (·) and B i (·) are the derivatives of the Airy function of the first and second kinds, respectively.
, and ε im 2 = ω 2 α 2 2 R m 1 Aθ im 2 R im 2 . s is the dimensionless Laplace transform parameter in respect to dimensionless time t D . The expressions for N 1 , N 2 , and N 3 are listed in Table 2.
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 1 × 10 −10 .

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 w,inj = 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 = 0. The solution of Wang et al. (2018) is a special case of this study when r s → ∞, t inj → ∞, and ω = 0.

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 scaledependent 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 2.4 × 1.07 × 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 scaledependent dispersivity, including asymptotic, parabolic, exponential, and linear functions, as summarized by Pickens and Grisak (1981b). In this section, the model of the scaleindependent 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 where r 0 is the distance (L), α 2 (r 0 ) = α 0 , k is a constant (dimensionless), and the modified solutions are where m = 1 2k , 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 1 , T 2 , T 3 , T 4 , T 5 , and T 6 are listed in Table 3. y 3 = (ε 1 ) 1/3 r D + 1 4ε 1 , y 4 = (ε 1 ) 1/3 r 0D + 1 4ε 1 , Table 3. Expressions of coefficients in solutions of Eqs. (17a)-(18c).
W 10 k mr m−1 sD I m (ε 1 r D ) + 0.5ε 1 r m sD I m−1 (ε 1 r D ) + I m+1 (ε 1 r D ) W 11 −kε 1 r m+2 0D K m−1 (ε 1 r 0D ) W 12 k mr m 0D I m (ε 1 r 0D ) + 0.5ε 1 r m+1 0D I m−1 (ε 1 r 0D ) + I m+1 (ε 1 r 0D ) W 13 0.5 exp rD 2 A i (y 4 ) + ε 1/3 1 exp rD 2 A i (y 4 ) W 14 r m 0D K m (ε 1 r 0D ) W 15 r m 0D I m (ε 1 r 0D ) W 16 exp r 0D 2 A i (y 4 ) 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. Substituting Eq. (16) into the dispersivity coefficient (D α ), one has where D 0 is the molecular diffusion coefficient (L 2 T −1 ). 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 0 , 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 0 , the mechanical dispersion is still dominant since the value of D 0 is generally much smaller than the mechanical dispersion term of kQ 2πBθ m 1 . For instance, the diffusion coefficient in water ranges from 1 × 10 −9 to 2 × 10 −9 m 2 s −1 , and it is much smaller in the porous media (Freeze and Cherry, 1979). When k = 0.01, Q = 0.1 m 3 s −1 , B =1 m, and θ m 1 = 0.3, one has kQ 2πBθ m 1 = 5.3 × 10 −4 m 2 s −1 . Therefore, it is reasonable to ignore the molecular diffusion effect when r ≤ r 0 . The values of kQ 2πBθ m 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.

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 , 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 solutefree 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 (1985Chen ( , 1986Chen ( , 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 In the aquitards, the solutions are where the letters "u" and "l" in the subscripts represent the upper and lower aquitards, respectively, ϕ w = r sD + 1 4λE 3 , and ϕ 2s = E 1/3 4 r sD + 1 4E 4 . The expressions for a 2 , b 1 , T 1 , T 2 , and T 3 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.
The solutions of Chen (1985Chen ( , 1986Chen ( , 1991 and Yates (1988) are special cases of this study when r s → ∞, t inj → ∞, ω = 0, and V w,inj = 0. When r s → ∞, t inj → ∞, and V w,inj = 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.

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 inj = Q cha = 100 mL s −1 , t inj = 300 s, α 1 = 2.5 cm, α 2 = 2.5 cm, θ m = 0.30, θ im = 0.01, ω = 0.001 d −1 , R m 1 = R im 1 = R m 2 = R im 2 = 1, B = 50 cm, µ m 1 = µ m 2 = µ im 1 = im 2 = 10 −7 s −1 , and h w,inj = h w,cha = B. These parameters are from the experimental applications of Chao (1999), Chen et al. (2017), and Wang et al. (2018Wang et al. ( , 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 Table 4. Expressions of coefficients in solutions of Eqs. (20a)-(23c). 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.

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 −1 . The injection time is 5 h (t inj = 300 min) for the tracer while maintaining the same injection rate of 9.9 mL min −1 . 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 scaledependent 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 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.
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: where C OBS and C COM represent the observed and computed concentrations, respectively, Er is the error, and N is the 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. number of sampling points. In this study, the genetic algorithm (GA) is employed to search the optimal parameter values, such as θ m 2 , α 1 , and ω 1 for CDM of Eqs. (14) and (15) and θ m 2 , α 0 , k, and ω 1 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 scaledependent 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. 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 2 ) to compute the errors of best fitness with two BTCs simultaneously in Figs. 2 and 3, and R 2 is defined by  (14) and (15) are without the scale effect and Eqs. (17) and (18) are with the scale effect, respectively.
where C 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 lo- The overall R 2 shows the same observation as the overall Er, where the overall R 2 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.

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): 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. I j represents any one parameter of interest, like the volume of water in the wellbore (V w ), k, θ m = θ m 1 = θ m 2 , ω = ω 1 = ω 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): where I j is a small increment of I j .
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 0 ), 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 0 = 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.

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.) re- 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. quired 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 trans- fer 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 A i (·), B i (·) Airy functions of the first kind and second kind, respectively A i (·), B i (·) Derivative of the Airy functions of the first kind and second kind, respectively α 0 Longitudinal dispersivity (L) in the formation zone at r > r 0 α 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 m 1 , C im 1 Resident mobile and immobile concentrations (M L −3 ) of the skin zone, respectively C m 2 , C im 2 Resident mobile and immobile concentrations (M L −3 ) of the formation zone, respectively C um , C uim Resident mobile and immobile concentrations (M L −3 ) of the upper aquitard, respectively C lm , C lim Resident mobile and immobile concentrations (M L −3 ) of the lower aquitard, respectively C inj (t), C cha (t) Concentrations (M L −3 ) of the tracer in the wellbore in the injection and chasing phases, respectively C 0 Concentration (M L −3 ) of the tracer injected into the wellbore C w Concentration (M L −3 ) of the tracer in the wellbore D u , D l Vertical dispersion coefficients (L 2 T −1 ) of the upper and lower aquitards, respectively D 0 Molecular diffusion coefficient (L 2 T −1 ) h Hydraulic head (L) h 0 Hydraulic head (L) at r e h w,inj , h w,cha Water level in the wellbore in the injection and chasing phases (L) k A constant (dimensionless) ranging from 0 to 1 K 1 , K 2 Hydraulic conductivities (L T −1 ) of the skin and formation zones, respectively K d Equilibrium distribution coefficient (M −1 L 3 ) for the linear sorption process I m (·), K m (·) The mth-order modified Bessel function of the first and second kinds, respectively Q Pumping rate (L 3 T −1 ) (negative for injection and positive for pumping) Q inj , Q cha Well flow rates (L 3 T −1 ) in the injection and chasing phases, respectively. r Radial distance (L) from the center of the well r s Distance (L) from the center of the well to the outer boundary of the skin zone r w Radius (L) of the well r e Radial distance (L) from the center of the well to the outer boundary r 0 Radial distance (L) for the linear