A fractional-order infiltration model to improve the simulation of rainfall/runoff in combination with a 2D shallow water model

Hydrologic modeling is one of the most important fields of research within environment-related disciplines. An accurate simulation of rainfall/runoff processes allows the prediction of flood risks in natural or urban environments and a more efficient water resources management.

The use of distributed models for hydrologic simulation provides detailed computation of the spatial variations of the variables of interest within the domain, such as water depth, flow velocity in two directions or even infiltration rate. Recent work on this topic can be found in Xia et al. (2017), Yu & Duan (2017), Bellos & Tsakiris (2016), Cea & Bladé (2015), Liang et al. (2015), Costabile et al. (2013, 2012) and Cea et al. (2010), where the applications cover a range of test cases and practical cases of different degrees of difficulty. This is of special relevance in abrupt terrain topography, characteristic of mountain river catchments, such as the Araguás and Arnás catchments that we will consider later in this paper. This information is averaged when using lumped simulation models (Caviedes-Voullième et al. 2012; Fernández-Pato et al. 2016).

The numerical stability in unsteady shallow water (SW) flow simulation has been a matter of recent research (Murillo & García-Navarro 2010) due to the required modifications of the basic scheme to real situations. The size of the allowable time step to ensure stability is crucial when applying explicit methods for hydrological purposes where dynamically wet/dry fronts are established. In overland hydrological problems, wet/dry frontiers may be found everywhere in combination with dominant friction and bed slope. These wetting and drying situations over all the cells imply that special care must be taken when calculating wet/dry fronts (Murillo & García-Navarro 2010).

In order to include the computation of water losses by infiltration in rainfall/runoff simulation models, where this is the dominant sink term, there is a wide range of infiltration models, such as Green-Ampt (GA) (Green & Ampt 1911), Horton (Horton 1933) and Philip (Philip 1969). The combination of a distributed surface flow model with infiltration models has been frequently applied to river catchment simulation (Esteves et al. 2000; Fiedler & Ramírez 2000; López-Barrera et al. 2011; Caviedes-Voullième et al. 2012; Costabile et al. 2013; Singh 2015; Bellos & Tsakiris 2016; Fernández-Pato et al. 2016). Some authors formulate the distributed surface flow model as a kinematic or diffusive wave simplification (Luo 2007; Cea et al. 2010; López-Barrera et al. 2011; Costabile et al. 2012; Langhans et al. 2014). Other authors (Caviedes-Voullième et al. 2012; Simons et al. 2014) use the fully dynamic wave model (a.k.a. shallow water (SW)) in order to avoid errors when simulating surface flows with non-negligible accelerations. In Fernández-Pato et al. (2016), a distributed two-dimensional (2D) fully dynamic SW model together with the Green-Ampt infiltration model is applied to real rainfall/runoff situations in the Arnás river catchment. The agreement between numerical simulations and overland flow data in that basin (moderate/high slopes and with a high infiltration rate) still has the capacity for further improvement. In this paper, we explore the use of a fractional-order Green-Ampt infiltration model (FOGA) as a means to improve the hydrograph fitting.

In recent years, fractional-order derivatives have been applied to hydrological modeling (Borthwick 2010; Su 2010, 2012; Benson et al. 2013), water movement in soils (Pachepsky et al. 2003; Su 2014) and solute transport modeling in both overland (Deng et al. 2004, 2006) and subsurface flows (Martinez et al. 2010; Sun et al. 2014), since the classical advection-dispersion equation seems to fail to capture some important solute transport features.

Fractional derivatives, unlike ordinary derivatives, provide an excellent instrument for the description of memory and hereditary complex processes. For this reason, although fractional calculus dates from the 17th century, it has recently become a research area of growing interest due mainly to the ever-widening range of applications in physics, engineering, chemistry, biology, and economics. We can refer to Kilbas et al. (2006), Machado et al. (2011), Metzler & Klafter (2000), Oldham & Spanier (1974) and Samko et al. (1993) and the references therein for a general overview of this research area. In particular, fractional calculus is a key tool in the study of anomalous diffusion: subdiffusion and superdiffusion processes. These processes are described in detail in the excellent report by Metzler & Klafter (2000) and the authors focus the discussion on the case that the mean square displacement of the growth of the particles has a power-law pattern in the course of time, unlike the Brownian motion.

Regarding the infiltration calculation, a fractional-order Green-Ampt infiltration model, which is a generalization of the classical model, was proposed in recent literature for a better fitting to experimental data (Gerasimov et al. 2010; Voller 2011). Namely, in Voller (2011), a FOGA model is proposed for predicting the infiltration rates into columns of non-homogeneous soil. This curve can exhibit a non-monotonic behavior, which facilitates soil infiltration recovery during dry periods. The author also reports a superdiffusion behavior of the infiltration front at the early times and this motivates the use of a fractional derivative in the infiltration model.

In the present work, the 2D shallow water equations (SWE), which govern the unsteady surface flow, are discretized with the finite volume numerical scheme designed in Murillo et al. (2007) and Murillo & García-Navarro (2010). The novelty in the present work is that the FOGA method by Voller (2011) is used as a mass source term for a better modeling of the infiltration process in complex real catchments where infiltration is the major sink term. Two basins in which the traditional formulation of infiltration methods does not allow a correct fitting to the observed outlet hydrographs (López-Barrera et al. 2011; Caviedes-Voullième et al. 2012; Fernández-Pato et al. 2016) are simulated using the FOGA model. In particular, a considerable delay in rising limbs of the numerical hydrographs with respect to the experimental data is observed when classical laws are used for estimating the infiltration losses. The main objective of this paper is to show how the modification of the infiltration Green-Ampt method by means of fractional calculus significantly improves both rising and falling limbs. We will show an enhancement in the numerical results for two real catchments if the order of the fractional derivative is variable in time and distributed in space.

The structure of the paper is as follows: first a brief overview of the surface flow equations and the classical and fractional-order Green-Ampt infiltration models are presented. Several numerical tests are then given for a better understanding of the FOGA model. The improvements of the fractional infiltration model compared to the classical model are then shown, by means of the calibration of three real storm events in two different river basins.

The model considered in this paper uses the 2D SWE for the overland flow and the Green-Ampt law for the infiltration processes. Both are described in the next sections and the latter in more detail since a generalized version of that law is obtained using a fractional Darcy's law.

Assuming dominant advection, (1) can be classified as a hyperbolic system. Its solution is approximated with the well-balanced explicit, first-order, upwind finite volume scheme described in Murillo & García-Navarro (2010). The wet/dry fronts are well tracked providing stable solutions due to the use of dynamical control of the time step size with a numerical mass error of the order of the machine precision. The use of a distributed surface flow model allows the calculation of all the hydraulic and hydrologic variables, such as the water depth h, the flow velocities or the infiltration rate f, in every cell of the computational mesh.

Once the function s has been obtained, the infiltration capacity comes from (13) (recall that ⁠) and the cumulative infiltration is evaluated as ⁠. It is worth noting here that the infiltration capacity and cumulative infiltration are computed at each cell of the computational domain so that and (Fernández-Pato et al. 2016).

In this section, the generalized Green-Ampt infiltration model is derived by using the fractional Darcy's law given in (15). Unlike Voller (2011), the variables and parameters have physical dimensions because the infiltration rate f will be used in the simulations of rainfall/runoff in real river catchments.

From (21), one observes that has dimensions and therefore the dimensions of must be so that the subsurface hydraulic flux q has the correct physical dimensions ⁠.

Note also that, as in the GA model, both the infiltration rate and the cumulative infiltration are computed per cell in our model and, therefore, they are spatially distributed in the horizontal plane.

Table 1 shows a summary of the intermediate equations for both GA and FOGA models. Note that if ⁠, both GA and FOGA models are governed by the same differential equation.

Several numerical experiments are presented to illustrate the effect of the order of the Caputo fractional derivative on the FOGA model. Different conditions of surface water availability are considered but the same parameters ⁠, ⁠, of this model are used in all cases, unless otherwise specified.

Figure 2 shows the temporal evolution of all the variables of interest: rainfall, infiltration capacity, infiltration rate, infiltration volume, rainfall volume and surface water. For illustration purposes, two choices of have been considered: and ⁠.

By comparing the plots corresponding to and in Figure 2, it is observed that when lowering the value of ⁠, the infiltration rate globally decreases for all times. This conclusion is consistent with the results obtained in Voller (2011) for a continuously ponded soil.

A plane with slope 0.005 is considered in this example and the dimensions of the domain are 2,000 m × 20 m × 10 m (see Figure 3). Manning's roughness coefficient is set to n = 0.02 sm^−1/3. The initial conditions for the surface water equations are zero water depth and zero discharge everywhere, i.e., dry surface conditions. Water enters the domain only through rainfall, which is assumed to be constant in space, hence there are no inlet boundaries. The only open boundary is at the outlet (downslope) and free outflow is assumed. In order to study how the order of the Caputo fractional derivative affects the shape of the outlet hydrograph, seven test cases have been designed (see Table 2). The values of the hydraulic conductivity and the rainfall rate R corresponding to Case 3.1 are taken as reference. In Cases 3.2 and 3.3, the values of and R are multiplied by 10. In all the cases, the simulations are performed for the values of ⁠.

Figures 4–6 show the numerical outlet discharge Q for all the proposed cases. This value is computed from the integration of the outlet unit discharges predicted by the surface model at the outlet boundary. In light of the results, several conclusions are reached. For Case 3.1, the lower the order ⁠, the lower the soil infiltration rate, since the outlet hydrographs are showing higher peak discharge values. The same conclusion is achieved in Case 3.2, where the value of is increased by a factor of 10. Nevertheless, a change in the trend of the outlet hydrographs is observed for certain values of ⁠, leading to non-monotonic curves. This behavior is a consequence of non-monotonic infiltration rate curves, generated by the FOGA model (Voller 2011), for certain values of and the classical parameters set ⁠.

In order to test the influence of the rainfall rate on the infiltration rate recovery predicted by the FOGA model, Case 3.3 is proposed. By maintaining the same value as in Case 3.2, the rainfall rate is increased by a factor of 10. As seen in Figure 4(e) and 4(f), the same trend as in Case 3.2 for the outlet hydrographs is reached. This shows that this change on the curves trend depends only on the value, not on the ratio between and R.

The other Green-Ampt parameters, and have been also modified in Cases 3.4 and 3.5 observing the same change in the trend, to a lesser extent, in the outlet hydrographs (Figure 5).

Figure 6 shows the results corresponding to Case 3.6 in which an unsteady rainfall pattern is considered and several values of have been tested, ranging from 0.35 to 1. As observed in the previous cases, the smaller the value of ⁠, the lower the soil infiltration capacity and, hence, the higher the peak discharges at the outlet. In general terms, it seems that, for this particular case, the order of the fractional derivative has a global effect, raising or lowering the full hydrograph.

In this section, two real catchments (Araguás and Arnás) of the Northern Spanish Pyrenees are simulated. The infiltration parameters are properly calibrated in order to fit observed outlet hydrographs, corresponding to three different storm events. The classical and fractional-order infiltration models are compared and the improvement on the numerical prediction of the fractional model is shown.

The Araguás catchment is located in the Central Pyrenees (Figure 7(a)) and it has an extension of 0.45 km² (García-Ruiz et al. 2010). Its altitude ranges from 780 to 1,100 m.a.s.l. and the mean slope varies from 20% to 43%. Due to the small size of the basin, a constant Manning's roughness coefficient of n = 0.025 sm^−1/3 is set. A triangular unstructured mesh of 7,728 cells is used for the spatial discretization of the catchment (Figure 7(b)).

In this catchment, a single event (which is referred to as Event 1) is considered. The flow discharge measurements were acquired at the outlet of the basin (see Figure 7(a)) with a 5-minute frequency. Rainfall was registered by a rain gauge also with the same frequency. Figure 8 shows the observed hyetograph and the outlet hydrograph for this particular storm event.

In this section, the numerical results obtained for the Araguás catchment are presented. In order to highlight the capacities of distributed models, Figure 9 shows the spatial distribution of the numerical values of water depth, flow velocity and cumulative infiltration at t = 13,500 s. It can be seen that the main channel cumulative infiltration values are significantly different from the ones computed in the hillsides. This detailed computation of the hydraulic and hydrologic variables leads to a better prediction of the outlet hydrographs and, therefore, to a better fit of the observed data.

Table 3 summarizes the simulated cases set with a constant distribution of ⁠. Case 4.1 is considered as the reference case. It corresponds to the parameter set that provides the best fit using the classical GA model ⁠. Cases 4.2 and 4.3 keep the parameters K, and as in Case 4.1 but the order of the fractional derivative is modified. Figure 10 shows the outlet hydrographs corresponding to Cases 4.1, 4.2 and 4.3. As seen on the previous numerical examples, the effect of decreasing corresponds to a global reduction of the soil infiltration capacity predicted by the GA model. Hence, the lower the value of the greater the value of the outlet discharge peak.

It is also interesting to explore the possibility of reproducing the same outlet hydrograph of Case 4.1 by means of several combinations of the infiltration parameters with ⁠. As shown in Figure 11, Cases 4.4 to 4.9 generate very similar hydrographs to Case 4.1 by modifying one single parameter each time for two different values of ⁠. In particular, we consider the values of α = 0.7, 0.9. Thus, the numerical results suggest that the fitting cannot be improved in this event using a single value for ⁠.

In light of the previous numerical results, the GA model seems to predict an excessive soil infiltration capacity at the beginning of the storm when the soil is almost dry.

The particular function proposed (24) discriminates between the presence or absence of surface water but it depends weakly on the surface water depth as it gets quickly truncated (Figure 12). The main purpose of this essentially phenomenological fitting model is to provide a way to determine the order of the fractional derivative at the earlier stages of infiltration. Hence, it does not formulate any physical process. Once the infiltration depth and the surface water depth exceed a threshold, the model reduces to the classical GA formulation. Table 4 gives the FOGA model parameters obtained for the storm Event 1 considered in this catchment, as well as the classical GA model parameters which provide the best fit. It also presents the values for both GA and FOGA models, computed as ⁠, where and are the computed and experimental discharges, respectively, and N is the number of discharge curve points. The relative difference between the error produced by both models is also shown. For this event, the FOGA model fitting error is 29.3% lower than the one produced by the GA model. Figure 13 shows the numerical hydrograph obtained with both the classical GA and FOGA models while the distributed values of h, F and for this event are given in Figure 14. The FOGA model reproduces the hydrograph better. The arrival time is better fitted than with the GA model although the whole physical process is not correctly simulated because the receding part of the curve is lower than the observed one.

The Arnás catchment (2.84 km², 900–1,340 m.a.s.l.) is located in the Northern Spanish Pyrenees (see Figure 15(a)). Geologically, the catchment and its land use have suffered several changes in recent years which have significantly modified the vegetation cover. This includes patches of forest, grassland meadows, dense bush areas and bare land. The soil types and vegetation mapping have been widely studied in Lana-Renault (2007), Lana-Renault et al. (2007) and Serrano-Pacheco (2009). All these maps conform an accurate hydrological characterization of the Arnás watershed.

Numerical simulations of the Arnás catchment was carried out in López-Barrera et al. (2011) by means of a 2D diffusion wave model for the surface flow and both the Horton and Green-Ampt models for water losses due to infiltration. In Caviedes-Voullième et al. (2012), a combination of the 2D SWE with Soil Conservation Service (SCS)-Curve Number model for the precipitation losses was used for rainfall/runoff simulation in this catchment. In both cases, the numerical results showed a poor agreement with the observed data. A more careful calibration was reported in Fernández-Pato et al. (2016), where a sensitivity analysis was performed in order to identify the influence of the topography and the storm characteristics on the Horton and Green-Ampt infiltration models. This led to a significant improvement in the fitting of the outlet hydrographs for several storm events.

In this watershed, the catchment topography meshing has been widely studied in Caviedes-Voullième et al. (2012), where the authors found the optimal mesh for solving the SWE in order to minimize the computational time without losing quality in the numerical results. All the simulated events for this catchment use this optimal mesh (see Figure 15(c)).

In this basin, two events (which are referred to as Event 2 and 3) are simulated and compared with the observed data. In both cases, discharge measurements were taken at the outlet (see Figure 15(a)) with a frequency of 5 minutes. On the other hand, rainfall intensity was registered by a rain gauge with a 5-minute frequency for Event 2 and 60-minute frequency for Event 3. Figure 16 shows the observed hyetographs and outlet hydrographs for both storm events.

The calibrations obtained for Events 2 and 3 are presented in this section. Tables 5 and 6 summarize the parameters for all the considered events in this catchment using the classical GA model and the FOGA infiltration model based on a variable ⁠. The hydrograph fittings are shown in Figures 17 and 18 for Events 2 and 3, respectively. In both cases, significant improvements over the classical GA model are observed. Hydrograph rising limbs show a better adjustment to the observed data in both events. Tables 5 and 6 also presents the values for both GA and FOGA model and the relative difference between the error produced by both models. In this catchment, the FOGA model fitting error is 16.7% lower than the one produced by the GA model for Event 2 and 25.8% lower for Event 3. The error values obtained for these events (together with the one obtained for the Event 1) quantify the improvement in the hydrograph fitting by the FOGA model. As in the previous section, spatial distributions of h, F and are presented in Figure 19 for Event 2 at t = 15,000 s.

In this work, a generalized Green-Ampt infiltration law has been combined with a 2D full shallow water model and applied to real catchment simulation. The water losses due to infiltration have been calculated in a spatially and temporally distributed way in order to take advantage of the detailed topographic information. The Green-Ampt infiltration law has been generalized by using a Caputo fractional derivative in Darcy's law for the subsurface hydraulic flux. The order of the fractional derivative has been treated as a new parameter of the infiltration model, which has been properly estimated.

In the light of the numerical results, the conclusions of this paper can be summarized in the following points:

• The numerical results suggest that the fitting cannot be improved in field basin events using a constant and uniform value for in the whole computational domain.

• The performance of the method in the presence of long and complex rainfall events has led to a formulation of the order of the fractional derivative as a combined function of the cumulative infiltration F and the surface water depth h in each cell of the computational domain.

• With a proper estimation of ⁠, a significant improvement in the fitting of the outlet hydrograph to the experimental data is observed in all the tested real cases. The rising and falling limbs of the hydrograph are successfully predicted with this model, removing the delays in the peak discharge time reported in previous publications.

In general terms, the application of a FOGA infiltration model to real catchments calibration significantly improves previous published results, leading to promising future research opportunities in cases where the sink term is the major mechanism to fit surface runoff. Nevertheless, a physical justification for this fractional model with variable order is still an open research question.

The present work was partially funded by the Aragón Government through the Fondo Social Europeo. This research has also been supported by the Research Projects CGL2015-66114-R funded by the Spanish Ministry of Economy and Competitiveness (MINECO) and MTM2016-75139-R funded by MINECO. The corresponding author also wants to thank the MINECO and Hydronia L.C.C. for his Research Grant DI-14-06987.