## 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

*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 where

*g*is the acceleration due to gravity.

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

*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 where is the hydraulic head and is the saturated hydraulic conductivity.

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

*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: 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 and 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.

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

*z*(see for example Diethelm (2010)). The solution of problem (19) is given by and therefore it is now a linear combination of and .

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

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

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.

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

Test case . | R (mm/s)
. | K (m_{α}/s)
. ^{α} | ψ (m) . | Δθ (m^{3}/m^{3})
. |
---|---|---|---|---|

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) . | Δθ (m^{3}/m^{3})
. |
---|---|---|---|---|

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 |

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.

## 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 km^{2} (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)).

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

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

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.

Test case . | α
. | K (m_{α}/s)
. ^{α} | ψ (m) . | Δθ (m^{3}/m^{3})
. |
---|---|---|---|---|

4.1 | 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) . | Δθ (m^{3}/m^{3})
. |
---|---|---|---|---|

4.1 | 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 |

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.

*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: 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.

Inf. model . | K (m_{α}/s)
. ^{α} | ψ (m) . | Δθ (m^{3}/m^{3})
. | a (m^{−1})
. | b (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) . | Δθ (m^{3}/m^{3})
. | a (m^{−1})
. | b (m^{−1})
. | . | . | . |
---|---|---|---|---|---|---|---|---|

GA | 0.02 | 3.0 | – | – | – | 0.047 | – | |

FOGA | 0.02 | 3.0 | 68 | 250 | 0.5 | 0.033 | −29.3% |

### Arnás catchment

#### Catchment description and meshing

The Arnás catchment (2.84 km^{2}, 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)).

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

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

Inf. model . | K (m_{α}/s)
. ^{α} | ψ (m) . | Δθ (m^{3}/m^{3})
. | a (m^{−1})
. | b (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) . | Δθ (m^{3}/m^{3})
. | a (m^{−1})
. | b (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) . | Δθ (m^{3}/m^{3})
. | a (m^{−1})
. | b (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) . | Δθ (m^{3}/m^{3})
. | a (m^{−1})
. | b (m^{−1})
. | . | . | . |
---|---|---|---|---|---|---|---|---|

GA | 0.025 | 2.0 | – | – | – | 0.361 | – | |

FOGA | 0.025 | 2.0 | 70 | 150 | 0.5 | 0.268 | –25.8% |

## 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

*Application of Fractional Calculus to Rainfall-Streamflow Modelling*

*PhD Thesis*

*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

*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*

*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*