MODELLING OF HILLSLOPE STORAGE UNDER TEMPORALLY VARIED RAINFALL RECHARGE

Abstract. Water storage inside hillslopes is a crucial issue of environment and water resources. This study separately built a numerical model and an analytical model employing a hillslope-storage equation to simulate the water storage in a sloping aquifer response to recharge. The variable width of hillslope was hypothetically represented by an exponential function to categorize the hillslope into three types: uniform, convergent, and divergent. An integral transform technique was introduced to derive the analytical solution whereas a finite difference method was employed for the numerical modelling. As a result, under the same scenario a gap existed between the two solutions to distinct forms of the water storage equation, and the gap decreases with a falling recharge rate for convergent hillslopes. Moreover, all outflows gradually approached one value based on different hillslopes under the same accumulative recharge amount for six typical rainfall recharge patterns. Particularly, while the recharge stops, the outflow decreases and then mildly rises for a long time for convergent hillslope because of the slow water release near the upstream boundary where the storage water is relatively abundant due to the widest width.


Introduction
In Taiwan, surface water flow down a hillslope cause considerably severe soil erosion on the ground. Consequently, the management of soil and water conservation in watersheds has been a crucial issue in Taiwan. In general, water transportation, hillslope shape, soils erosion and transport, and aquifer configuration are the significant factors changing a watershed [3,9,10]. This study is focused on the hill storage of different hillslope forms including convergent, uniform, and divergent hillslopes.
In the past, some lab experiments and in-situ observations had explored the water storage under the ground. For instance, Anderson and Burt [1] measured soil moisture content through an automatic detection system and concluded that topography is the major influencing factor. Moseley [12] surveyed surface flow and subsurface flow in a forest catchment and found that subsurface flow is highly reduced on gentle slopes. Subsequently, O'Loughlin [13] proposed a method of topographic analysis to estimate the water-saturated zone of a catchment area. Later, McDonnell [11] revealed that the flow speed of groundwater is influenced by the permeability and slope of an aquifer using an isotope analysis in a watershed. Genereux et al. [8] traced surface water flow from distinct upstream areas to the outlet in a watershed using a chemical method, and found the travel time of water flow can be estimated from the topography of catchment. Likewise, Woods and Rowe [26] revealed that the discharge of subsurface flow greatly changes with environmental conditions and topography. Meanwhile, Woods et al. [27] proposed a new topographic index along with the field data to calculate the spatial distribution of groundwater flow and the thickness of the saturated zone.
On the other hand, some researchers have investigated groundwater flow and subsurface flow through numerical methods and analytical approaches. Childs [4] first considered the bottom slope of an aquifer and presented a generalized Boussinesq equation for groundwater level in an unconfined sloping aquifer. Then, a bivariate quadrature function was presented by Evans [6] to signify distinct types of topography in a watershed and the terrain analysis was incorporated with slope mapping. Later, the nonlinear Boussinesq equation was linearized and solved analytically by Brutsaert [2] for groundwater table in an unconfined aquifer. This analytical solution suggests a significant framework to explore slope features and hydraulic characteristic of a sloping aquifer. Fan and Bras [7] combined the principle of mass conservation and Darcy's formula to acquire the continuity equation for groundwater flow, and then analytically solve the equation using the method of characteristics (MOC). Referring to Evans [6], Troch et al. [21] categorized the conventional hillslope patterns into nine hillslopes and also used MOC to solve the hillslope-storage equation via the kinematic wave model. Troch et al. [22] again solved many versions of the water storage equation after linearizing and simplifying it. They employed the multistep solver to deal with time and the finite difference method to discretize the space. Shortly, Troch et al. [23] considered the width variation by an exponential function in the linearized Boussinesq equation for different hillslopes and derived the analytical solutions employing Laplace transforms; then, they verified the analytical results with simulated results from a numerical model assuming a constant recharge rate. Dralle et al. [5] studied the impact of spatiotemporally varying recharge on groundwater levels. Their results show that spatially varying recharge may significantly alter groundwater storage variability. In addition, the interactive effects of different recharge patterns and slope morphology on water storage and baseflow recession characteristics were explored. Overall, they revealed that non-uniform recharge can reduce or enhance hydrological influences of hillslope morphology. Recently, Sahoo et al. [16] discussed the dynamics of subsurface water storage in semi-gauged and ungauged catchment areas by building a finite difference model using the hillslope-storage Boussinesq equation using a multistep ordinary differential equation solver in time. In their study, the bedrock leakage flux was considered and estimated by Darcy's equation in an unconfined aquifer. Ranjram and Craig [15] built a Proxy model employing the nonlinear hillslope-storage equation to investigate the aggregate hillslope subsurface flow under uniform recharge in the Dennis Creek watershed. When compared with the numerical solutions, the quality of the Proxy model is proved to be excellent and it could be an efficient tool for further research under quite general conditions.
By and large, all the historical studies have discovered geology is a crucial factor in hill storage, whereas most studies on groundwater storage variability have only considered constant rainfall recharge except that the temporally-varying rainfall recharge was considered in the work of Wu and Hsieh [28]. Nevertheless, Wu and Hsieh [28] didn't consider any type of rainfall recharge and soil texture effect on groundwater storage. In their research, only an analytical solution for groundwater level change was presented and the fixed width of a groundwater aquifer was assumed to be one unit, whereas both analytical and numerical solutions for the change of groundwater storage capacity in variable width aquifers under different scenarios will be proposed in this research. Hence, this study used the hill-storage Boussinesq equation to account for water storage in hillslopes, including uniform, convergent, and divergent hillslopes, but considered temporally varied recharge sources including six patterns of recharge distribution to fit natural rainfall recharge conditions. Receiving a variable, rather than constant, rainfall recharge with time on hillslope storage is more feasible and realistic. Further, the water storage variability in the hillslopes for five soils was also discussed. The presented analytical solution to the linearized hill-storage equation was derived employing the integral transform method. The numerical model using the nonlinear hill-storage equation was built through the method of finite differences. The division in time was carried out by the Runge-Kutta scheme, whereas discretization of space was performed by the upwind scheme and central difference.

Mathematical formulation
A schematic of an unconfined sloping aquifer with surface recharge overlying an impervious base with an inclined angle of the hillslope θ is shown in Figure 1. The sloping aquifer with variable width and a constant thickness was assumed to be isotropic, homogeneous, and saturated. The only drainage of groundwater flow is a ditch at the outlet. Besides, there is no vegetation on the ground.

Governing equation
According to the principle of mass conservation, the following continuity equation for groundwater flow considering any temporally varied recharge distribution and a variable hillslope width can be obtained , and R is temporally varied but spatially uniform rainfall recharge [LT −1 ] in a small catchment. The variable hillslope width was considered as an exponential function, following the work of Troch et al. [23] to signify three different types of hillslopes.
in which a is a coefficient [L −1 ] and c is the hillslope width at the outlet [L]. Note that the hillslope pattern is uniform for a = 0, convergent for a> 0, and divergent for a< 0. According to Darcy's law, the flow discharge was derived: where s ≈ n · h · w, n and θ are effective or drainable porosity and inclined angle of the hillslope, h is average water depth [L], and k p is saturated hydraulic conductivity The abovementioned h and R are defined as where R k is the amount of a recharge rate accompanied with the Heaviside function U (−), being the k th digital value of collected data within the time step ∆t = t k − t k−1 , and M is the total number of collected data.
To solve the nonlinear equation (2.4) analytically is difficult. Thereby, referring to Troch et al. [2], the following linearization relation was introduced:

Initial condition
The initial variation of hillslope storage was presumed: in which γ is initial constant water depth [L] and 0 ≤ γ ≤ D.

Boundary conditions
No influx at the upstream boundary: Besides, no store water at the outlet:

Analytical model
To solve the governing equation (2.8) associated with the initial condition (2.9), and the boundary conditions (2.12) and (2.13), the generalized integral transform (G.I.T.) method proposed byÖzisik [14] was employed Integral transform: Inverse transform: where H is the transformed function of H and K (β m , x) is the kernel. Before solving the foregoing boundary value problem, equation (2.8) was rewritten as where A = kpbD n cos θ, B = kp n sin θ − akpbD n cos θ. Next, during the processes of analytical solution derivation, the first-order differential term of space needs to be eliminated to conform to the form of the heat equation which can be solved by the generalized integral transform with the reference ofÖzisik [14] so that the following change-of-variable technique, i.e. Equation (2.17), was presented.
in which s * is the transformed variable of s, and inserting equation (2.17) into equations (2.16), (2.9), (2.12), and (2.13) gives 4A t . Taking the integral transform (2.14) of equations (2.18) to (2.21) results in where β m is the positive root of the following characteristic equation with m being a natural number: Specifically, the eigenvalue β m is influenced by the average depth D, the parameter a, the fitting parameter b, and the slope θ.
Solving the first-order ordinary differential equation ( Taking inverse transform (2.15) of s yields Then, s can be obtained via the relation (2.17): After acquiring the water storage, the following variables such as average water depth (h), flow rate (Q), outflow discharge (q), and relative water storage (s r ) were sequentially derived: Figure 2. Schematic of mesh for the numerical model.
Here, the analytical solutions were obtained by employing the integral transforms rather than the Laplace transforms because the convergence of the present solutions is better as discussed in Wu and Hsieh [28].

Numerical model
Besides analytically solving the linearized equation (2.8), we solved the nonlinear equation (2.4) by building a finite difference model. In the numerical model, the third-order total variation diminishing (TVD) Runge-Kutta method reported in Shu [18] was used to discretize time, and the upwind scheme and central difference were applied to perform the space discretization with reference to Swanson and Turke [20]. Figure 2 illustrates that the domain was divided into m segments with an equal interval of ∆x by m + 1 nodes along the x direction. Note that the first node (i = 1) and the last node (i = m+1) are virtual outside the domain. Accordingly, the difference equation of hillslope storage for space discretization becomes where s α is the solution of order α, i is the node number, and j is time.
In addition, the initial condition (2.9) becomes Eliminating s α j (0) and s α j (0) with boundary condition s α j (0) = 0, we can obtain With respect to the time discretization, the TVD Runge-Kutta scheme results in the following: where s 0 is the initial condition, and s 1 , s 2 , and s 3 are the solutions of each order. Φ (s) = ∂s ∂t . The final difference equations are appended in the Appendix.

Validation of both mathematical models
The parameter values L = 100 m, θ=0.05 (in radian), b = 1, n = 0.3, k p = 1 mh −1 , D = 2 m, γ = 0.4 m, and R = 10 mmd −1 adopted by Troch et al. [23] were employed to verify the proposed analytical solution. The parameters controlling the shape and width for three hillslope types were (1) uniform: a= 0 and c= 21.627 m, (2) convergent: a= 0.02 m −1 and c= 6.77 m, and (3) divergent: a=−0.02 m −1 and c= 50.024 m. All three types of hillslopes have the same projected areas and received the consistent rainfall recharge. The initial water storage was also presumed to be identical for all the three types of hillslopes. Figure 3 generally shows the spatial distribution of groundwater tables in the different hillslopes for 1, 5, and 20 days. The figure reveals that the present analytical results using the integral transform are consistent with that obtained by Laplace transform, i.e., equation (33) in Troch et al. [23], thus validating the present analytical model. Nevertheless, the Laplace transform method cannot be used directly in this study for two causes. First, finding the function of the inverse Laplace transform is challenging, and even found, it is usually too complicated. Second, the solution convergence derived by our analytical approach was better than that acquired by Laplace transforms. The detail was discussed in the work of Wu and Hsieh [28] and thus omitted herein. As reported in Verhoest and Troch [24], the first 999 terms of their series expressions need to be summed, whereas the present solution only needs summation of the first (usually less than) 50 terms to reach convergence. The convergence speed of the present solution is faster. The simulation time of each case, number of terms to convergence and the criteria of convergence of the present analytical solutions for three types of hillslopes are listed in Table 1.
To verify the present numerical model, the following parameters were adopted for uniform hillslope: a = 0 and c = 50 m, L = 100 m, D = 2 m, n = 0.3, k p = 1 mh −1 , θ = 0.05 (in radian), γ = 0, and R = 10 mmd −1   with reference to the work of Troch et al. [22]. Figure 4 illustrates the variation in relative storage and the results are very consistent with the previous results, which validates the present finite difference model.

Comparison of both solutions
Because the analytical and numerical solutions for the same linearized equation (2.8) are completely overlapped, the comparison of water storage between both solutions' results under three types of hillslopes was omitted. Figure 5 demonstrates the modelling results for convergent hillslope for various durations with D = 2m and γ = 1m. In the figure, the parameter b was adjusted for better modelling results for each case. When (R, b) = (50, 0.5 − 0.7) in Figure 5, an apparent gap was observed between these two solutions. The gap was estimated through an error analysis denoted by an averaged absolute relative percentage difference (RPD). The RPDs were 3.93% as b = 0.7 and 2.78% as b = 0.2 in Figure 5a, 8.72% as b = 0.7 and 16.09% as b = 0.5 in Figure 5b, and 23.76% as b = 0.7 and 35.49% as b = 0.5 in Figure 5c. However, it should be noted that the RPDs were small except the cases in Figure 5c. These results signify the gap increased with an increasing duration even while an optimal value of b was searched. This implies that the parameter b must be tuned for different Table 2. Elapsed time of performance of calculating relative water storage (Fig. 5) Table 2 with the computer specification: Intel(R) Core (TM) i7-4790K CPU @ 4.00GHz, 8 GB of RAM by using MATLAB.
Although Table 2 indicates that the present numerical solution dominates over the analytical solution, there is slight difference in computing time between both solutions. The linearized analytical solution was derived for the linearized Boussinesq equation, whereas the numerical solution was obtained by solving the nonlinear Boussinesq equation. To improve the above shortcoming of analytical solutions, we may employ the perturbation technique to analytically solve the nonlinear Boussinesq equation in future research, which is a very tough challenge. Nevertheless, the comparison between the present analytical solution and previous analytical solutions agrees very well in this study, thus validating the present analytical solution. Further, analytical solutions can help clarify the steps in the research processes and benefit the verification of new developed numerical models. In this study, the difference of computing time between analytical solutions and numerical solutions is almost negligible by means of modern computers.
In summary, the foregoing results reveal that the proposed analytical solution is strongly sensitive to the parameter b which pertains to hill storage, hillslope thickness, and hillslope width. Thereby, tuning b for distinct cases can yield the analytical results closer to the numerical ones. The optimal value of b was determined through the trial and error technique.

Variation of remaining hill storage
To calculate the remaining water storage left in the hillslope after drainage in response to the temporally varied recharge for the three types of hillslopes with different slopes, a parameter of dimensionless remaining content of stored water denoted by s was defined as  Figures 6-8. As expected, the storage water drained fastest on the steepest divergent hillslopes. These figures also display that the remaining hill storage decreased when the simulation duration and slope increased. In short, the reduction rate of water storage rose for steep slopes, particularly for divergent hillslope due to the quick and massive drainage initially.
Using the values of hydraulic conductivity and drainage porosity of the work of Sahoo and Sahoo [17], the variation of s at various durations for different hillslopes and soil types/aquifer settings is illustrated in Figure 9.  -----,,----,--..,....----,.--...... -"T""-.-----r-    When the soil composition of the aquifer is gravel (k p = 60 m/d, n = 0.33), the relative water storage capacity is mostly below 50% because of quick drainage, and the storage capacity of the divergent hillslope decreases most fast, especially for steep slope as shown in Figure 9a. When the aquifer is composed of sand (k p = 25 m/d, n = 0.4), the storage capacity is more than 30% at the duration of 1 day for all hillslopes. However, at the duration of 3 days, the storage capacity decreases by 10-25%. At the duration of 5 days, the storage capacity is reduced by about 25-40% when compared to the first day as shown in Figure 9b. When the aquifer soil is composed of sandy loam (k p = 12 m/d, n = 0.32), the storage capacity is mostly more than 50% at the duration of 1 day for all hillslopes. At the duration of 3 days and the θ below 55%, the storage capacity will be reduced by 30-50%. At the duration of 5 days, the storage capacity is reduced by about 50-70% when compared to the first day as shown in Figure 9c. When the aquifer is composed of silt and clay, k p = 0.7 and 0.01 m/d and n = 0.5 and 0.62, respectively, water is extremely difficult to flow in the aquifer. When θ is greater than 55%, there is a slight impact on the storage capacity, especially for the divergent hillslope after 3 days. The storage capacity can reach 85-95% for both silt and clay as shown in Figure 9d and e, respectively. Based on the aforementioned results, the storage capacity for the same hillslope, duration, and slope of different aquifer compositions follows the order: clay > silt > sandy loam > sand > gravel, and thus the different drainage flow rates follow the order: gravel > sand > sandy loam > silt > clay, implying that hydraulic permeability and slope are two major factors affecting the water storage capacity. Such a result is consistent with that of Sahoo and Sahoo [17].

Effect of variable recharge rates
In nature the rainfall recharge is usually nonuniformly distributed, so temporal variation of it was considered. Referring to the works of Wartalska et al. [25] and Shih [19], six recharge patterns over time shown in Figure 10 were introduced to discuss their effects on hillslope storage. Assuming the same aquifer and groundwater conditions and setting D = 5 m, θ = 5%, and γ = 0 for the simulated scenarios, Figure 11 shows the temporal variation of the water level at various durations under six recharge patterns for convergent hillslope. In the initial stage of the simulation, groundwater level will vary according to the amount of rainfall recharge. After a period of continuous recharge, water level and outflow will eventually reach an equilibrium, and this balance condition will vary with different soil textures and configurations. As illustrated in Figure 11, the distribution of water level change within a recharge period of 24 h was highly affected by the recharge pattern, but the difference of water level distribution among all the recharge types became very slight after the cease of recharge. Similar effects were discovered for uniform and divergent hillslopes, and thus they were omitted. Figure 11 further shows the distribution of water level for the type of peak at the first section was the largest at the early stage because of the initial abundant recharging water, and conversely the water level for the type of peak at the last section was the smallest in the early stage. When the simulation time is more than 24 hours, the distribution of all recharge patterns almost overlapped because of no more recharge. Figure 12 demonstrates the hydrograph of outflow under six recharge patterns with the same accumulative recharge amount for the three types of hillslopes, in which three more recharge patterns, i.e., peak in the first section, peak in the last section, and double peak were added. Figure 12a illustrates that for the convergent hillslope, the maximum peak among these outflows occurred at the pattern of peak in the last section, and finally all outflows increasingly approached one value. Similar results in Figure 12b and c were also observed for uniform and divergent hillslopes separately, except that all outflows decreasingly approached one value. Specifically, at the outlet, the cross-sectional area was relatively small for the convergent shape, so the outflow was the lowest. In contrast, the outflow was the highest because the cross-sectional area was relatively large for the divergent shape. These hydraulic effects are signified in Figure 12a and c. Figure 12 also indicates that for convergent hillslope, while the recharge terminates, the outflow reduces and then gradually rises for a long time owing to the slow water release near the upstream boundary where the storage water is relatively abundant because of the widest width. Moreover, the outflow decreases gradually while the recharge ceases for uniform hillslope, but it drops more quickly because of the rapid release of stored water for divergent hillslope.  . The remaining storage ratio versus slope at various durations for different hillslopes and different soil types/aquifer settings. Note that con., uni., and div. mean convergent, uniform, and divergent, respectively.

Conclusions
Both analytical and numerical models were developed to simulate the variation of water storage response to any temporally variable recharge for three types of hillslopes. The numerical model was built based on the nonlinear equation but the analytical solution was derived for the linearized equation.
The present analytical solution was derived using the generalized integral transforms instead of Laplace transforms to avoid the difficulty which inverse Laplace transform might encounter. The analytical results were consistent with the results of Troch et al. [21] for convergent, uniform, and divergent hillslopes. On the other hand, our numerical results agreed well with the results obtained through a numerical integration technique in the work of Troch et al. [7].
The results showed that hillslope pattern and slope greatly affect the drainage of hillslope water storage. It is easily deduced that drainage of the storage was the slowest for convergent hillslope and the fastest for divergent hillslope for the same recharge distribution and slope. Comparison of both solutions reveals that the discrepancy between the results decreases with a falling recharge rate for convergent hillslope. Further, the distributions of groundwater level and outflow discharge were highly affected by the time-varying recharge pattern within the    recharging period, but they almost became unaffected by all the recharge types after the cease of recharge for all the three types of hillslopes. The water levels of the types of peak at the first section and peak of the first quarter section were the largest at the early stage because of the initial abundant recharging water.
To sum up, the present analytical model for groundwater storage and water level discussed the effect of six different patterns of temporally varying rainfall recharge and five soils, respectively, which was not considered in previous research. The convergence speed of the analytical solutions derived by the generalized integral transforms was faster than that of previous research by Laplace transforms. The present analytical solutions can help clarify the steps in the research processes and benefit the verification of more complex analytical and numerical models developed in the future.