Abstract

In this work, a distributed two-dimensional (2D) shallow water (SW) flow model is combined with a fractional-order version of the Green-Ampt (FOGA) infiltration law to improve rainfall/runoff simulation in real catchments. The surface water model is based on a robust finite volume method on triangular grids that can handle flow over dry bed and multiple wet/dry fronts. When supplied with adequate infiltration laws, this model can provide useful information in surface hydrology. The classical Green-Ampt law is generalized by using a Caputo fractional derivative of order less than or equal to 1 in Darcy's law. The novelty of this combination is that, on the one hand, the distributed SW simulation provides a detailed surface water distribution and, on the other hand, the FOGA model offers the possibility to model infiltration rates not monotonically decreasing. In order to obtain the best results, a non-uniform order of the fractional derivative depending on the cumulative infiltration and the existence of available surface water is proposed for realistic cases. This allows significant improvement of previous published numerical results in the literature for several storm events in catchments where the infiltration process occurs.

INTRODUCTION

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.

MATHEMATICAL MODEL

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.

Surface flow model

In this paper, the surface flow is formulated using the 2D SWE (Vreugdenhil 1994):  
formula
(1)
where  
formula
(2)
is the vector of conserved variables and the superscript T denotes transpose. Here h represents the water depth and and are the unit discharges, with u and v the depth averaged components of the velocity vector along the - and -axes, respectively. The fluxes of the conserved variables can be written as  
formula
(3)
where g is the acceleration due to gravity.
We now describe the three terms in the right-hand side of (1). The term corresponds to friction and it is defined as  
formula
(4)
where are the friction slopes in the - and -directions, respectively. They are written in terms of Manning's roughness coefficient :  
formula
(5)
The term in (1) is defined by  
formula
(6)
and accounts for the pressure force variation along the bottom in - and -directions, formulated in terms of the bed level z slopes.
Finally, the term represents the mass sources/sinks due to rainfall/infiltration:  
formula
(7)
where R is the rainfall intensity and f is the infiltration rate.
When there is available surface water for infiltration , the infiltration rate f coincides with the infiltration capacity , which is predicted by the infiltration model. Otherwise, f is defined by comparing R and and it is given by:  
formula
(8)

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.

Fractional-order Green-Ampt infiltration model

Basic concepts

The classical Green-Ampt model (Green & Ampt 1911; Mein & Larson 1973; Te Chow et al. 1988) assumes a sharp wetting front at the position (see Figure 1) separating the saturated soil region with a water content equal to the porosity from the unsaturated region with an initial water content . Additionally, the water suction at the wetting front, denoted by , is assumed to remain constant. Under these considerations, the vertical hydraulic flux per unit area in the saturated area q is given by Darcy's law  
formula
(9)
where is the hydraulic head and is the saturated hydraulic conductivity.
Figure 1

Graphical representation of the variables of the Green-Ampt infiltration wetting front.

Figure 1

Graphical representation of the variables of the Green-Ampt infiltration wetting front.

The main assumption in the saturated region is:  
formula
(10)
The surface water depth is used to define the upper boundary condition . The lower boundary condition, at the wetting front position, is set as . Hence, the solution of (10) can be easily obtained and it is given by  
formula
(11)
On the other hand, (10) also leads to an integral mass balance equation where the subsurface hydraulic flux q equals the infiltration capacity , which is defined by  
formula
(12)
where . Hence,  
formula
(13)
By combining (9), (11) and (13), the classical Green-Ampt equation is obtained:  
formula
(14)

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 Voller (2011), the classical Green-Ampt infiltration model is generalized by considering the subsurface hydraulic flux as a fractional-order derivative of H. This is justified by means of several empirical results that highlight some deviations of this theoretical model from the field infiltration measurements in heterogeneous media and a generalized Darcy's law is used to improve the numerical results. The derivation provided in Voller (2011) for a 1D vertical problem is repeated here for the sake of completeness in order to show how it has been incorporated into our spatially distributed model. Focusing in every single computational cell, the vertical water movement is characterized by means of the flux:  
formula
(15)
where is the hydraulic conductivity with the proper dimensions for the fractional model , allowing the correct physical dimensions for the hydraulic flux . In addition, denotes the Caputo fractional derivative (Diethelm 2010) of order with , which is defined by  
formula
(16)
and  
formula
(17)
is the Riemann-Liouville fractional integral operator of order , and denotes the Euler's Gamma function. It is well known (Diethelm 2010) that if , then the operator coincides with the partial derivative . Thus, the fractional-order Green-Ampt infiltration law that is described in the next section is a generalization of the classical version of this law.
As the Caputo fractional derivative of the function is given by (Diethelm 2010)  
formula
(18)
it is possible to derive the generalized Green-Ampt infiltration law.

Fractional-order GA model

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.

The balance of the subsurface flux q leads to a simple governing equation for the saturated region , and the function H is the solution of the following two-point boundary value problem with Dirichlet boundary conditions  
formula
(19a)
 
formula
(19b)
 
formula
(19c)
It is worth noting that in general the Caputo fractional derivative lacks the property and thus the solution of problem (19) is not a linear combination of the functions and z (see for example Diethelm (2010)). The solution of problem (19) is given by  
formula
(20)
and therefore it is now a linear combination of and .
Hence, the subsurface hydraulic flux is given by  
formula
(21)
where we used (15) and (18) with and . Replacing the expression of the flux (21) into the mass balance Equation (13), the following initial-value problem for the wetting front is obtained  
formula
(22)
Note that if , then the initial-value problem (22) that governs the wetting front s at each time level coincides with the classical Green-Ampt infiltration model (see (14)) and one would recover the solution of that model. In this sense, the FOGA model generalizes the classical Green-Ampt infiltration model.

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.

Table 1

Summary of GA and FOGA models

GA FOGA 
  
  
  
  
  
  
  
  
  
GA FOGA 
  
  
  
  
  
  
  
  
  

NUMERICAL TESTS

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.

Test 1: unsteady rainfall over horizontal soil

A horizontal soil with closed walls is assumed with unsteady rainfall pattern for the generation of the surface water ( at ) in order to examine the influence of the order of the Caputo fractional derivative on the infiltration curves. This numerical experiment represents a starting point for the application of the FOGA infiltration method to natural storms in real catchments. In this experiment, the total rainfall volume is defined by  
formula
(23)

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 .

Figure 2

Unsteady rainfall: temporal evolution of the hydrological variables for (a) and (b).

Figure 2

Unsteady rainfall: temporal evolution of the hydrological variables for (a) and (b).

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.

Test 2: rainfall on a slope

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 .

Table 2

Parameter setting for Example 3

Test case R (mm/s) Kα (mα/s) ψ (m) Δθ (m3/m3
3.1 0.0825  0.05 0.38 
3.2 0.0825  0.05 0.38 
3.3 0.825  0.05 0.38 
3.4 0.0825  0.005 0.38 
3.5 0.0825  0.05 0.038 
3.6 Unsteady  0.05 0.038 
3.7 Unsteady  0.05 0.38 
Test case R (mm/s) Kα (mα/s) ψ (m) Δθ (m3/m3
3.1 0.0825  0.05 0.38 
3.2 0.0825  0.05 0.38 
3.3 0.825  0.05 0.38 
3.4 0.0825  0.005 0.38 
3.5 0.0825  0.05 0.038 
3.6 Unsteady  0.05 0.038 
3.7 Unsteady  0.05 0.38 
Figure 3

Example 3: topography.

Figure 3

Example 3: topography.

Figures 46 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 .

Figure 4

Rainfall on a slope: outlet hydrographs for Cases 3.1 (a), (b), 3.2 (c), (d), and 3.3 (e), (f). The figures on the right represent a close-up view.

Figure 4

Rainfall on a slope: outlet hydrographs for Cases 3.1 (a), (b), 3.2 (c), (d), and 3.3 (e), (f). The figures on the right represent a close-up view.

Figure 5

Rainfall on a slope: outlet hydrographs for Cases 3.4 (a), (b) and 3.5 (c), (d). The figures on the right represent a close-up view.

Figure 5

Rainfall on a slope: outlet hydrographs for Cases 3.4 (a), (b) and 3.5 (c), (d). The figures on the right represent a close-up view.

Figure 6

Rainfall on a slope: outlet hydrographs for Case 3.6.

Figure 6

Rainfall on a slope: outlet hydrographs for Case 3.6.

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.

APPLICATION TO REAL CATCHMENTS

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.

Araguás catchment

Catchment description and meshing

The Araguás catchment is located in the Central Pyrenees (Figure 7(a)) and it has an extension of 0.45 km2 (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)).

Figure 7

Araguás catchment characteristics. (a) Catchment location and hypsometry map. (b) Computational mesh.

Figure 7

Araguás catchment characteristics. (a) Catchment location and hypsometry map. (b) Computational mesh.

Event description

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.

Figure 8

Experimental hyetograph and hydrograph for Araguás basin.

Figure 8

Experimental hyetograph and hydrograph for Araguás basin.

Numerical results

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.

Figure 9

Event 1: spatial distribution of the water depth , flow velocity and cumulative infiltration at t = 13,500 s.

Figure 9

Event 1: spatial distribution of the water depth , flow velocity and cumulative infiltration at t = 13,500 s.

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.

Table 3

Set of cases in Araguás catchment assuming a constant distribution

Test case α Kα (mα/s) ψ (m) Δθ (m3/m3
4.1  0.02 3.0 
4.2 0.95  0.02 3.0 
4.3 0.9  0.02 3.0 
4.4 0.9  0.02 3.0 
4.5 0.7  0.02 3.0 
4.6 0.9  0.035 3.0 
4.7 0.7  0.12 3.0 
4.8 0.9  0.02 6.0 
4.9 0.7  0.02 43.0 
Test case α Kα (mα/s) ψ (m) Δθ (m3/m3
4.1  0.02 3.0 
4.2 0.95  0.02 3.0 
4.3 0.9  0.02 3.0 
4.4 0.9  0.02 3.0 
4.5 0.7  0.02 3.0 
4.6 0.9  0.035 3.0 
4.7 0.7  0.12 3.0 
4.8 0.9  0.02 6.0 
4.9 0.7  0.02 43.0 
Figure 10

Araguás catchment, Cases 4.1 to 4.3.

Figure 10

Araguás catchment, Cases 4.1 to 4.3.

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 .

Figure 11

Araguás catchment, Cases 4.4 to 4.9.

Figure 11

Araguás catchment, Cases 4.4 to 4.9.

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.

Therefore, we propose to use a variable-order of the fractional derivative which is given in terms of the cumulative infiltration F and the water depth h. This will lead to an order of the derivative variable in time and distributed in space (horizontal plane) with a different value for each computational cell. The dependency on h is necessary for long and multiple rainfall events. For this reason, and taking into account that the infiltration rate is also controlled by surface water availability, the variable is formulated as a combined function of the cumulative infiltration and the available surface water as follows:  
formula
(24)
where a and b are constant values to calibrate and is the minimum value that the order of the derivative can reach. Then, it holds that .

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.

Table 4

Event 1: Infiltration parameters for the Araguás catchment

Inf. model Kα (mα/s)  ψ (m) Δθ (m3/m3a (m−1b (m−1   
GA  0.02 3.0 – – – 0.047 – 
FOGA  0.02 3.0 68 250 0.5 0.033 −29.3% 
Inf. model Kα (mα/s)  ψ (m) Δθ (m3/m3a (m−1b (m−1   
GA  0.02 3.0 – – – 0.047 – 
FOGA  0.02 3.0 68 250 0.5 0.033 −29.3% 
Figure 12

Graphical representations of the two terms of Equation (24): (a) with , , and (b) with , . (c) 3D representation of Equation (24).

Figure 12

Graphical representations of the two terms of Equation (24): (a) with , , and (b) with , . (c) 3D representation of Equation (24).

Figure 13

Hydrograph fitting for Event 1 .

Figure 13

Hydrograph fitting for Event 1 .

Figure 14

Water depth h, cumulative infiltration F and for Event 1 in Araguás catchment at t = 15,000 s.

Figure 14

Water depth h, cumulative infiltration F and for Event 1 in Araguás catchment at t = 15,000 s.

Arnás catchment

Catchment description and meshing

The Arnás catchment (2.84 km2, 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.

Figure 15

Arnás catchment characteristics. (a) Catchment location and hypsometry map. (b) Manning's roughness coefficient map. (c) Catchment refined mesh.

Figure 15

Arnás catchment characteristics. (a) Catchment location and hypsometry map. (b) Manning's roughness coefficient map. (c) Catchment refined mesh.

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)).

Events description

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.

Figure 16

Observed hyetographs and hydrographs for Arnás basin. (a) Storm Event 2: outlet discharge. (b) Storm Event 3: outlet discharge.

Figure 16

Observed hyetographs and hydrographs for Arnás basin. (a) Storm Event 2: outlet discharge. (b) Storm Event 3: outlet discharge.

Numerical results

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.

Table 5

Event 2: Infiltration parameters for the Arnás catchment

Inf. model Kα (mα/s)  ψ (m) Δθ (m3/m3a (m−1b (m−1   
GA  0.01 3.5 – – – 0.447 – 
FOGA  0.01 3.5 53 90 0.5 0.371 –16.7% 
Inf. model Kα (mα/s)  ψ (m) Δθ (m3/m3a (m−1b (m−1   
GA  0.01 3.5 – – – 0.447 – 
FOGA  0.01 3.5 53 90 0.5 0.371 –16.7% 
Table 6

Event 3: Infiltration parameters for the Arnás catchment

Inf. model Kα (mα/s)  ψ (m) Δθ (m3/m3a (m−1b (m−1   
GA  0.025 2.0 – – – 0.361 – 
FOGA  0.025 2.0 70 150 0.5 0.268 –25.8% 
Inf. model Kα (mα/s)  ψ (m) Δθ (m3/m3a (m−1b (m−1   
GA  0.025 2.0 – – – 0.361 – 
FOGA  0.025 2.0 70 150 0.5 0.268 –25.8% 
Figure 17

Hydrograph fitting for Event 2 .

Figure 17

Hydrograph fitting for Event 2 .

Figure 18

Hydrograph fitting for Event 3 .

Figure 18

Hydrograph fitting for Event 3 .

Figure 19

Water depth h, cumulative infiltration F and for Event 2 in Arnás catchment at t = 15,000 s.

Figure 19

Water depth h, cumulative infiltration F and for Event 2 in Arnás catchment at t = 15,000 s.

CONCLUSIONS

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.

ACKNOWLEDGEMENTS

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.

REFERENCES

REFERENCES
Benson
,
D.
,
Meerschaert
,
M.
&
Revielle
,
J.
2013
Fractional calculus in hydrologic modeling: a numerical perspective
.
Advances in Water Resources
51
,
479
497
.
35th Year Anniversary Issue
.
Borthwick
,
M.
2010
Application of Fractional Calculus to Rainfall-Streamflow Modelling
.
PhD Thesis
.
University of Plymouth
.
Caviedes-Voullième
,
D.
,
García-Navarro
,
P.
&
Murillo
,
J.
2012
Influence of mesh structure on 2D full shallow water equations and SCS Curve Number simulation of rainfall/runoff events
.
Journal of Hydrology
448–449
,
39
59
.
Costabile
,
P.
,
Costanzo
,
C.
&
Macchione
,
F.
2012
Comparative analysis of overland flow models using finite volume schemes
.
Journal of Hydroinformatics
14
(
1
),
122
135
.
Costabile
,
P.
,
Costanzo
,
C.
&
Macchione
,
F.
2013
A storm event watershed model for surface runoff based on 2D fully dynamic wave equations
.
Hydrological Processes
27
,
554
569
.
Deng
,
Z. Q.
,
Singh
,
P.
&
Bengtsson
,
L.
2004
Numerical solution of fractional advection-dispersion equation
.
Journal of Hydraulic Engineering
130
(
5
),
422
431
.
Deng
,
Z. Q.
,
de Lima
,
J.
,
de Lima
,
M.
&
Singh
,
V.
2006
A fractional dispersion model for overland solute transport
.
Water Resources Research
42
,
W03416
.
Diethelm
,
K.
2010
The Analysis of Fractional Differential Equations. An application-oriented exposition using differential operators of Caputo type. Lecture Notes in Mathematics. Volume 2004. Springer-Verlag, Berlin
.
Fernández-Pato
,
J.
,
Caviedes-Voullième
,
D.
&
García-Navarro
,
P.
2016
Rainfall/runoff simulation with 2D full shallow water equations: sensitivity analysis and calibration of infiltration parameters
.
Journal of Hydrology
536
,
496
513
.
Fiedler
,
F. R.
&
Ramírez
,
J. A.
2000
A numerical method for simulating discontinuous shallow flow over an infiltrating surface
.
International Journal for Numerical Methods in Fluids
32
,
219
240
.
García-Ruiz
,
J.
,
Lana-Renault
,
N.
,
Beguería
,
S.
,
Lasanta
,
T.
,
Regueés
,
D.
,
Nadal-Romero
,
E.
,
Serrano-Muela
,
P.
,
López-Moreno
,
J.
,
Alvera
,
B.
,
Martí-Bono
,
C.
&
Alatorre
,
L.
2010
From plot to regional scales: interactions of slope and catchment hydrological and geomorphic processes in the Spanish Pyrenees
.
Geomorphology
120
,
248
257
.
Gerasimov
,
D.
,
Kondratieva
,
V.
&
Sinkevich
,
O.
2010
An anomalous non-self-similar infiltration and fractional diffusion equation
.
Physica D: Nonlinear Phenomena
239
,
1593
1597
.
Green
,
W.
&
Ampt
,
G.
1911
Studies on soil physics: 1. Flow of air and water through soils
.
Journal of Agricultural Science
4
,
1
24
.
Horton
,
R.
1933
The role of infiltration in the hydrologic cycle
.
Transactions American Geophysical Union
14
,
446
460
.
Kilbas
,
A.
,
Srivastava
,
H.
&
Trujillo
,
J.
2006
Theory and Applications of Fractional Differential Equations
.
North-Holland Mathematics Studies. Volume 204
.
Elsevier Science B.V.
,
Amsterdam
.
Lana-Renault
,
N.
2007
Respuesta hidrológica y sedimentológica en una cuenca de montaña media afectada por cambios de cubierta vegetal: la cuenca experimental de Arnás, Pirineo Central (Hydrological and Sedimentological Response in an Average Mountain Basin Affected by Changes in Vegetation Cover: the Arnás Experimental Basin)
.
PhD Thesis
.
Universidad de Zaragoza
.
Langhans
,
G.
,
Nearing
,
G.
,
Diels
,
J.
,
Stone
,
J. J.
&
Diels
,
M. A.
2014
Modeling scale-dependent runoff generation in a small semi-arid watershed accounting for rainfall intensity and water depth
.
Advances in Water Resources
69
,
65
78
.
Liang
,
D.
,
Özgen
,
I.
,
Hinkelmann
,
R.
,
Xiao
,
Y.
&
Chen
,
J. M.
2015
Shallow water simulation of overland flows in idealised catchments
.
Environmental Earth Sciences
74
,
7307
7318
.
Machado
,
J.
,
Kiryakova
,
V.
&
Mainardi
,
F.
2011
Recent history of fractional calculus
.
Communications in Nonlinear Science and Numerical Simulation
16
,
1140
1153
.
Martinez
,
F.
,
Pachepsky
,
Y.
&
Rawls
,
W.
2010
Modelling solute transport in soil columns using advective-dispersive equations with fractional spatial derivatives
.
Advances in Engineering Software
41
,
4
8
.
Civil-Comp Special Issue
.
Mein
,
R.
&
Larson
,
C.
1973
Modeling infiltration during a steady rain
.
Water Resources Research
9
,
384
394
.
Murillo
,
J.
,
García-Navarro
,
P.
,
Burguete
,
J.
&
Brufau
,
R.
2007
The influence of source terms on stability, accuracy and conservation in two-dimensional shallow flow simulation using triangular finite volumes
.
International Journal for Numerical Methods in Fluids
54
,
543
590
.
Oldham
,
K.
&
Spanier
,
J.
1974
The Fractional Calculus. Theory and applications of differentiation and integration to arbitrary order
.
With an annotated chronological bibliography by Bertram Ross, Mathematics in Science and Engineering, Vol. 111. Academic Press (A subsidiary of Harcourt Brace Jovanovich, Publishers), New York-London
.
Pachepsky
,
Y.
,
Timlin
,
D.
&
Rawls
,
W.
2003
Generalized Richards’ equation to simulate water transport in unsaturated soils
.
Journal of Hydrology
272
,
3
13
.
Philip
,
J.
1969
Theory of Infiltration
. In:
Advances in Hydroscience
(
Chow
,
V. T.
, ed.).
Volume 5
.
Elsevier
,
New York
, pp.
215
296
.
Samko
,
S.
,
Kilbas
,
A.
&
Marichev
,
O.
1993
Fractional Integrals and Derivatives. Theory and applications
.
Gordon and Breach Science Publishers, Yverdon. Translated from the 1987 Russian original, revised by the authors
.
Serrano-Pacheco
,
A.
2009
Simulación numérica bidimensional de procesos hidrológicos e hidráulicos sobre lecho irregular deformable (Two-Dimensional Numerical Simulation of Hydrological and Hydraulic Processes on A Deformable Irregular bed)
.
PhD Thesis
.
Universidad de Zaragoza
.
Simons
,
S.
,
Busse
,
T.
,
Hou
,
J.
,
Özgen
,
I.
&
Hinkelmann
,
R.
2014
A model for overland flow and associated processes within the Hydroinformatics Modelling System
.
Journal of Hydroinformatics
16
(
2
),
375
391
.
Sun
,
H.
,
Zhang
,
Y.
,
Chen
,
W.
&
Reeves
,
D. M.
2014
Use of a variable-index fractional-derivative model to capture transient dispersion in heterogeneous media
.
Journal of Contaminant Hydrology
157
,
47
58
.
Te Chow
,
V.
,
Maidment
,
D.
&
Mays
,
L.
1988
Applied Hydrology
.
McGraw-Hill Civil Engineering Series, McGraw-Hill Higher Education
,
New York
.
Vreugdenhil
,
C.
1994
Numerical Methods for Shallow Water Flow
.
Kluwer Academic Publishers
,
Dordrecht
,
The Netherlands
.
Yu
,
C.
&
Duan
,
J.
2017
Simulation of surface runoff using hydrodynamic model
.
Journal of Hydrologic Engineering
22
(
6
),
04017006
.