Preserving the dominant physical processes in a lumped hydrological model

Journal of Hydrology 399 (2011) 121–131

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Preserving the dominant physical processes in a lumped hydrological model M.L.V. Martina a,⇑, E. Todini a, Z. Liu b a b

Department of Earth and Geo-Environmental Sciences, University of Bologna, Piazza di Porta San Donato, 1, Bologna 40126, Italy Bureau of Hydrology, Ministry of Water Resources, 2 Lane 2, Baiguang Road, Beijing 100053, China

a r t i c l e

i n f o

Article history: Received 25 May 2009 Received in revised form 13 December 2010 Accepted 29 December 2010 Available online 13 January 2011 This manuscript was handled by K. Georgakakos, Editor-in-Chief Keywords: Lumped hydrological model Physically meaningful parameterisation Process and parameter lumping Geo-morphological lumping

s u m m a r y This paper analyses the relationship between the Saturated Area Ratio (SAR) and the Relative Water Content (RWC) within a hydrological catchment. The SAR–RWC relationship is fundamental in lumped modelling in correctly representing the dominant physical processes. According to the simulations provided by a distributed model (TOPKAPI) on several catchments, the SAR–RWC relationship is affected by hysteresis, which is also supported by observations and numerical simulations by other researchers. Due to the strong linearity, this relationship is difficult to derive analytically. In this paper we propose a physical interpretation of the hysteresis and present a methodology to model it. By incorporating this mechanism in a lumped model, we demonstrate improvement in discharge reproduction as evaluated by comparative streamflow simulations. Ó 2011 Elsevier B.V. All rights reserved.

1. Introduction The physically based models are the models that use physically based equations to describe the internal hydrological processes, which control the catchment response. They generally require a detailed discretization in time and space and the solution is obtained in ‘‘distributed’’ mode as opposed to the ‘‘lumped’’ mode, which is a characteristics of data driven and conceptual models (Todini, 2007). At present, physically based distributed models in contrast to the lumped models, which focus more on the input–output relationship, are recognized to be the most appropriate tools for simulating and investigating hydrological processes (Todini, 2007). Nonetheless the power of synthesis still makes lumped models potentially very attractive (OKane and Flynn, 2007; Sivapalan, 2003). The interest in lumped models lies not only in using simpler models but to a greater extent in identifying the ‘‘dominant processes’’, namely the essential hydrological features, to be preserved in the lumping process. Another reason for the current interest in lumped models is their use in Global Circulation Models (GCMs). Given the computational limits, GCMs use very simplified models to represent hydrological processes such as the surface runoff. These types of hydrological models often do not correctly preserve the non-linearity of the hydrological processes. In this case lumped models that correctly represent the physical phenomena may drastically improve the overall performance of GCMs. Another important aspect in this regard is its application on ungauged basins ⇑ Corresponding author. Tel.: +39 0512094535; fax: +39 0513370734. E-mail address: [email protected] (M.L.V. Martina). 0022-1694/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2010.12.039

where calibration is difficult if not impossible. Physically based models with measurable parameters, on the other hand, can be more easily applied. On one hand the historical evolution of hydrological models proceeded from simple conceptual models to more complex and physically based ones: detailed equations have been gradually introduced in an attempt to better reproduce the complex reality (Singh, 1998). On the other hand, various parsimonious lumped models have been developed with the aim of representing reality with widely different parameterisations of infiltration, soil saturation, drainage and runoff formation processes. However, the basic question still remains as to whether it is possible to directly set up a lumped hydrological model that encapsulates the physical properties and processes that operate at different scales, without having set up distributed models. The aim of the present paper is to show that the physical properties of the basic processes can only be retained at small spatial scales, while, due to inherent topological non-linearity, physically based lumped models need to be derived through an averaging process. This operation requires a correct representation of non-linear phenomena, such as soil water filling and depletion hysteresis and the exfiltration at the end of a rainfall event, which can be obtained via a distributed modelling approach. Starting from the point-scale equations and their physically meaningful parameterisation, the present work follows three lumping procedures: (1) a ‘‘process and parameter lumping’’ from the point equation to the pixel scale equation; (2) a ‘‘geo-morphological lumping’’ from the pixel equations to macrocatchment scale equations; (3) an ‘‘empirical lumping’’ of new dominant processes detectable by the distributed

122

M.L.V. Martina et al. / Journal of Hydrology 399 (2011) 121–131

model, but generally not described by the catchment scale equations. The process and parameter lumping makes it possible to define a correct formulation of point processes at a finite spatial dimension, and provides an indication on the limiting dimensions of a pixel that can still be considered as physically representative, as shown by Martina and Todini (2008). Geo-morphological lumping, as it will be shown, aims at aggregating all the grid cells (pixels) of the catchment in one single cell with equivalent hydrological properties. This requires a ‘‘non-trivial’’ aggregation of the distributed model parameters and processes, defined in each single cell, into a unique lumped parameter value. The empirical lumping, which makes it possible to preserve the description of internal dynamics at the lumped scale, exploits the diagnostic skills of the distributed model to infer internal relations representing the dominant processes that are now based on average quantities and is achieved by deriving a set of functions via simulation with distributed model. However the most interesting aspect of the methodology is the representation and interpretation of the relationship between the saturated area extent, referred to herein as Saturated Area Ratio (SAR) and the soil water content, also here called Relative Water Content (RWC), of the catchment. The SAR–RWC relationship is fundamental to lumped models to compute correctly the runoff. According to the experiments conducted, the SAR–RWC relationship is affected by an evident hysteresis which is also supported by observations and numerical simulations by other authors (Mishra and Seth, 1996; Niedzialek and Ogden, 2004; OKane and Flynn, 2007; Norbiato and Borga, 2008). The strong non-linearity of this relationship makes it difficult to be derived. Here, we offer a physical interpretation of the hysteresis and propose a methodology to model it. An example of its application to a real catchment will complement the theoretical description later in the paper.

The problem is that none of the existing lumped models is actually representing it: in general, when precipitation stops, surface flow stops as well in these models and calibrated parameters will compensate for this lack of physical representativeness. 3. An overview of the lumped TOPKAPI model The TOPKAPI model was developed in two different forms: the distributed and the lumped (Todini and Ciarapica, 2001; Ciarapica and Todini, 2002; Liu et al., 2005). There has been several applications of the distributed model (Bartholmes and Todini, 2005; Martina et al., 2006; Diomede et al., 2008a,b; Liu et al., 2008; Vischel et al., 2008a) which showed the capability of the model to correctly reproduce not only the discharge but also other internal variables such as the soil moisture (Vischel et al., 2008b). The lumped TOPKAPI model has also been applied on several catchments showing good hydrological capabilities (Liu and Todini, 2002). However some limitations are recognized on the derivation of the relationship linking the soil water content of the catchment to the extent of the saturated area. This work improves this important component in lumped modelling and proposes a physical interpretation of it, but in order to illustrate this aspect it is necessary to provide a background. In the following paragraphs we present a brief overview of the main important features of the lumped TOPKAPI model such as the integration of the governing equations and their aggregation to the catchment scale, the resulting non-linear reservoir equations representing the main model components. For a detailed description of the model and related methodologies please refer to the previously cited works (Ciarapica and Todini, 2002; Liu and Todini, 2002; Martina and Todini, 2008). 3.1. Integration of the governing equations

2. Existing lumped hydrological models and effects of lumping Lumped hydrological models, from the simple impulse response to the more complex conceptual models have been in use for decades for the design of hydraulic defences, water resources allocation and flood forecasting and warning (Todini, 2007). The more physically oriented of these lumped models aim at reproducing the non-linear water balance occurring at a finite scale in the soil. This water balance is then extended to the entire catchment either by applying the model for different areas (sub-catchments), or by using a sort of statistical distribution of water storages as in PDM, TOPMODEL or HBV model (Moore and Clarke, 1981; Beven and Kirkby, 1979; Bergstrom and Graham, 1998). Nonetheless, no model lumped at the catchment scale was conceived and formulated by up-scaling to the catchment scale the finite pixel equations used in the distributed version, via lumping or in other words via an integration process. Lumping the pixel equations produced a lumped analogue of the distributed ones. The comparison of the results produced by both the distributed and the lumped versions, revealed a marked discrepancy in terms of soil moisture deficit: the soil moisture depletion in the lumped version was significantly smaller than what could be explained by evapotranspiration and other losses. A closer analysis revealed that a peculiar phenomenon which was evident in the distributed model was not retained by the lumped model. The amount of water reaching the channel network through the soil at the end of a rainfall event was significantly smaller than the contribution of the upper edges of the saturated areas: the longitudinal development and the terrain slopes of the latter being much larger than the ones relevant to the terrains close to the river network. It is this imbalance that creates the phenomenon of exfiltration, which depletes the soil much more effectively.

In a previous work, Martina and Todini (2008) illustrated the conditions and limitations for the integration of the governing equations. These are usually partial differential equations (PDEs) at ‘‘point-scale’’, while models refer to finite scale (i.e. the grid scale). The integration of the governing equations from the point to the finite dimension of a grid cell was shown to generate, in a moderate scale range, models which preserve the physical meaning (although as averages) of processes as well as of model parameters and correctly reproduce phenomena such as the soil moisture balance and the consequent saturated areas dynamics which is probably the most important driving mechanism in the formation of surface runoff. Two important topics are related to the integration of governing equations: (1) the numeric method adopted for the integration and (2) the, unavoidable, parameters averaging over the grid cell (or figuratively ‘‘pixel’’). With regard to the first topic, the system of governing equations can be solved numerically by merely replacing the partial derivatives by finite differences on a discrete numerical grid. In other words, finite difference numerical schemes provide the solutions of the discretized (algebraic) version of the original system of partial differential equations. This approach despite its simplicity presents many disadvantages and a more successful approach to integration is to exploit the characteristic structure of the equations. This is the case, for instance, of the characteristic lines method, where the original equations are integrated along the characteristic lines and the resulting system is solved on the model space–time grid (Chow, 1959; Eagleson, 1970; Singh, 1996, 1998; Liu and Todini, 2005). An alternative method is, according to a finite-element approach, to integrate first the simplified equations over the grid cell of the model and then to discretize the ‘‘pixel-equations’’ obtained. For the second aspect, the reproduction of the vertical water profile based on the integration of the governing equations

123

M.L.V. Martina et al. / Journal of Hydrology 399 (2011) 121–131

seemed, to a certain degree, not necessary for the overall representation of the horizontal water flow and it is actually possible to lump the horizontal flow equations at the grid cell resolution, while still preserving the physical meaning of the model parameters. The experiments conducted by Martina and Todini (2008) showed that the model, structured according to the proposed integration method, produces acceptable results upon a grid scale of the order of the kilometre. Beyond that scale, greater divergences from the correct solution (i.e. considering the parameters spatial variability and using a very fine grid resolution) are found on the sub-surface flow, showing an evident dependency of the solution on the spatial distribution of parameters within the grid cell. However, within the scale range, a good phenomena representation is achieved keeping the physical meaning of the model parameters since they are simple average than the real ones. 3.2. Lumping the governing equations at the catchment scale This paragraph describes the methodology used, starting from the grid cell equations to obtain the (mass and momentum) conservative equations of the TOPKAPI lumped model on larger scales, based on the following assumptions:  The first assumption is that the overall behaviour of the basin does not depend on the actual position of individual slopes, but on the overall distribution of slopes. In the transition to the lumped structure of the model, single slopes are no longer described according to their actual position, but rather according to an overall distribution function.  The second assumption is that the time variation of the water content is assumed to be constant in space. This is the fundamental assumption on which the development of the lumped model is based and its validation is the starting point of the present study. In the lumped TOPKAPI model (Todini and Ciarapica, 2001; Liu and Todini, 2002), a catchment is regarded as a dynamic system which comprises three reservoirs: the soil reservoir, the surface reservoir and the channel reservoir. The precipitation over the catchment is partitioned into direct runoff and infiltration, on the basis of the soil condition and actual evapotranspiration by using an appropriate curve (described in the sequel) which reflects the non-linear relationship between the soil water storage and the saturated contributing area in the basin. Infiltration and direct runoff are then entered into the soil reservoir and the surface reservoir, respectively. The outflows from the two reservoirs, namely the interflow and the overland flow, are finally drained into the channel reservoir to form the channel flow (Fig. 1). This conceptual scheme is widely used for many lumped models, but the lumped version of TOPKAPI includes the following essential features: (a) The lumped equations for the three reservoirs are directly derived from the distributed model geomorphologic characteristics as described in previous works such as Todini and Ciarapica, 2001; Liu and Todini, 2002; Liu et al., 2005); (b) The estimation of non-linear functions which make it possible to (i) partition the precipitation into infiltration and surface runoff and (ii) to compute exfiltration directly from the results of the distributed model. The main aim in the development of the lumped model is in fact to preserve the physical representation of the phenomena without introducing new parameters and losing the physical meaning of the distributed version parameters. With regard to issue (a) – derivation of the lumped equations – a solution has been

Precipitation

[P] Infiltration

Direct Runoff

[φP]

[(1-φ)P]

Soil

[α s , β s ]

Subsurface flow [Qs]

Exfiltration [Qexf]

Surface [αo , β o ] Overland flow [Qo]

River [αc , βc ] Outflow Discharge

[Q] Fig. 1. Conceptual scheme of the TOPKAPI lumped model.

already proposed and successfully used (Todini and Ciarapica, 2001; Liu and Todini, 2002; Liu et al., 2005). Nonetheless a problem remains in terms of issue (b), that is, the estimation of the non-linear functions. In particular, in the referenced papers, the partitioning function between infiltration and exfiltration was based on an empiric law resulting from the simulations of the distributed model, which produced a saturated area vs. basin water storage content curve similar to those used in most of the variable saturated area models such as the PDM (Moore and Clarke, 1981), the Xinanjiang (Zhao, 1977), the ARNO (Todini, 1996, 2002a), the VIC (Liang et al., 1996a, Liang et al., 1996b). In the case of the lumped TOPKAPI model, this curve is a two parameters Beta function fitted on the results of the distributed model simulation. Nevertheless, this approach does not completely explain the physical dynamics at the soil level. The present work proposes a new partitioning function, along with its physical interpretation, by generalising the results obtained in the case of a simple plane hillslope. 3.3. Non-linear reservoir equations for the lumped model A detailed description of the non-linear reservoir formulation for the lumped model at the catchment scale is given in Liu and Todini, 2002; Todini and Ciarapica, 2001; Liu et al., 2005. The authors recall here only the final results of the derivation. The kinematic wave equation is integrated over all the grid cells composing the catchment in order to yield the lumped formula of the non-linear reservoirs for the watershed. This integration is performed under the assumption that during the transient (unsteady) phase the variation in time of the water content is not ignored as it is in TOPMODEL but assumed to be uniform in space over the pixel. This is the fundamental assumption used for the development of the TOPKAPI lumped model, which leads to an ordinary differential equation expressing the variation in time of the total water volume stored in the soil, on the surface and in the river network.

124

M.L.V. Martina et al. / Journal of Hydrology 399 (2011) 121–131

Under a set of simplifying assumptions, the resulting equations for all the three components composing the model still preserve the form of a non-linear reservoir equation as the distributed model does for each cell. Moreover, the parameters of the non-linear watershed reservoirs are directly derived from the parameters of the single non-linear cell reservoirs of the distributed model. The equations obtained respectively for the soil, surface and river components are:

dV sT ¼ u  P  Q exf  bs V csTs dt

ð1Þ

dV oT ¼ ð1  uÞ  P þ Q exf  bo V cooT dt

ð2Þ

dV cT ¼ Q s þ Q o  bc V cccT dt

ð3Þ

where the subscript T stands for the total volume of the water in the catchment, while the subscripts s, o and c, refer respectively to the soil component, overland component and channel component. The derivation and the meaning of the parameters b and c can be found in Annex 1 of Liu and Todini, 2002. They essentially depend on the topographic characteristics of the catchment, such as slope, and on hydraulic properties, of all the cells of the catchment, such as hydraulic conductivity and porosity, which are already used in the distributed model. P is the precipitation, u 6 1 is the precipitation fraction which infiltrates, Q exf is the exfiltration flow, Q s and Q o are the input discharge contributions to the channel component from the sub-surface flow and overland flow respectively. In order to solve the non-linear reservoir equations for the three components, it is necessary to compute at each time step both u and Q exf because they do not explicitly depend on the averaged soil water content of the catchment, but on its spatial distribution. This

(a) Isochrone lines of the water table (b) Isochrone lines of the water table profile in the wetting up transition phase profile in the wetting up transition phase staring from complete dry soil. staring from medium saturated soil.

(c) Isochrone lines of the water table profile in the drying down transition phase staring from complete saturated soil without Evapotranspiration.

(d) Isochrone lines of the water table profile in the drying down transition phase staring from complete saturated soil with Evapotranspiration.

(e) Water table profiles relative at the (f) Comparison of the water table profiles steady state of some precipitation event of for a drying down phase and for a steady different intensity state with the same saturated areas extension Fig. 2. Water table profiles produced by the TOPKAPI distributed model on synthetic plane hillslope. The upslope is on the left side and the distance in x-axis refers to this point.

M.L.V. Martina et al. / Journal of Hydrology 399 (2011) 121–131

is a problem of crucial importance for the accuracy of the lumped model. However, as it will be shown in the sequel, it is possible to consider the fraction of the saturated area, the Saturated Area Ratio, neglecting its actual spatial distribution. 4. The proposed methodology to improve the physical representation within the lumped model The following sections describe the methodology proposed to estimate the saturated area extent on the basis of the catchment soil water content necessary for the computation of u and Q exf values. 4.1. The estimation of the saturated area extent The governing equations, namely the water balance and the momentum equation, can be written in aggregated form (Eqs. (1)-(3)) and make it possible to directly compute the solution at the catchment scale provided that the input to each non-linear reservoir is known. Consequently, one has to add a mechanism for splitting the precipitation into the two non-linear reservoirs representing the soil (infiltration) and surface (direct runoff) components. This can be done by defining the percentage u of water that infiltrates into the soil, as a function of the saturated area As:

u¼1

As ¼ 1  SAR AT

ð4Þ

where AT is the total area of the catchment and SAR stands for Saturated Area Ratio. This assumption derives directly from the hypothesis of a saturation excess runoff mechanism (Dunne, 1978). This mechanism was shown to be the dominant hydrologic process at the catchment scale, particularly in humid climates, and has already been successfully used in many other models such as PDM (Moore and Clarke, 1981), Xinanjiang (Zhao, 1977, 1984, 1992), ARNO (Todini, 1996, 2002a), VIC (Liang et al., 1996a,b) and TOPMODEL (Beven and Kirkby, 1979). The first four models are very similar and use a parametric function, without any link to physical properties, to relate the SAR to the mean soil water content (Zhao,

125

1992; Todini, 1996). In the TOPMODEL (Beven and Binley, 1992) the authors attempted to provide a physically based methodology for the identification of the saturated area through the definition of the ‘‘topographic index’’. Nonetheless, the steady state hypothesis they used is rarely realistic (Seibert et al., 2003), and the model parameters, if calibrated at different space scales, can vary hugely and make a physical interpretation rather questionable (Franchini et al., 1996) even when the original formulation of TOPMODEL is modified to correct the resulting mass balance error (Sauliner and Datin, 2004). It is thus fundamental to estimate the saturated area extension under more physically realistic assumptions. In the TOPKAPI lumped model, the parameters of the aggregated formulation derive directly from the parameters of the TOPKAPI distributed model, which have been shown to preserve their physical meaning for cell sizes raging from 10 m to 1 km (Martina and Todini, 2008). This range is large enough for most applications of distributed hydrological models. This result is of crucial importance since it is a necessary condition for the preservation of the physically based properties. In order to analyse the saturated area dynamics it is useful to refer to the results obtained by several authors. Ogden and Watts (2000) studied the saturated area genesis into non-convergent hillslopes by means of VS2D, an unsteady two-dimensional model for the saturated and unsaturated flow (Healy, 1990; Lappala et al., 1993). They showed that the time to reach a steady state required by the system is generally longer than a single precipitation event. It is therefore incorrect to estimate the saturated area extension on the basis of the steady state assumption and to use one-to-one relationship with average variables (such as the mean soil water content). Niedzialek and Ogden (2004) by means of Sandbox, a physically based distributed model, while looking for a correlation between physical properties of the catchment and response time of the saturated area extension, highlighted a characteristic hysteresis in the relationship between saturated area and mean soil water content on a synthetic V-shaped watershed. The same results were found by Martina (2004) while simulating a real catchment of medium-size (4000 km2) with the distributed TOPKAPI model (Fig. 3).

Fig. 3. The Reno catchment with the different stations simulated in this work.

M.L.V. Martina et al. / Journal of Hydrology 399 (2011) 121–131

To summarize the referenced findings based on both numerical simulations and field observations, the relationship between the saturated area of the catchment and its average soil water content must be: (1) non-linear (as other conceptual models already correctly interpreted), (2) non-univocal (one-to-one), but (3) hysteretic, and (4) time variant. This is why it is so difficult to determine such a relation on physical grounds at the catchment scale. In the present work a possible approach for its assessment is proposed, using the diagnostic power of the distributed model. Before describing the approach proposed, we will discuss the physical interpretation of the hysteretic relation between the SAR and the average soil moisture content, and to propose a method for estimating the limiting curves from the cloud of dots representing the possible states of the catchment in terms of soil moisture and saturated area fraction. 4.2. SAR–RWC relationship on a hillslope A numerical experiment was set up on a plane hillslope 100 m long and 1 m wide with a 1 m soil layer over an impervious bed. The soil layer was divided into 100 cells of 1  1 m and characterised by a hydraulic conductivity of 105 m/s and a slope of 1.0%. The dynamics of the soil moisture content during and after precipitation events of differing intensity and duration were simulated by means of the distributed TOPKAPI model. The soil moisture content is converted into the saturation percentage a [0, 1] for each cell, which corresponds to a soil water content a L, where L is the soil thickness. What is interesting to observe for this case of simple plane hillslope is the relation between the saturated area extent and the water content profile. The latter can be analytically derived on the basis of the unsaturated flow equations (Martina, 2004), for three different conditions: (a) steady state, (b) wetting-up transient phase and (c) draining-down transient phase. Following Fig. 2a–f, it is interesting to note that: (1) In the wetting phase, when the soil fills under constant precipitation in time and space, the water content profile rises parallel to itself, no matter the initial state condition, while the profile reaches its equilibrium state from uphill to downhill. (2) The steady state profile shows the maximum saturated area for a given soil water storage content. (3) In the drying phase, e.g. after the end of a precipitation event, the saturated area decreases more rapidly than the corresponding soil water storage. Consequently, owing to (1), the relation between saturated area extension and soil water storage increases linearly during the wetting phase. This is not so for more complex slopes (concave or convex). However, it is theoretically possible to determine this law by analysing the topography. The consequence of (2) is that the locus of all the steady states is a limiting boundary for all the transient states (see Fig. 2e). The observation number (3) is seen to justify the unsteady assumptions of TOPKAPI as opposed to assuming that the water table instantaneously reaches the steady state, as in TOPMODEL. In the lumped TOPKAPI it was in fact assumed (Todini and Ciarapica, 2001; Liu and Todini, 2002) that the water table rises parallel to itself with a rate equal to the precipitation intensity, i.e.:

dg ffip dt

ð5Þ

where g is the water table depth or, equivalently in our definition, the relative saturation level for each cell. 4.3. SAR–RWC relationship on natural catchments The presence of the above-described complex saturated area dynamics and the corresponding water table profile is also seen in application to natural catchments. The distributed TOPKAPI model was implemented on the Reno catchment, approximately 4000 square kilometres in a hilly-mountainous area of mild and wet climate (see Fig. 3). The average soil Relative Water Content (RWC) and the Saturated Area Ratio (SAR) have been computed for different sub-basins of various sizes (see Figs. 4 and 5) using a 5-year time series of observed hourly precipitation available at 12 gauging stations. In order to verify the steady-state results described in Section 5.1, the RWC and the SAR have been computed at equilibrium, i.e. after the steady state was reached, for several synthetic precipitation events of differing intensity, assumed constant in space and time. These assumptions, while generally not necessary, make very simple the application of the methodology aimed here just at the estimation of the boundaries of the curves. Fig. 4 shows the results on the area–volume diagram and the Steady States (SS) limiting curve which envelops the cloud of dots, which were generated by the distributed model using the observed hydrological time series, as its lower boundary. The second boundary of the RWC–SAR cloud of dots, relevant to the drying phase, can be found by analysing the drying process of the catchment from the complete saturation. The drying process is primarily controlled by the evapotranspiration, which strongly depends on the solar radiation, on the temperature and on the soil moisture. In order to estimate the boundary curve, the authors decided to simulate a drying phase produced by a constant evapotranspiration rate equal to the maximum monthly value for the Reno catchment. As Figs. 4 and 5 show, these conditions seem to be the ones that better describe the upper boundary of the cloud of dots, referred to as the Drying Curve (DC) in the sequel. Following the results of the simulations, we propose the conceptual scheme of Fig. 6 for the saturated area dynamics. A

[Calcara] 1

DC

0.8

RWC

126

SS

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

SAR Fig. 4. Saturated Area Ratio (SAR) vs. Relative Water Content (RWC) for a 5 year simulation period of the Reno catchment at the Calcara river section (grey dots). The steady state curve (SS) was obtained as a set of SAR–RWC values for different simulations relative to precipitation events of different intensity and after the equilibrium state has been reached. Drying Curve (DC) represents the SAR and RWC values for the drying down transition phase only.

127

M.L.V. Martina et al. / Journal of Hydrology 399 (2011) 121–131

[Casalecchio]

[Calcara]

0.5 0 0

0.5

0.5 0 0

1

SAR [Castenaso]

0.5

0.5

0.5

1

SAR [Vergato]

0.5 0 0

0.5

1

SAR

1

RWC

1

1

1

0.5 0 0

1

0.5

SAR [PonteBastia] RWC

RWC

0.5

SAR [Mordano] RWC

0 0

1

1

0 0

0.5

SAR [SestoImolese]

1

RWC

1

RWC

1

RWC

RWC

1

[Forcelli]

0.5 0 0

0.5

0.5 0 0

1

SAR

0.5

1

SAR

Fig. 5. Same as previous figure for the other sub-basins of the Reno catchment.

1

se

ha

which has the minimum admissible SAR value, i.e., in Fig. 6, to the intersection between the DC curve with the straight horizontal line starting from the current state. Once this point has been reached, the draining process will follow along the DC curve until a new precipitation event occurs. This simplified dynamics would explain in fact the hysteresis described in the previous paragraphs. Based on the model simulations, the following important observation may be made.

Drying phase

p ng

i

RWC

y Dr

g Wettin

phase

Steady States curve Drying Boundary curve

0 0

1

SAR Fig. 6. Conceptual scheme of the saturated area extension dynamics for the wetting up and drying down phases on the SAR–RWC plan.

‘‘typical’’ cycle can be thought as a sequence of a storm – no storm – storm. A dot in the SAR and RWC area–volume domain represents an aggregate state of the catchment. Positioning it on the SAR–RWC domain, as soon as a new storm event begins, it is possible to compute the equilibrium state which the system tends to, given the observed precipitation. In fact this state is on the SS boundary curve (Fig. 6). In the conceptual scheme this phase has been sketched on the SAR–RWC plan by means of a straight line.. Note that if the precipitation rate amount changes, the equilibrium state to which the system tends will be different, but still computable. Until the steady state is reached, the system state will be represented by a dot over the ‘‘wetting phase’’ line. Once the steady state is reached, the system will stay on the same SAR–RWC state by definition. As soon as the precipitation ends, the catchment will almost instantaneously tend to a state with the same RWC value

(1) While the SS curve depends on the precipitation intensity amount, the DC curve depends on the evapotranspiration rate. (2) Although the horizontal line assumed for the drying phase (Fig. 6) is a simplification, nonetheless, it is a realistic interpretation of what happens. Such a result has been obtained in the synthetic case on the hillslope and in the numerical experiments using the distributed TOPKAPI model over the Reno catchment and it is also confirmed by Niedzialek e Ogden (2004). Similarly in their work, the drying phase was characterised first by a quick decrement of the SAR at constant volume, with practically no reduction of RWC, followed by a slower decrement of the SAR along a so called ‘‘envelope area’’. On the SAR–RWC domain this area coincides with the DC curve relevant to for the drying phase. 4.4. SAR–EXF relationship on a natural catchment In addition to what has been found for the SAR–RWC relationship, the authors have used the results from the distributed TOPKAPI model to analyse another important process, namely the exfiltration (also called return flow) (Liu and Todini, 2002). The term exfiltration is here used to indicate the amount of water returning to the surface from the soil caused by an unbalance between the incoming sub-surface flow amounts and the smaller ones that can flow downstream due to a limited draining capacity in the soil. Exfiltration contributes to the overland flow, and has an

128

M.L.V. Martina et al. / Journal of Hydrology 399 (2011) 121–131

linear reservoir of the soil component, in two different conditions. (1) In the wetting-up phase SAR is computed referring to the straight line on the SAR–RWC plan (Fig. 6), connecting the current state of the catchment (SAR0, RWC0) with that of the steady state relative to the current precipitation intensity. (2) For the time interval subsequent the end of a precipitation event, SAR is computed according to the DC curve (Fig. 6). Both SS and DC curve resulted from a simulation with the distributed TOPKAPI model. The return flow Q exf is computed, using the SAR–EXF relationships resulting from the simulation of a draining phase with a distributed model. Also note that Q exf is computed only during the dry periods since its contribution during a storm is negligible.

important role in depleting the soil moisture with the consequent modification of the initial soil moisture conditions in case of a successive rainfall event Indeed, the distributed TOPKAPI model under very wet catchment conditions clearly shows the presence of the exfiltration phenomenon, which is not usually reproduced in lumped model formulations, and in particular in the soil moisture accounting conceptual models. In order to include the exfiltration process in the lumped model, the distributed TOPKAPI model was again used to find a relationship with the SAR. Fig. 7 shows the discharge produced for exfiltration (EXF) during a drying phase starting from the catchment at full saturation. As can be seen from the diagram this relation is strongly non-linear and in a number of catchments is close to a threshold-like behaviour. However if (a) the SAR–EXF relation as simulated during the drying phase only is compared to (b) the SAR–EXF dots from the continued 5-years simulation, the curve in (a) seems to be just the lower boundary of the clouds of dots in (b). Based on a post analysis of the SAR–EXF relationship, the empirical function used in the model only represents the lower boundary of the cloud of SAR–EXF states showed by the catchment simulation.

6. Application on a real case The lumped TOPKAPI model has been implemented on the Reno catchment in the river section of Casalecchio in order to compare the simulated discharges with the observed ones as well as with those resulting from the distributed model. The adopted grid cell resolution of the distributed model was 1 km, and the time step was an hour. The precipitation was recorded by a gauging network of 12 stations (shown in Fig. 3). The calibration has been performed just on the distributed model, while for the application of the lumped model re-calibration was not needed but the derivation of aggregated parameters has been performed following the procedure described above and in the original paper by Liu and Todini (2002). It was also interesting to assess the benefits of the introduction of the SAR–RWC relationship in the new formulation which considers the hysteresis dynamic. As was previously discussed, Liu and Todini (2002) already provided an algorithm for model lumping but they used a simple parametric relationship to represent the SAR–RWC curve, thus neglecting the important physical features embedded in the dynamics. In order to show the contributions of the new SAR–RWC formulation, a simulation was also performed using the lumped model with a parametric function

5. Integration of the SAR–RWC dynamics in the TOPKAPI lumped model As the conceptual scheme of Fig. 6 clearly shows, the application of the TOPKAPI lumped model requires the solution of the three non-linear reservoir equations (Eqs. (1)-(3)) which depends on the determination of the partition factor u the return flow Q exf . In the previous paragraph was proposed a method to dynamically compute these factors which can now be implemented in the model. It is useful to summarize the findings. The partition factor u is assumed equal to 1 – SAR, i.e. direct runoff occurs only on saturated area. The estimation of SAR is based on the soil water content (RWC) which is explicitly computed by the water balance equation embedded in the non-

[Calcara]

1 0 0

0.5

−6

1 0

1

0

SAR

0

1

1 0.5 0 0

0.5

1

SAR

[Vergato]

2

−6

EXF (m/s 10 )

EXF (m/s 10−6)

0.5 SAR

[Mordano]

1

[PonteBastia]

0

1

0.5 SAR

0.5

SAR

1 0 SAR

0

EXF (m/s 10−6)

0

0.5

0

1

1

−6

1

0

1

[SestoImolese] EXF (m/s 10 )

EXF (m/s 10−6)

[Castenaso]

0.5

0.5

2

SAR

2

0

[Forcelli]

2

EXF (m/s 10 )

−6

EXF (m/s 10 )

EXF (m/s 10−6)

[Casalecchio] 2

1

2 1 0 0

0.5

1

SAR

Fig. 7. Runoff per unit area produced by exfilitration (EXF) vs. Saturated Area Ratio of the catchment (SAR). Grey dots are the data from the 5 year hydrologic simulation, black dotty line is the empirical function used by the TOPKAPI lumped model.

129

M.L.V. Martina et al. / Journal of Hydrology 399 (2011) 121–131

a

1200

observed distributed new lumped original lumped

800

3

Discharge [m /s]

1000

600

400

200

0 01/01/98

01/01/99

01/01/00

01/01/01

Date

b

1200

observed distributed new lumped original lumped

Discharge [m3/s]

1000

800

600

400

200

0 01/10/00

01/11/00

01/12/00

01/01/01

Date Fig. 8. Hydrographs for the Reno catchment at the Casalecchio river section as simulated by the distributed, the original (Liu and Todini, 2002) and new (this work) lumped TOPKAPI models compared with the observed. In (a) is shown the time series for the whole simulation period of 5 years, in (b) is shown only the event of November 2000.

1200

100 distributed new lumped original lumped

90

800

80 70

RMSE [m3/s]

Simulated discharge [m3/s]

1000

600

400

60 50 40 30 20

200

10

0

0

0

200

400

600

800

1000

1200

3

100

300

500

700

900

1100

3

Observed discharge ranges [m /s]

Observed discharge [m /s] Fig. 9. Simulated discharge by the new lumped model vs. observed discharge.

Fig. 10. RMSE relative to the different models at different ranges of observed discharge.

130

M.L.V. Martina et al. / Journal of Hydrology 399 (2011) 121–131

(Beta function) for SAR–RWC originally used by Liu and Todini (2002), instead of the complete formulation described in the present paper. Therefore, there are three models under comparison: (a) the distributed TOPKAPI, (b) the new lumped TOPKAPI (proposed by this work), (c) the original TOPKAPI (in the formulation originally proposed by Liu and Todini (2002)). As it can be seen from Fig. 8 the hydrograph is well reproduced both by the distributed and the lumped model but the difference in terms of computation time is huge. Particular attention should be paid to the quality of the peak discharge which is heavily by the estimation of the direct runoff as a function of the extent of the saturated area. Fig. 8 also shows that there are still some deficiencies, due to the exfiltration estimation in the quality of the recession limb. In Fig. 9 the scatter plot shows the discharge simulated by the lumped model vs. the observed discharge. The simulated discharge is very close to the observed both in the lower and in the higher values. Indeed Fig. 10 reports the RMSE (Root Mean Square

Errors) at different ranges of the observed discharge. It can be seen how the new lumped model clearly improves the results of the original and produces errors much closer to the distributed model. Besides the capability of reproducing the outlet discharge, it would be interesting to assess also the capability of the model to simulate the RWC and the SAR, but there are not observed data available for this purpose. For the only purpose of comparing the two formulations of the lumped model, we therefore took as reference the results from the distributed model which is not the ‘‘reality’’ but we think it represents well enough the complexity of the RWC–SAR dynamic. Figs. 11 and 12 shows respectively as the RWC and the SAR are better simulated both in the wet and dry periods by the new lumped model as compared to the lumped model originally proposed by Liu and Todini (2002). This is due to the new formulation of the RWC–SAR which preserves the important physical features of the dynamics such as the non-linearity and the hysteresis by means of cycles of wetting and drying phase.

1 0.9 0.8

RWC

0.7 0.6 0.5 0.4 0.3

Distributed New lumped Original lumped

0.2 0.1 01/01/98

01/01/99

01/01/00

01/01/01

Date Fig. 11. Simulated RWC by the models under comparison.

0.8

0.7

0.6

SAR

0.5

0.4

0.3

0.2

Distributed New lumped Original lumped

0.1

0 01/01/98

01/01/99

01/01/00

Date Fig. 12. Simulated SAR by the models under comparison.

01/01/01

M.L.V. Martina et al. / Journal of Hydrology 399 (2011) 121–131

7. Conclusions The aim of this paper was to show that the physical properties of the basic processes can only be retained at small spatial scales, while, due to inherent topological non-linearity, physically based lumped models need to be derived through an averaging process. This operation requires a correct representation of non-linear phenomena, such as soil water filling and depletion hysteresis and the exfiltration at the end of a rainfall event, which can be obtained via a distributed modelling approach. In order to preserve the essential features of the rainfall–runoff process from the point-scale to the grid cell scale or larger scales, mathematical formulations that integrate the point equations to the scale of the cell are essential to correctly reproduce the physical phenomena. Moreover a correct integration of the differential equations from the point to the finite dimension of a grid cell, and from the grid cell to larger scales, can actually generate relatively scale independent models, which preserve the physical meaning (although as averages) of the model parameters. However, the correct derivation of the governing equations and of the model parameters is not sufficient for a lumped model to preserve the physics of the phenomena simulated. An important dynamic to be captured is that of the saturated area (SAR) and its variability as a function of the Relative Water Content (RWC). It has been shown that the SAR–RWC relationship is: (1) non-linear (as other conceptual models already correctly interpreted), (2) non-univocal, (3) hysteretic, and in general, (4) time variant. The present work proposed, in the framework of the lumped model, the conceptualization of the SAR–RWC relationship by means of cycles of wetting and drying phase. The new lumped model can reproduce not only the outlet observed discharge but also the time variability of the SAR and RWC as simulated by the fully distributed physically based model. References Bartholmes, J., Todini, E., 2005. Coupling meteorological and hydrological models for flood forecasting. Hydrol. Earth Syst. Sci. 9, 333–346. Beven, K., Binley, A., 1992. The future of distributed models: model calibration and uncertainty prediction. Hydrol. Process. 6 (3), 279–298. Bergstrom, S., Graham, L.P., 1998. On the scale problem in hydrological modelling. J. Hydrol. 211, 253–265. Beven, K.J., Kirkby, M.J., 1979. A physically based, variable contributing area model of basin hydrology. Hydrol. Sci. 24, 1–3. Ciarapica, L., Todini, E., 2002. TOPKAPI: a model for the representation of the rainfall–runoff process at different scales. Hydrol. Process. 16 (2), 207–229. Chow, V.T., 1959. Open-channel Hydraulics. International Student Edition. McGraw Hill Book Company Inc. Diomede, T., Nerozzi, F., Paccagnella, T., Todini, E., 2008a. The use of meteorological analogues to account for LAM QPF uncertainty. Hydrol. Earth Syst. 12, 141–157. Diomede, T., Davolio, S., Marsigli, C., Miglietta, M.M., Moscatello, A., Papetti, P., Paccagnella, T., Buzzi, A., Malguzzi, P., 2008b. Discharge prediction based on multi-model precipitation forecasts. Meteorol. Atmos. Phys. 101, 245–265. Dunne, T., 1978. Field studies of hillslope flow processes. In: Kirkby, M.J. (Ed.), Hillslope Hydrology. Wiley, New York, pp. 227–293. Eagleson, P.S., 1970. Dynamic Hydrology. McGraw-Hill, New York. Franchini, M., Wendling, J., Obled, C., Todini, E., 1996. Physical interpretation and sensitivity analysis of the TOPMODEL. J. Hydrol. 175, 293–338. Healy, R.W., 1990. Simulation of Solute Transport in variably Saturated Porous Media with Supplemental Information on Modifications to US Geological

131

Survey’s Computer Program VS2D, US Geol. Surv. Water Resour. Invest. Rep., 904025. Lappala, E.G., Healy, R.W., Weeks, E.P., 1993. Documentation of Computer Program VS2D to Solve the Equations of Fluid Flow in variably Saturated Porous Media, US Geol. Surv. Water Resour. Invest. Rep., 83-4099. Liang, X., Wood, E.F., Lettenmaier, D.P., 1996. Surface soil moisture parameterization of the VIC-2L model: evaluation and modifications. Global Planet Change 13, 195–206. Liu, Z., Martina, M.L.V., Todini, E., 2005. Flood forecasting using a fully distributed model: application of the TOPKAPI model to the Upper Xixian catchment. Hydrol. Earth Sys. Sci. 9 (4), 347–364. Liu, Z., Todini, E., 2002. Towards a comprehensive physically-based rainfall–runoff model. Hydrol. Earth Syst. Sci. 6 (5), 859–881. Liu, Z., Todini, E., 2005. Assessing the TOPKAPI non-linear reservoir cascade approximation by means of a characteristic lines solution. Hydrol. Process. 19 (10). Liu, Z., Tan, B., Tan, X., Xie, A., 2008. Application of a distributed hydrologic model to flood forecasting in catchments of different conditions. J. Hydrol. Eng. 13 (5), 378–384. Martina, M.L.V., 2004. The Distributed Physically-based Modelling of Rainfall– Runoff Processes, Ph.D. Dissertation, The University of Bologna. Martina, M.L.V., Todini, E., 2008. Watershed hydrological modelling: toward physically meaningful processes representation, vol. 63. In: Hydrological Modelling and the Water Cycle, Springer, pp. 229–241. Martina, M.L.V., Todini, E., Libralon, A., 2006. A Bayesian decision approach to rainfall thresholds based flood warning. Hydrol. Earth Syst. Sci. 10, 413–426. Mishra, S.K., Seth, S.M., 1996. Use of hysteresis for defining the nature of flood wave propagation in natural channels. Hydrol. Sci. J. 41 (2), 15370. Moore, R.J., Clarke, R.T., 1981. A distribution function approach to rainfall–runoff modelling. Water Resour. Res. 17, 1367–1382. Niedzialek, J.M., Ogden, F.L., 2004. Numerical investigation of saturated source area behavior at the small catchment scale. Adv. Water Resour. 27, 925–936. Norbiato, D., Borga, M., 2008. Analysis of hysteretic behaviour of a hillslope-storage kinematic wave model for subsurface flow. Adv. Water Resour. 31 (1), 118–131. Ogden, F.L., Watts, B.A., 2000. Saturated area formation on nonconvergent hillslope topography with shallow soils: a numerical investigation. Water Resour. Res. 36 (7), 1795–1804. OKane, J.P., Flynn, D., 2007. Thresholds, switches and hysteresis from the pedon to the catchment scale: a non-linear system theory. Hydrol. Earth Syst. Sci. 11 (1), 44359. Sauliner, G.-M., Datin, R., 2004. Analytical solution to a bias in the TOPMODEL framework balance. Hydrol. Process. 18, 1195–1218. Seibert, J., Bishop, K., Rodhe, A., McDonnell, J., 2003. Groundwater dynamics along a hillslope: a test of the steady state hypothesis. Water Resour. Res. 39 (1), 1014. Singh, V.P., 1996. Kinematic Wave Modelling in Water Resources Surface-water Hydrology. John Wiley & Sons Inc., New York. Singh, V.P., 1998. Hydrologic Systems: Rainfall–Runoff Modelling, vol. 1–2. Prentice Hall, A Division of Simon & Schuster, Englewood Cliffs, New Jersey 07632. Sivapalan, M., 2003. Process complexity at hillslope scale, process simplicity at the watershed scale: is there a connection? Hydrol. Process. 17 (5), 1037–1041. Todini, E., 1996. The ARNO rainfall–runoff model. J. Hydrol. 175, 339–382. Todini, E., 2002. The ARNO model. In mathematical models of large watershed hydrology. In: Singh, V.P., Frevert, D.K., Meyer, S.P. (Eds.). Water Resources Publications, Littleton, Colorado, USA, pp. 687–716 (Chapter 16). Todini, E., 2007. Hydrological catchment modelling: past, present and future. Hydrol. Earth Syst. Sci. 11, 468–482. Todini, E., Ciarapica, L., 2001. The TOPKAPI model. Mathematical Models of Large Watershed Hydrology. In: Singh, V.P. et al. (Eds.). Water Resources Publications, Littleton, Colorado. Vischel, T., Pegram, G., Sinclair, S., Parak, M., 2008a. Implementation of the TOPKAPI model in South Africa: initial results from the Liebenbergsvlei catchment. Water SA 34, 1–12. Vischel, T., Pegram, G.G.S., Sinclair, S., Wagner, W., Bartsch, A., 2008b. Comparison of soil moisture fields estimated by catchment modelling and remote sensing: a case study in South Africa. Hydrol. Earth Syst. Sci. 12, 751–767. Zhao, R.J., 1977. Flood Forecasting Method for Humid Regions of China. East China College of Hydraulic Engineering, Nanjing, China. Zhao, R.J., 1984. Watershed Hydrological Modelling. Water Resources and Electric Power Press, Beijing. Zhao, R.J., 1992. The Xinanjiang model applied in China. J. Hydrol. 135, 371–381.