Combining analytical solutions of Boussinesq equation with the modified Kozeny- Carman equation for estimation of catchment-scale hydrogeological parameters

Man Gao 1, 2, Xi Chen 1, 2 *, and Jintao Liu 3, 4 1 Institute of Surface-Earth System Science, Tianjin University, Tianjin, China 2 Tianjin Key Laboratory of Earth Critical Zone Science and Sustainable Development in 5 Bohai Rim, Tianjin University, Tianjin, China 3 State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing, China 4 Department of College of Hydrology and Water Resources, Hohai University, Nanjing, China 10

Abstract: Saturated hydraulic conductivity (K), drainable porosity (f), and effective aquifer thickness (D) are essential hydrogeological parameters for hydrologic modelling and 15 predicting. Streamflow recession analysis using analytical solutions of Boussinesq equation can yield estimated values for two of these three hydrogeological parameters when one is known a priori. In this study, we improved the inverse method for parameters estimation by combining the modified Kozeny-Carman equation with analytical solutions of Boussinesq equation to express the three hydrogeological parameters (K, f, and D) in 20 relation to catchment characteristics and recession constants in a sloping aquifer. Here, the three parameters can be estimated simultaneously from streamflow recession analysis.
Results of the estimated parameters are compared with the field measurements and the soiltexture based estimations in four small experimental catchments. It shows that our estimated values of these catchment-scale parameters can represent equivalent values in 25 the measured aquifer profiles/sites. In hilly areas, the slope aquifer takes a vital effect on the estimates of K and f. Neglecting the sloping effect can lead to overestimation of K and underestimation of f in 1 ~ 2 orders of magnitude in the study catchments. However, even in the hilly catchments, the estimated aquifer thickness D is much greater than that from measurements on hillslopes while it approaches riparian thickness, indicating that the 30 riparian zone takes a vital role on flow recession and the parameter estimations.

Introduction
Saturated hydraulic conductivity (K), drainable porosity (f), and aquifer thickness (D) are basic physical attributes of catchment and essential parameters for physically based 35 hydrologic models as well as for the hydrologic analysis of ungauged basins (Ali et al., 2014;Bogaart et al., 2016;Carrillo et al., 2011;Li et al., 2014;Vannier et al., 2016).
Extensive measurements of these parameters to represent their intrinsic variations in a catchment are costly and time-consuming. Thus, indirect methods, based on other more easily measurable soil properties (such as soil porosity or texture) and/or regular 40 hydrologic observation data (such as streamflow) are sometimes preferred.
The Boussinesq equation is a general differential equation for groundwater flow in unconfined aquifers of hillslope areas. The analytical and numerical solutions of the Boussinesq equation illustrates streamflow behaviors in relation to catchment properties like hillslope gradients, shapes, and hydrogeologic parameters of the aquifer (Troch et al., 45 2013). Such analytical solutions have been widely used for estimating catchment-scale hydrogeological parameters since the parameter estimation only needs the regular streamflow observations (Brutseart and Nieber, 1977;Brutseart and Lopez, 1998;Zhang et al., 2009).
The analytical solutions of Brutseart and Nieber (1977) for hydrogeological parameter 50 estimations were derived based on assumptions of specific initial and boundary conditions of an ideal aquifer (Szilagyi and Parlange, 1998) Mendoza et al. (2003) to derive the non-linear solutions in a drier and steeper landscape, where they used it for estimating hydraulic parameters at semi-arid mountainous watersheds in Mexico (Mendoza et al., 2003). Rupp 60 and Selker (2005) derived an analytical solution for the vertically heterogeneous aquifer on assumption of a power-law decreasing of K from aquifer surface to bottom. Gao et al. (2013) found that considering the high variability of vertical K could increase the reliability of the estimated values of f and vertical mean K at catchment scale. As an important component of hydraulic gradient, the slope of the aquifer is also crucial to accurately 65 express the groundwater flow behavior. Extended analytical solutions from horizontal to sloping aquifer were derived by Brutsaert (1994) and Hogarth et al. (2014) in different ways. Pauritsch et al. (2015) claimed that analytical models for sloping aquifers can narrow the ranges of the estimated values of K and D compared to that for horizontal aquifers in complex and heterogeneous aquifer systems in an alpine catchment. 70 Streamflow recession analysis for determining hydrogeological parameters is popularly based on two-stages flow recessions (Brutsaert, 2005;Troch et al., 2013) in two distinct flow regimes: the early-time and the late-time regimes. Correspondingly, two equations are derived to describe the flow recessions. Their solutions could be only used to estimate two of the three unknown parameters (K, D, and f). It means that one parameter must be given 75 a priori. For example, Zhang et al. (2009) estimated catchment-scale K and f with given D. Among the three parameters, K and f are dependent on soil properties. The drainable 80 porosity (f), defined as volume of water that a saturated soil will yield by gravity to the total volume of the rock or soil (Bear, 1972), is directly dependent on soil pore connected, soil particle-size distribution, and the soil pore structures (Ding et al., 2016). Saturated hydraulic conductivity (K) can be expressed as a function of the hydraulic radius of soil pore (defined as soil pore volume over the pore-solid surface area), indicating that K 85 depends on soil particle size, soil porosity, and pore shape factor and tortuosity (Zhang and Schaap, 2019). Evidently, both K and f can be estimated from easy-to-measure soil properties, like soil texture, bulk density, and organic content (Saxton and Rawls 2006;Stephens et al., 1998;Zhang et al., 2016), popularly using the pedotransfer functions (PTFs) (Saxton and Rawls 2006). 90 As both K and f depend on soil physical properties, K and f are significantly correlated.
Such hypothesis has been validated from field experiments (Nyberg, 1995 andEspeby, 1990) and physically based equations, including the Kozeny-Carman equation (Carman, 1937) and its modified equations (Ahuja et al., 1984(Ahuja et al., , 1989Rawls et al. 1993), which relate K to effective porosity (equal to saturated water content minus field capacity) that is viewed 95 identical to f for shallow homogenous unconfined aquifers (Bear, 1972). Therefore, if the analytical solutions of Boussinesq equation combine the modified Kozeny-Carman equation, the three unknown parameters (K, D, and f) can be determined simultaneously by using streamflow recession analysis.
The aim of this study is to extend the estimations of catchment-scale hydrogeological 100 parameters (K, f, and D) from two to three parameters by using assessable information of

Analytical solutions of Boussinesq equation for sloping aquifer 110
In hilly areas, a catchment can be viewed as an assembly of a series of sloping hillslopes along with river networks (Fig. 1a). To analytically estimate catchment-scale hydrogeological parameters, the non-uniformly distributed sloping aquifers can be simplified as an equivalently homogeneous aquifer (Fig. 1b).
where " is the flow rate per unit width of the aquifer, is water table, K is the saturated hydraulic conductivity, is the hillslope gradient, and x is the distance from the river to hillslope ridge. 120 Assuming that aquifer is isotropic and homogeneous, substituting Eq. (1) into the continuity equation yields: (2) where t is time and f is the drainable porosity.  Brutsaert (1994) derived the analytical solution of flow rate in the hillslope with specific 125 conditions, including: initially fully saturation of the aquifer, a sudden drawdown of the water table (hydraulic head η) to the channel as outflow starts for the flow condition at x=0 (Fig. 1b), and the infinitely wide aquifer that ensures negligibility of draining flow influence on the upper boundary of the aquifer. With that, the analytical solution of Eq. (2) for unit width hillslope is 130 where "8 is fast flow in the early-time recession, p is a constant equal to 0.3465, and D is the aquifer thickness. (4) where = tan /( ); @ = /2 for small in gentle hillslopes and thick aquifers, while @ = for large in steep hillslopes and thin aquifers.

Recession analysis
In the flow recession period, the relationship between streamflow ( ) and change of 140 streamflow (− / ) was proposed by Brutsaert and Nieber (1997) The function ( ) is often expressed by power-law equation as: where a and b are constants for a specific catchment. As streamflow in catchment outlet 145 can be viewed as the collection of flow from hillslopes along the river, Q is the multiple of for the early-time recession are obtained as: Inserting Eq. (4) into Eq. (6) Additionally, if values of recession segments in contiguous streamflow data are equal, only 170 the latest data are involved in the calculations using Eqs. (12) and (13).

Modified Kozeny-Carman equation relating to f
Kozeny (1927) derived a power-law function that relates the saturated hydraulic conductivity K to soil porosity based on the Hagen-Poiseuille's equation. It was later modified as the Kozeny-Carman equation after Carman (1937Carman ( , 1956, which accounts for 175 tortuosity in tube flow. Ahuja et al. (1984) proposed a modified Kozeny-Carman equation relating to effective porosity, which is identical to drainable porosity f in this study. The where is equal to 5.36 × 10 CW for K in a unit of m/s, obtained by fitting the measured data for 26 soil texture/porosity classes. The value λ ranges 0.165 ~ 0.694, estimated by fitting a log-log curve between water content and pressure head using the -33 and -1500 185 kPa water contents (Rawls et al., 1993), according to USDA soil texture classes and over 900 measurements. Eq. (14) has been fitted by a variety of soils, making the equation universally applicable.

Estimation of catchment-scale hydrogeological parameters
Substituting Eq. (14) into Eq. (7) gives Eq. (16), K and D can be obtained from implicit equations as follows According to Eq. (14), drainable porosity (f) can be expressed as Hence, f can be obtained easily when K is estimated from Eq.  (19)).

205
The proposed approach for catchment-scale hydrogeological parameter estimation was tested in four experimental catchments in the northern hemisphere. These include, the relict Schöneben Rock Glacier (SPG) in Austria, the Hemuqiao experimental catchment (HMQ) in China, the Panola Mountain Research Watershed (PMRW) and the WS10 from HJ Andrews Experimental Forest in the USA (Fig. 2) book. As to SPG, data published in a study by Pauritsch et al. (2015) are directly used in our analysis.
All the four catchments have field measures of K and D, which are used for validating the estimates of our proposed methods. In PMRW, soil K values were measured at two soil pits using constant head permeameters method by McIntosh et al. (1999), while K values 225 for saprolite and bedrock were measured by falling head permeameter (White et al., 2002).
In HMQ, the K values were measured at three typical soil pits by using falling head permeameter (Han et al., 2016). In WS10, vertical K values were measured at ten soil pits by using constant head permeameter (Harr, 1977). For SPG, the catchment equivalent values of K, f, and D were obtained from Pauritsch et al. (2015) according to field 230 measurements (Winkler et al., 2016).
Although the aquifer deposits in each of these catchments are dominated by one type of soil textures (Table 1), the point scale hydrogeological parameters observations show great heterogeneity, especially for saturated hydraulic conductivity K and soil/saprolite thickness D (Table 2). K shows highly vertical variations. The value of K at the upper soil layer can 235 be one order larger than the bottom one or the underlain saprolite (Table 2). Across catchments, the soil thickness D differs but the riparian areas much thicker than that the hillslopes for all catchments. For example, the soil thickness D for the typical soil profiles in PMRW ranges 0.6 ~ 1.6 m on hillslope but reaches 5 m at the riparian area (Peters et al., 2003). The average soil thickness on a typical hillslope in HMQ is 0.8 m (Han et al., 2018) 240 while the soil in riparian with colluvium and residuum deposits can reach up to 6 m (Han, et al., 2016) spatial variability of the saprolite thickness (0 ~ 5 m) over granodiorite (Peters et al., 2003) is much greater than that of the upper soil thickness (0.6 ~ 1.6 m).
Additionally, values of K and f for each catchment are estimated in terms of soil texture by using pedotransfer function proposed by Saxton and Rawls (2006) (Table 2).

Estimated recession constants
Plots of recession data in the form of − ef e3~ for each of the four catchments are shown in Fig. 3. In PMRW, HMQ, and WS10, the lower envelopes with 10% of data points are excluded in order to remove the effect of outliers on the envelope lines (Brutsaert and Lopez, 1977). For SPG, the recession data are adopted directly from Pauritsch et al. (2015). 255 Here, only winter data are used to fit the lower envelop lines as suggested by Pauritsch et al. (2015). The lower envelop lines with slopes of 3 and 1are used to derive the recession intercepts for early-and late-time recessions ( 8 and D ), respectively. Compared to this range, the late-time recessions are relatively fast in our catchments except for SPG where the soil deposits are much thicker.

Estimated catchment-scale hydrogeological parameters
Given catchment physical properties (area, river length, and slope in Table 1), the estimates 270 of af and as as well as the soil texture-based λ (Table 3) The estimated catchment-scale K is 75% larger than the geometric mean of the measured K in the vertical profile for these three catchments. An exception is that the estimated K value (1.75×10 -6 m/s) in WS10 is much below the lower limit of the measured 295 values in a shallow soil profile (4.4×10 -5 ~ 1.1×10 -3 m/s). The two orders of magnitude smaller K value estimated in WS10 could be attributed to the fact that the measurements were only taken from the upper soils (1.5 m in maximum) with abundant macropores (Harr, 1977), which could be much greater than that in the soil below (e.g. ~5.77 m thickness in our estimation). Fig. 4b shows the estimated catchment-scale K in this study vs. the 300 estimated K from soil texture. Both estimated K values are close to the 1:1 line, indicating that both estimates are in agreement. The f values from our study method are close to the estimated f from soil texture with slight overestimation by 4% ~ 47% (Fig. 4c,  The estimated aquifer thickness D values in our study are comparable to the measured 305 values in SPG (Fig. 4d). But PMRW and HMQ showed significantly larger than soil thickness on the hillslope, the magnitudes are within the range of mean soil thickness on hillslopes and the maximum soil thickness in riparian. In WS10, the estimated D is close to the measured D when both soil and saprolite are considered as effective groundwater aquifers. 310

Sensitivity analysis for the hydrogeological parameters
As 5 TU n 6 A in Eq. (9) is proven to be much larger than 4 (Table 3) α, λ, and γ. Among all the independent variables, the catchment slope α is the most sensitive variable to the three hydrogeological parameters (Fig. 5). The next sensitive variable is late-time recession rate ( D ) for K and D, and γ for f. The most insensitive variables are λ and 8 for f and D, and, λ and γ for K. Among the three parameters, K is the most sensitive 330 parameter to the independent variables (α, D , 8 , and L).
Overall, accurate representation of catchment aquifer slope α and analysis of late-time recession could increase the reliability of the three estimated parameters, especially for K.

Comparison of parameter estimates from analytical solutions in horizontal and
Combing Eqs. (7) with (9), K and f for sloping aquifer can be expressed as  Gao et al (2018). Whilst the measured K values are in the order of magnitude of 10 -5 and 10 -6 m/s in these four catchments. Thus, K could be extremely overestimated by horizontal aquifer assumption.

Effect of vertical heterogeneity of soils on the determination of K and f
For natural soils, K usually decreases with soil depth (Brooks and Boll, 2004) and field 365 measurements also indicate that f could decrease with soil depth as well (Harr, 1977;McGuire et al., 2007). The decrease of K with depth could reduce the late-time recession rate as the water level becomes lower. The depth dependent K affects estimates of hydrogeological parameters from streamflow recessions. According to our previous study (Gao et al., 2013) (e.g. PMRW, HMQ, and WS10) shows that K at the upper soil can be 1 ~ 2 orders of magnitude greater than that in the deep deposits. The estimated K from our proposed method can be regarded as an equivalent value that represents K in a specific range of the 375 soil profile. For example, in PMRW, the estimated value of K is more likely to represent the equivalent value between soil and saprolite. In HMQ, the estimated value of K is closer to the measures within the soil profile at 30 ~ 90 cm depths even though the soil profiles on the hillslope can be as deep as 1.2 ~ 4.4 m (Fig. 4a). The ranges of the measured K cover our estimated values in three catchments, indicating that the measured ranges of K can be 380 used to restrict parameter variations in model calibration. However, it should be mentioned that the selected profiles and sites for measures should be representative of a catchment.
This study shows that the field measurements only taken at the upper soils (1.5 m in the maximum depth) cannot capture the whole profile values of 5.77 m in depth in WS10. Burns et al. (2001) demonstrated that the hillslopes contribute flow mainly during large storm event while the riparian zone in the gentle areas with thick deposits of sediment takes a vital role on runoff generation and flow regulation (Clark et al., 2009;McGlynn and McDonnell, 2003;Klaus et al., 2015). A riparian zone can generate a vast majority of runoff during small event and between events (recession processes) (Burns et al., 2001;McGlynn 390 and McDonnell, 2003), and contribute a primary flow component throughout the recessions (Clark et al., 2009). Our estimated catchment-scale aquifer thickness ranges between the measured thicknesses in hillslope and riparian zone but is much thicker than the thickness in hillslope areas. It implies that thicker deposits in the riparian zone have involved to slow down flow recessions since the thicker aquifer reduces hydraulic conductivity K as shown 395 in Fig. 6. Thus, detail investigations on the riparian zone properties, such as the deposit thickness and hydrologic connectivity with hillslopes and streams are vital for appropriate analysis of flow recession behaviors.

Conclusions
Catchment-scale hydrogeological parameters (K, f, and D)  thickness is much greater than the measured thickness on hillslopes but close to the thickness in the riparian zone, indicating that riparian zone could take an important role on 420 flow recessions even in the hilly catchments. ** The SPG surface sediment characteristics is relative coarse and is even descripted as "covered by coarse-grained, blocky material consisting of gneissic rocks ranging from cubic decimeters to a few cubic meters" (Pauritsch et al., 2015). Its soil texture is set as loamy sand combining this description with comparison of the hydraulic properties (K and f) measured with related soil.