Forchheimer flow to a well-considering time-dependent critical radius

Previous studies on the non-Darcian flow into a pumping well assumed that critical radius ( RCD) was a constant or infinity, whereRCD represents the location of the interface between the non-Darcian flow region and Darcian flow region. In this study, a two-region model considering time-dependent RCD was established, where the non-Darcian flow was described by the Forchheimer equation. A new iteration method was proposed to estimate RCD based on the finite-difference method. The results showed that RCD increased with time until reaching the quasi steady-state flow, and the asymptotic value of RCD only depended on the critical specific discharge beyond which flow became nonDarcian. A larger inertial force would reduce the change rate of RCD with time, and resulted in a smaller RCD at a specific time during the transient flow. The difference between the new solution and previous solutions were obvious in the early pumping stage. The new solution agreed very well with the solution of the previous two-region model with a constant RCD under quasi steady flow. It agreed with the solution of the fully Darcian flow model in the Darcian flow region.


Introduction
Darcy's law indicates a linear relationship between the fluid velocity and the hydraulic gradient (Bear, 1972), which is a basic assumption used to handle a great deal of problems related to flow in porous and fractured media.However, many evidences from the laboratory and field experiments show that this linear law may be invalid in some situations, especially when the groundwater flow velocity is sufficiently high or sufficiently low, where non-Darcian flow prevails (Basak, 1977;Bordier and Zimmer, 2000;Engelund, 1953;Forchheimer, 1901;Izbash, 1931;Liu et al., 2012;Soni et al., 1978) Darcy's law considers kinematic forces but excludes inertial forces of flow.However, the inertia forces become significant in respect to the kinematic forces when the velocity is great, leading to non-Darcian flow (Engelund, 1953;Forchheimer, 1901;Irmay, 1959;Izbash, 1931).Forchheimer (1901) proposed a heuristic Forchheimer law describing the non-Darcian flow, which is an extension of Darcy's law by adding a second-order velocity term, representing the inertial effect.To verify the applicability of the Forchheimer law, many approaches were introduced, such as the dimensional analysis (Ward, 1964), the capillary model (Dullien and Azzam, 1973), the hybrid mixture theory (Hassanizadeh and Gray, 1987), and the volume averaging method (Whitaker, 1996).Recently, Giorgi (1997) and Chen et al. (2001) analytically derived the Forchheimer law from the Navier-Stokes equation.Another widely used model describing the non-Darcian flow is the Izbash equation (Izbash, 1931).This equation is a fully empirical power-law function obtained through analyzing experimental data.The Izbash equation is preferred for modeling purpose, since the power index in the Izbash equation can be parameterized depending on flow conditions (Basak, 1977).Some researches demonstrated that the Forchheimer and Izbash equations were identical for some cases (George and Hansen, 1992).
Due to the high velocities, non-Darcian flow is likely to occur near pumping/injecting wells (Yeh and Chang, 2013;Wen et al., 2008b).Several studies showed that the non-Darcian effect had significant influence on hydraulic parameter estimations.For instance, Theis solution cannot be used to explain the pumping test data in the Chaj-Doab area near Gujrat water distributory in Pakistan (Ahmad, 1998), while Birpinar and Sen (2004) and Wen et al. (2011) found that the Forchheimer law worked very well.Quinn et al. (2013) demonstrated that non-Darcian flow effect increased as the initial applied head differential increased in a series of slug tests.Specifically, Quinn et al. (2013) showed that the hydraulic conductivity was underestimated by Darcy's law when the initial applied head differentials were greater than 0.2 m.They pointed out that Darcian flow conditions can be maintained in fractured dolostone and sandstone when the initial applied head differentials were less than 0.2 m (Quinn et al., 2013).Mathias Introduction

Conclusions References
Tables Figures

Back Close
Full and Todman (2010) showed that the Jacob method, based on Darcy's law, cannot fit the step-drawdown tests of van Tonder et al. (2001) when the pumping rate is greater 10 m 3 h −1 .However, the Forchheimer law fitted the step-drawdown tests data very well (Mathias and Todman, 2010).In this study, we will focus on the non-Darcian flow into a pumping well by the Forchheimer law.
Although many efforts have been devoted to study the non-Darcian flow around the well, the exact solutions have not been obtained due to the non-linearity of the problem (Mathias et al., 2008;Yeh and Chang, 2013).For example, Sen (1990Sen ( , 2000) ) employed the Boltzmann transform to analytically solve the problems related to the non-Darcian flow.This method was showed to be problematic, since both initial and boundary conditions cannot be simultaneously transformed into a form only containing the Boltzmann variable (Camacho and Vasquez, 1992;Wen et al., 2008a).Wen et al. (2008a, b) obtained the semi-analytical solutions of the non-Darcian flow model by combining the linearization procedure and the Laplace transform method (LL method).This LL method assumed that transient flow in the non-Darcian flow region can be taken as a quasi-steady state flow.Wen et al. (2008aWen et al. ( , 2008b) ) pointed out that solutions by the Boltzmann transform and the LL method coincided at late time.To test the accuracy of the semi-analytical solutions (Wen et al., 2008a;Sen, 2000), Mathias et al. (2008) and Wen et al. (2009) employed the finite-difference method to study the non-Darcian flow problems, and their results showed that the semi-analytical solution only agreed very well with the numerical solution at late pumping stage.
All above-mentioned investigations assume that the non-Darcian flow occurs over the entire domain, which is called a fully non-Darcian flow (F-ND) model hereinafter.In fact, the regime of the flow to the pumping well can be divided into two regions: non-Darcian flow occurs within a narrow region around well, due to the relatively high velocity of flow there, and Darcian flow prevails over the rest domain.One may think that such two-region flow could be described by the Forchheimer law, which would automatically reduce to the Darcy's law at the location far from the well (because the second-order velocity term in the Forchheimer law will be negligible if velocity approaches zero).Introduction

Conclusions References
Tables Figures

Back Close
Full However, Mackie (1983) demonstrated that the two-region model could fit the experimental data in the laboratory better than the F-ND model.Huyakorn and Dudgeon (1976) employed a two-region model to study flow near a pumping well.Basak (1978) presented analytical solutions of the two-region model for steady-state flow to a fully penetrating well.Sen (1988) and Wen et al. (2008b) derived the analytical solutions of the two-region model for transient flow to a pumping well.All researches mentioned above implied that the critical radius is a constant, where the critical radius represents the location separating the non-Darcian and Darcian flows (Sen, 1988;Wen et al., 2008b).For example, the critical radius is infinity for the F-ND model and is zero for the fully Darcian flow model, while it is a finite constant for the two-region model, and its value is computed under the quasi-steady state flow condition (Sen, 1988;Wen et al., 2008b).Actually, the critical radius changes continuously with time for the transient flow, and cannot be determined straightforwardly.This is particularly true for the variable-rate pumping tests (Bear, 1972;Mishra et al., 2012), slug tests (Quinn et al., 2013) or the step-drawdown tests (Louwyck et al., 2010;Mathias and Todman, 2010), since the regime of the flow around the well drastically changes with time.
In this study, we will investigate non-Darcian flow into a fully penetrating pumping well considering a time-dependent critical radius using the finite-difference method.A new iteration procedure will be proposed to estimate the moving critical radius.This new model reduces to the F-ND model when the critical radius is infinite and it becomes the fully Darcian flow model when the critical radius is 0. post-linear laminar flow, and (D) non-Darcy post-linear turbulent flow (Basak, 1977;Bear, 1972).For radial flow to a pumping well, the velocity in the aquifer decreases with the distance from the well.Therefore, the radial flow might experience all four-flow regimes.To simply the problem, we use a two-region model that considers a non-Darcian flow region near the well and a Darcian flow region away from the well.
A unique feature of the two-region model used in this study is that the critical radius is allowed to vary with time whereas it was assumed to be constant in previous studies (Dudgeon et al., 1972a, b;Huyakorn andDudgeon, 1976, 1978;Mackie, 1983;Sen, 1988;Wen et al., 2008b).
Generally, the start of the non-Darcian flow can be determined by the critical Reynolds number (Re C ), where the Reynolds number is defined as where ν is the kinematic viscosity of the fluid (L 2 T −1 ); D p is the characteristic grain diameter (L); q(r, t) is specific discharge (LT −1 ) at distance r (L) and time t (T); Re is Reynolds number which depends on time and space (dimensionless).The critical Reynolds number (Re C ) refers to Re at the start of non-Darcian flow.Up to present, there is still considerable debate on Re C for the initiation of non-Darcian flow in porous media.Bear (1972) pointed out that Re C varied between 3 to 10; Scheidegger (1974) gave Re C to be 0.1 to 75; Zeng and Grigg (2006) suggested the range of Re C as 1 to 100.Re C will be set as 10 in this study.According to Eq. ( 1), one can see that the specific discharge has a linear relationship with Re.Therefore, the critical specific discharge (q C ) can also be used to determine the start of the non-Darcian flow, since one can calculateq C for a given Re C .When the specific discharge is less than or equal to q C (or Re ≤ Re C ), the flow is considered as Darcian.When the specific discharge is greater than q C (or Re > Re C ), the flow is taken as non-Darcian.Denoting R C (t) as the critical radius at which q = q C (or Re = Re C ), then it is non-Darcian flow when r ≤ R C (t) and Darcian flow when r > R C (t), as shown in Fig. 1.Introduction

Conclusions References
Tables Figures

Back Close
Full For the quasi-steady state flow, one has where B is the thickness of the aquifer (L); Q is the well discharge (L 3 T −1 ).For the case with a constant pumping rate, R C is also a constant for a specific Re C .This constant R C was used in previous two-region models of transient non-Darcian flow (Sen, 1988;Wen et al., 2008b).Actually, R C is not a constant for transient flow, and it cannot be determined directly since the velocity distribution changes with time.In this study, a new iteration method will be proposed to determine R C as described below.

Mathematic model
Figure 1 shows the physical model investigated in this study, where a pumping well fully penetrates a confined aquifer.The origin of the cylindrical coordinate system is at the center of the well.The r axis is horizontal and outward from the well, and the z axis is upward vertical.Two assumptions are made in this study.First, the non-Darcian and Darcian flow may coexist and the critical radius is time-dependent, and the non-Darcian flow is governed by the Forchheimer law.Second, the system is hydrostatic before the pumping starts, so R C (t = 0) = 0.These assumptions, although quite idealized, are standard in well hydraulic study (Papadopulos and Cooper, 1967;Sen, 1988;Wen et al., 2008b).Based on these assumptions, the governing equations of the two-region flow model can be described as follows where s Y (r, t) and s N (r, t) are drawdowns (L) at distance r and time t in Darcian flow and non-Darcian flow regions, respectively; S is the aquifer storage coefficient (dimensionless).
The outer boundary condition is Assume that the pumping rate is large enough to induce non-Darcian flow near the well.Here we will consider the wellbore storage with a finite diameter well and the boundary condition at the wellbore can be written as where Q is positive for the pumping rate; r w is the radius of the well (L); s w is the drawdown inside the well (L).Notice that well loss is not considered so the drawdown is continuous across the well screen The drawdown and the discharge from the Darcian flow region into the non-Darcian flow region are continuous at the critical radius In the non-Darcian flow region, we use the Forchheimer law to describe the flow (Forchheimer, 1901) in which β (TL −1 ) and K β (LT −1 )are empirical constants depending on the properties of the medium (Sidiropoulou et al., 2007).β is called the inertial force coefficient.K β is called the apparent hydraulic conductivity and it reduces to the hydraulic conductivity when β = 0 (Chen et al., 2001;Sidiropoulou et al., 2007).
In the Darcian flow region, one has Equations (3-12) compose of the mathematical model of the two-region model with a time-dependent critical radius R C (t).This new model is an extension of the previous model by Sen (1988).When R C (t) → ∞, this model becomes the F-ND model.When R C (t) = 0, it reduces to the fully Darcian flow model.

Dimensionless transformation
Defining the dimensionless variables in Table 1, Eqs. (3-12) can be rewritten as Notice that a negative sign has been used for defining q D in Table 1.(Eq.7) in the dimensionless form is The dimensionless Forchheimer law becomes where β D is the dimensionless inertial force coefficient.When r D > R CD , groundwater flow follows the Darcy's law in the dimensionless format as where λ is the ratio of the hydraulic conductivity and apparent hydraulic conductivity, and it is usually taken as unity (Sidiropoulou et al., 2007).

Numerical solution
Because of the non-linearity of the problem, it is not easy to obtain the analytical solution of drawdown even if R CD (t D ) is constant.In this study, we will employ the finitedifference method to investigate the problem considering a time-dependent R CD (t D ).
Due to the axisymmetric nature of the problem, the numerical simulation will be conducted with a non-uniform grid system, where the radial steps are smaller near the well and become progressively greater away from the well.Similar to previous studies (Mathias et al., 2008;Wen et al., 2009), we discretize the dimensionless space r D logarithmically.The dimensionless space domain [r wD , r eD ] is discretized into N nodes excluding the two boundary nodes r wD and r eD , where r eD is a relatively large dimensionless distance used to approximate the infinite boundary (Mathias et al., 2008;Wen Introduction

Conclusions References
Tables Figures

Back Close
Full  , 2009).For any node of r i , r wD < r i < r eD , i = 1, 2, . .., N, one has where r i +1/2 is calculated as follows log 10 (r i +1/2 ) = log 10 (r wD ) + i log 10 (r eD ) − log 10 (r wD ) After the spatial discretization, Eqs. ( 13) and ( 14) become where q YD,i and s YD,i are the dimensionless specific discharge q YD and dimensionless drawdown s YD at node i for the Darcian flow, respectively; q ND,i and s ND,i are the dimensionless specific discharge q ND and dimensionless drawdown s ND at node i for the non-Darcian flow, respectively.In terms of the Forchheimer equation of Eq. ( 20), one can obtain and

Conclusions References
Tables Figures

Back Close
Full where node N s means the location of R CD (t D ).At the well-aquifer boundary, one has where s wD is the dimensionless drawdown inside the well.Considering Eq. ( 19), s wD can be approximated as follows When r D > R CD , the finite-difference scheme of the specific discharge can be obtained from Eq. ( 21): As for the boundary at the infinity, the finite-difference scheme is Now one obtains a set of ordinary differential equations.It is notable that R CD or N s which is related to the index i in Eqs. ( 26) and ( 27) and Eqs. ( 30) and ( 31) is timedependent.In the following section, a new iteration method will be proposed to determine the values of R CD or N s .Introduction

Conclusions References
Tables Figures

Back Close
Full Before introducing the new iteration method, the relationship between R CD and the velocity distribution will be investigated first, based on the two-region model with a constant R CD .The values of the constant R CD are set as 0, 0.02, 0.04, 0.08 and 0.50, respectively.The other parameters are r wD = 1 × 10 −5 , β D = 20, λ = 1.The mathematic model with a constant R CD will be solved by the finite-difference method.
Figure 2a shows the specific discharge distributions with different R CD of 0, 0.02, 0.04, 0.08 and 0.50.The curve of R CD = 0 represents the fully Darcian flow model.One can find that the specific discharge decreases with increasing R CD at a given r D , starting from its maximum at R CD = 0 (Darcian flow).This observation is understandable.The increasing R CD implies a stronger contribution of the inertial effect, which also means a larger resistance to flow, thus it leads to a smaller specific discharge.After trying many different sets of aquifer parameters, numerical simulation indicates that this observation is universally valid.This observation will serve as the basis for the new iteration method to seek the location of R CD (t D ).
Similar to the use of Re C to determine the start of the non-Darcian flow, one can use q CD for the initiation of the non-Darcian flow, where q CD is the dimensionless critical specific discharge defined in Table 1.We denote r j CD as the newly computed critical radius at the j th step of the new iteration method, where j = 1, 2, 3, . . .Since the aquifer system is initially hydrostatic, the initial critical radius r 0CD is set as 0. For a given dimensionless time t 1D , the detailed procedures of the iteration method for searching R CD (t 1D ) will be introduced as follows.Firstly, the specific discharge distribution in the aquifer can be calculated using Eqs.(24-32) with R CD (t 1D ) = r 0CD , as shown in Fig. 2b.Based on the computed specific discharge distribution, one can find the new critical radius r 1CD according to a given constant q CD .Secondly, the new specific discharge distribution can be similarly calculated using Eqs.(24-32) with R CD (t 1D ) = r 1CD , and the new critical radius r 2CD can be obtained according to q CD .It is notable that r 1CD and r 2CD serve as the upper and lower limits for searching R CD (t 1D ), as illustrated in Introduction

Conclusions References
Tables Figures

Back Close
Full Fig. 2b.Similarly, one can estimate the new critical radius r 3CD using r 2CD , where r 3CD is located somewhere between r 1CD and r 2CD .Following the same procedures, a new critical radius r 4CD can be calculated based on r 3CD , and r 4CD is between r 2CD and r 3CD .One can repeat above computations until the new critical radius finally converges.For the actual problems, we define a convergence criterion R old CD − R new CD ≤ ξ, where R old CD and R new CD are the critical radius for the previous step and present step, respectively; ξ is a small positive value such as 0.001.If this criterion is satisfied, the new critical radius r jCD is thought as the estimation ofR CD (t 1D ).We develop a MATLAB program named as Two-Region Model with Moving critical radius (MTRM) to facilitate the computation.Figure 3 represents the flow chart of the MTRM algorithm, where t k is the time at time step k; k max is the total number of the time steps; dt D is the dimensionless time step; s D,i and q D,i are the dimensionless drawdown and dimensionless specific discharge at node i in the aquifer respectively.

Comparison with the previous solutions
To test the new solution, the fully Darcian flow solution of Papadopoulos and Cooper (1967), the fully non-Darcian flow solution of Mathias et al. (2008) and the two-region model of Sen (1988) will be introduced.Using the notations of this study, the original solutions of Sen (1988) in the Darcian flow region (r D ≥ R CD ) become where , and erf is the error function sign.
In the non-Darcian flow region (r D < R CD ), In the early stage, the differences among three previous solutions and the new solution of this study are obvious, as shown in Fig. 4a.Firstly, the solution of Papadopoulos and Cooper (1967) is smaller than the others near the well.This is because the inertial Introduction

Conclusions References
Tables Figures

Back Close
Full forces of the non-Darcian flow increase the resistance for flow, thus resulting in drawdown greater than those for the Darcian flow near the well.The second is that the F-ND solution agrees with the new solution near the well.The third is that the solution of Sen (1988) does not agree with the new solution at the early time.This is probably because of the Boltzmann transform method used by Sen (1988) to deal with the non-Darcian flow at the early time, which has been discussed in several previous studies (Camacho and Vasquez, 1992;Wen et al., 2008b).
In the late pumping stage, the transient flow approaches the quasi-steady state, and the specific discharge distribution is invariant with time according to Eqs. ( 3) and (4) or Eqs. ( 13) and ( 14), regardless of the Darcian flow or non-Darcian flow.Under the quasi-steady state flow condition, the critical radius obtained by this new solution becomes a constant which is the same as the one used by previous two-region models such as Sen (1988) and Wen et al. (2009).Therefore, it is not surprise to see that the new solution agrees with that of Sen (1988) very well at late time (see Fig. 4b).Another observation is that the new solution agrees with the solution of Papadopoulos and Cooper (1967) in the Darcian flow region, and with the solution of the F-ND model in the non-Darcian flow region near the well.

Effect of the inertial force coefficient to the critical radius
For the non-Darcian radial flow near the pumping well investigated here, the dimensionless inertial force coefficient (β D ) is of primary concern.To test the effect of β D to the groundwater flow, the representative parameters are chosen as follows: the value of time until the flow approaching the quasi-steady state condition.In the early pumping stage, the specific discharge is very large near the well and decreases quickly with the distance from the well, so R CD is very small.With time going, the cone of depression will expand along the radial direction and the slope of the cone of depression becomes flatter, so R CD becomes greater.Secondly, a larger β D would reduce the rate of change R CD vs. time, thus result in longer time to approach its asymptotic value, and consequently leads to a smaller R CD at a specific time in the transient state (see Fig. 5).This is because that a larger β D implies a stronger inertial force, which increases the resistance of flow.The third interesting observation is that the asymptotic value of R CD is the same for differentβ D .This can be explained using Eq. ( 2).Based on the definition of the dimensionless parameters defined in Table 1, Eq. ( 2) becomes Therefore, the value of R CD has no relationship with β D under the quasi-state state flow condition, while it only reciprocally depends on the critical specific discharge.

Effect of the critical specific discharge to the critical radius
The criterion to judge the initiation of the non-Darcian flow is an important factor of concern.Up to now, there is still considerable debate on what value of Re C to use for the start of non-Darcian low.The recommended values of Re C range from 0.1 to 100 for porous media flow (Bear, 1972;Scheidegger, 1974;Zeng and Grigg, 2006).
To check the influence of Re C on R CD during the transient flow, the values of q CD are chosen as 1, 2 and 5 considering the direct relationship of q CD and R CD in Eq. ( 2).The other parameters are β D = 1, and r wD = 1 × 10 −4 .Figure 6 represents the effect of q CD on R CD .It is obvious that the asymptotic value of R CD is equal to 1/q CD , as reflected in Eq. ( 37).Another interesting observation is that R CD decreases with increasing q CD , and it takes shorter time for R CD to approach its asymptotic value.The reason can be also explained using Eq. ( 37).Type curves are a series of curves that reveal the functional relationship between the well functions (or drawdown) and the dimensionless time factors (Sen, 1988;Wen et al., 2011).Type curve is one of the common approaches to identify the aquifer parameters or to predict the drawdown (Sen, 1988;Wen et al., 2011).Sen (1988) presented different type curves in the Darcian flow region and non-Darcian flow region based on a two-region model.In that model (Sen, 1988), R CD was a fixed value which only depends on the rate of pumping but independent of time.In this study, R CD changes with time, and the type curves might be different from the ones generated by Sen (1988).To investigate the behaviors of the type curves of the new solution, the value of q CD is set to be 2, and the two observation locations will be chosen, r D = 0.1 and 1.0.According to Eq. ( 37), the maximum of R CD is 0.5 at the quasi-steady state, so the flow at r D = 0.1 will experience both Darcian flow (at the early time) and non-Darcian flow (at late time), while the flow at r D = 1.0 is always Darcian.
Figure 7 shows the time-drawdown at r D = 0.1 for different dimensionless inertial force coefficients in the log-log scale.Two interesting observations can be seen from this figure.The first observation is that there is a deflection point in the curve, and the time of this deflection point becomes longer with increasing β D .This is because a larger β D implies a stronger inertial effect, which leads to a larger drawdown and longer time to approach the quasi-steady state condition.This observation is not found in the F-ND model (Wen et al., 2011) and in the two-region model (Sen, 1988).The second observation is that the drawdown in the quasi-steady state increases with increasing β D , and the reason for this has been explained in previous studies (Wen et al., 2011).Figure 8 represents the time-drawdown at r D = 1 in the semi-log scale.One notable point is that flow at r D = 1 is always Darcian, so there is no deflection point in the type curves.The differences among the curves with different β D are very small at the beginning, then they become larger with time going, and finally approach the same value at the quasi-steady state.This is because the specific discharge far from the well Introduction

Conclusions References
Tables Figures

Back Close
Full

Conclusions References
Tables Figures

Back Close
Full an empirical constant in the Forchheimer law (T L −1 ), named as inertial force coefficient in this study ν kinematic viscosity of the fluid (L 2 T −1 ) q ND , q YD dimensionless specific discharges defined in Table 1 in the non-Darcian flow and Darcian flow regions, respectively q CD dimensionless critical specific discharge defined in Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | of the critical radius of the two-region model Previous researches showed that the porous media flow may be divided into four regimes, such as (A) non-Darcy pre-linear laminar flow, (B) Darcy flow, (C) non-Darcy Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | The subscript D means the dimensionless variables.The boundary condition with the wellbore storage HESSD Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | et al.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 4 Iteration method to determine R CD or N s Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

)
Figures 4a and b show the distance-drawdown curves in the early and late pumping stages, respectively.In these two figures, Papadopoulos and Cooper (1967) represents the analytical solution of the fully Darcian flow model, Sen (1988) is the analytical solution of the two-region model by the Boltzmann transform method, and Mathias et al. (2008) represents the numerical solution of the fully non-Darcian flow model.The deflection point of the curve is the location of the critical radius.
Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | According to Eqs. (1)-(2) and the definition of the dimensionless critical specific discharge, one has q CD = 2 and R C = 2.4 m.The values of β D of this study are the same as the ones used by Wen et al. (2011), and they are 1, 10, and 50.Dimensionless radius of the well is r wD = 1 × 10 −4 .Figure 5 shows the critical radius (R CD ) changes with time for different dimensionless inertial force coefficients.Several observations can be seen.Firstly, R CD increases with HESSD Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 5.4 Type curves in the non-Darcian flow region and Darcian flow region Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |is very small at the beginning, so the drawdown approaches 0. With time going, the inertial force affects the groundwater flow at r D = 0.1.At the quasi-steady state, the drawdown in the Darcian flow region does not change with β D , where the reason has been explained in Sect.5.1, as shown in Fig.4b.6 Summary and conclusionsIn this study, a new two-region flow model considering the time-dependent critical radius (R CD ) is established to investigate the groundwater flow into a pumping well, and a new iteration method is proposed to estimate R CD , based on the finite-difference method.The convergence of this iteration method has been verified.In the non-Darcian flow region, the flow is governed by the Forchheimer equation, and the start of the non-Darcian flow is determined by the critical specific discharge, which is calculated by the critical Reynolds number.The new solution is compared with previous solutions, such as the fully Darcian flow model, the two-region model with a constant critical radius, and the fully non-Darcian flow model.The impacts of the dimensionless inertial force coefficient (β D ) and dimensionless critical specific discharge (q CD ) on the critical radius and flow field have been analyzed.Several findings can be drawn from this study:1.In the early stage, the new solution agrees with the fully non-Darcian flow solution near the well, differs with the fully Darcian flow model ofPapadopoulos and  Cooper (1967)  and the two-region model ofSen (1988).2. In the quasi-steady flow stage, the new solution agrees with the solution of Sen (1988) very well.It agrees very well with the solution of the fully Darcian flow model (Papadopulos and Cooper, 1967) in the Darcian flow region, and with the solution of the fully non-Darcian flow model (Mathias et al., 2008) in the non-Darcian flow region near the well.3. R CD increases with time until reaching the quasi-steady state flow, and the asymptotic value of R CD only depends on q CD .A larger β D would reduce the rate of Discussion Paper | Discussion Paper | Discussion Paper | change of R CD with time, and result in a smaller R CD at a specific time during the transient flow state.4.There is a deflection point in the type curve when the observation well location is within the non-Darcian flow region in the quasi-steady state, and the time associated with this deflection point becomes longer with a larger β D .Discussion Paper | Discussion Paper | Discussion Paper | the aquifer (L T −1 ) K β apparent hydraulic conductivity, an empirical constant in the Forchheimer law (L T −1 ) q specific discharge in the aquifer (L T −1 )q C critical specific discharge (L T −1 ) q Y , q N specific dischargesfor Darcian flow and non-Darcian flow (L T −1 ), respectively Q well discharge (L 3 T −1 ) s drawdown for aquifer (L) s Y , s N drawdowns for Darcian flow and non-Darcian flow (L), respectively s w drawdown inside well (L) S storage coefficient of the aquifer (dimensionless) r distance from the center of the well (L) r w radius of the well screen (L) R C critical radius for non-Darcian flow (L) Re Reynolds number (dimensionless) Re C critical Reynolds number (dimensionless) t pumping time (T) β

Fig. 1 .
Fig. 1.The schematic diagram of the non-Darcian flow into a fully penetrating pumping well considering the time-dependent critical radius.

Table 1 .
Dimensionless variables used in this study.

Table 1 r
D dimensionless distance defined inTable 1 r wD dimensionless radius of the well screen defined in Table 1 R CD dimensionless critical radius defined in Table 1 s ND , s YD dimensionless drawdown s defined in Table 1 in the non-Darcian flow and Darcian flow regions, respectively s wD dimensionless drawdown inside the well defined in Table 1 t D dimensionless time defined in Table 1 β D dimensionless inertial force coefficient defined in Table 1 λ ratio of the hydraulic conductivity and apparent hydraulic conductivity defined in Table 1