Description
Introduction
Continuum is a continuous distributed hydrological model, usually named HMC, that strongly relies on a morphological approach, based on a novel way for the drainage network components identification. The model is a compromise between models with a strong empirical connotation, which are easy to implement but far from reality, and complex physically based models which try to reproduce the hydrological processes with high detail but which introduce a hard parameterization and consequent uncertainty and lack of robust parameters identification.
The HMC aims at an equilibrium between simplicity and rigorous physical modeling while maintaining comparable performances to existing models; the reduced complexity of the schematizations and the relatively small number of parameters leads to a considerably lower calibration effort, increasing model robustness and portability to data scarce environments. The resulting increased computational efficiency leads to an easier utilization of the model in ensemble mode so that a direct quantification of prediction uncertainties is possible. All the main hydrological phenomena are modeled in a distributed way. Nowadays satellite products of meteorological, hydrological and “vegetation” variables are becoming promising and widely available all over the world. Examples of meteorological satellite products are: TRMM evapotranspiration estimations, H-SAF soil moisture. These kind of data are now available for a high percentage of the earth’s surface. This leads to a panorama where the possibility of running a hydrologic model using only satellite information is real. Along these lines develops one of the most important features of the HMC: it can be calibrated using only satellite data, e.g. surface temperature and/or soil moisture; results in reproducing observed discharge in this framework are quite promising. This makes the model suitable for application in data scarce environments.
The model is able to reproduce the spatial-temporal evolution of soil moisture, energy fluxes, surface soil temperature and evapotranspiration. Moreover, it can account for seasonal vegetation variability in terms of interception and evaporation. Deep flow and water table evolution are modelled with a simple scheme that reproduces the main physical characteristics of the process and a distributed interaction between water table and soil surface with a low level of parameterization. The introduction of the so-called “force-restore” equation for the surface energy balance allows the LST estimation and makes the model feasible to exploit remote sensing data. These latter can be used for calibration or, more appropriately, in a data assimilation framework. Referring to already tested calibration methodologies and making some basic assumptions, the calibration task can be reduced to just six parameters at catchment scale that are then spatially distributed by means of simple assumptions on the physical processes that they describe. Consequently, the parameter space is really small for a distributed continuous model, and HMC can be implemented with easily accessible data and territorial information (digital elevation model, basic soil and vegetation parameters). If more detailed territorial information were available, the parameterization methodology could be approached reducing the number of assumptions and linking the parameters more tightly to territorial characteristics. The sensitivity analysis has been carried out on five parameters: two parameters regulate the overland flow, the shape of the hydrograph and response time; two parameters are related to the soil characteristics and affect infiltration and relative soil humidity; and one parameter takes into account soil porosity and the ratio between vertical and horizontal saturated soil conductivity. A nice separation can be found between parameters in the fact that they influence quite distinct response of the model, with the result of further simplifying the calibration procedure. Standard calibration and validation based on streamflow data have been carried out on different periods and on different outlet sections (in terms of soil use, slope and response time). The model produces good results in terms of discharge. Further work is needed to introduce the modelling of the snow cover evolution and of the snowmelt in order to carry out multi-annual simulations where data are available. Model initialization influences the simulation for quite a long period; in particular the definition of initial water table level is a sensitive choice. The ideal situation is to extend as far as possible the warm-up period, including a long period without rainfall events.
The water table initialization methodology here is based on a possible spatial distribution of water table that exploits morphological constraints. Further validation analysis has been carried out comparing LST estimated by the model and satellite measurements at both pixel and basin scale. The results provide evidence that HMC reliably models the LST dynamics at various temporal scales, with some periods of overestimation, particularly during the warmer hours of summer. During the cold season the modelled LST has a lower variability with respect to the satellite estimates, but here the percentage of reliable data is quite scarce because of the more frequent cloud covering, and this makes the comparison more uncertain. The approach followed in the design of HMC proposes concentrating the efforts in augmenting the number of state variables that are predicted by the model and those that are also observables by using classical or remote instruments of measure. Specific attention is paid to distributed variables (e.g. LST fields) that offer very different information when compared to integral measures (e.g. discharge time series). The LST comparison showed potential for additional constraints to be used in the calibration phase or to be exploited in a more comprehensive assimilation framework. The distributed nature of the LST in comparison to traditional calibration time series (e.g. discharge data series) can add important information for a better estimation of state variables and parameter patterns. A demonstration of this potential is carried out by calibrating a sub-set of the parameters referring to LST satellite estimation and to morphologic information derived by the DEM. The results are comforting, and the proposed methodology led to a parameter set that well reproduces both satellite LST and streamflow observations, the latter used only in the validation phase. LST was successfully used to drive the calibration procedure. This was possible thanks to the model structure and the way parameters are treated and distributed in time and space. This could have a strong application impact in environments where reliable streamflow data are not available, given the worldwide availability of LST data.
Overland channel flow
The surface flow schematization distinguishes between channel and hillslope flow. In channels the momentum equation per unit of width is derived from the kinematic schematization with a nonlinear dependence between discharge and water velocity. The water depth for the i-th channel cell is evaluated combining the momentum equation and the mass balance equation:
where \(I_i\) represents the input per unit of area (the sumof runoff, saturation excess and inow discharge from upstream) to the grid cell [L T −1 ]. On the hillslopes the overland ow has a linear equation for the motion:
where \(u_h\) parameterizes the main morphologic characteristics of the hillslopes [T −1 ] (slope, roughness, etc.). The nal schematization is equivalent to a linear reservoir model. The parameters \(u_h\) and \(u_c\) need calibration at basin scale (i.e. one value for the entire catchment). In both hillslopes and channels, the re-infiltration process is accounted for: the input to the i-th cell must exceed its infiltration capacity; otherwise, it infiltrates the soil. Exfiltration is also possible.
Vegetation interception
Interception includes the portion of rainfall that is caught by tree leaves, grass and vegetation cover in general, and is evaporated before it touches the ground.
Ponding effects are also included in this initial abstraction. Interception is modelled by a simple empirical equation. A maximum retention capacity \(S_v\) is introduced, and it is estimated as a function of the leaf area index (LAI) by the relationship:
The water in the reservoir with capacity \(S_v\) is evaporated at the evaporation rate derived by the latent heat flux estimation without affecting the infiltration computation; the input is the precipitation. The advantage of using a LAI-dependent expression is that the model takes into account vegetation variability in space and time. LAI is usually updated every 15 days from satellite optical sensor data (e.g. from MODIS).
Infiltration and Subsurface Flow
The infiltration methodology is a modification of the Horton equation based on physically interpretable parameters. It accounts for soil moisture evolution even in condition of intermittent and low-intensity rainfall (namely lower than the infiltration capacity of the soil). The soil is schematized as a reservoir with capacity \(V_{max}\) [L], and a selective filter \(g(t)\) [L T−1 ] manages the inflow:
where \(f_0\) is the maximum infiltration rate for completely dry soils and \(f_1\) is the asymptotic minimum infiltration rate for saturated soils considered as a function of \(f_0\):
The method proposed has been further modified by introducing the field capacity of the soil, defined as the water content that can be held by capillarity against the force of gravity:
with the parameter \(c_t\) ∈ [0, 1]. In this configuration, the dynamic mass-balance equation for the soil can be written for each cell:
where:
The infiltration scheme has four parameters: the initial infiltration rate \(f_0\), the maximum soil retention capacity \(V_{max}\), and the parameters to define soil field capacity \(c_t\) and final infiltration rate \(c_f\).
The parameters \(f_0\) and \(V_{max}\) are related to the soil type and land use through the curve number (CN) parameter. They can be easily derived by soil use and soil type maps and they vary spatially in the catchment; \(c_t\) and \(c_f\) are calibration parameters and are assumed to be constant for the whole basin. In this way the pattern of \(f_1\) and \(V_{fc}\) is spatially modulated by the pattern of \(V_{max}\). The percolation rate separates into two components: a contribution to subsurface flow \(r_{Hy}\) and one to deep flow \(r_d\) or recharging water table defined as:
where the angle \(\alpha\) is such that \(\tan(\alpha)\) is the downslope index; \(sin(\alpha)\) is a decomposition term that increases with the terrain slope, reproducing the major proneness of the high slope areas to subsurface flow due to gravity. The downslope index seems to be less sensitive than local surface slope to changes in DEM resolution, and it is able to capture dominant controls on local drainage regimes, especially in cases where profile curvature exerts a strong control on the drainage pattern. The subsurface flow is propagated between cells following the surface drainage network directions, and the soil moisture state of each cell is updated by considering both the infiltration, estimated by the modified Horton method, and the inflow from the upper cells. Therefore a cell can reach saturation because of the percolation from upper cells causing saturation excess. The water in the reservoir with capacity \(S_v\) is evaporated at the evaporation rate derived by the latent heat flux estimation without affecting the infiltration computation; the input is the precipitation. The advantage of using a LAI-dependent expression is that the model takes into account vegetation variability in space and time.
Deep Flow and Water Table
Several approaches are possible to describe the dynamics of both the deep flow and the water table, with examples from Darcy’s law applications to conceptual reservoir models. However, it is often difficult to have the data necessary for the correct implementation and parameterization of water table dynamics. In HMC the water table evolution is modelled with a simplified approach that maintains a physical and distributed description of the process.
Above all we are interested in the water table interaction with the subsurface flow and soil surface and in its effects on surface flow and soil moisture spatial pattern; the adopted scheme allows also the reproduction of the baseflow far from rainfall events with a parsimonious parameterization. The layer of soil containing the aquifer is schematized as a unique homogeneous layer bounded by the lower impervious (bedrock) surface and the bottom of the root zone. The thickness of this layer is expressed in terms of maximum volume of water content of the aquifer, and it is estimated using the surface slope as a proxy. The maximum water content in every cell (i) of the basin is given by:
where \(VW_{max}\) is the absolute maximum water content of the aquifer on the whole investigated area; this sets a limit that is basically a calibration parameter. The reservoir is fed with rd (see previous section). The effect of porosity is considered as a multiplicative factor in the Darcy’s equation used to estimate the flux per unit area between two contiguous cells (i and j):
where x is the DEM spatial resolution, f1i the final infiltration rate estimated, \(h_{w}\) the water table level and \(R_{f}\) a factor that also takes care of differentiating the saturated vertical and horizontal conductivity. Each cell can drain towards all the neighbouring cells following the 2-D water table gradient that depends on the elevation and on the water content of each cell. When the water table reaches the surface (\(VW_{i} = VW_{mi}\)), the deep percolation term is inhibited, while the condition \(VW(t) ∼ 0\) is a limit that can only be reached after a very long and anomalous dry period.
Energy Balance and Evapotranspiration
The representation of surface mass and energy turbulent fluxes requires the solution of a conservation equation for mass and energy driven by temperature and moisture content. Since the vertical gradient of such variables is quite large, a high-resolution multiple layer model would be required to estimate soil surface temperature and moisture content with accuracy. Such an approach demands substantial amounts of computing resources to solve the balance equations.
An alternative approach makes use of computationally efficient parameterization of soil heat and moisture flux terms. In other studies is showed that the heat flux into the soil could be parameterized by the sum of a temperature-derivative term and the difference between ground surface and deep soil temperature. This approach is known as the “force-restore” method, because the forcing by net radiation is modified by a restoring term that contains the deep soil temperature. Since then the “force-restore” method has been widely used in land surface modelling. Some studies demonstrated that the “force-restore” equation is the solution of the heat diffusion equation, with purely sinusoidal forcing assuming that the thermal properties are constant with depth and the surface forcing term is also nearly independent of air temperature and has a strong periodic behaviour in time. The HMC solves a complete and explicit energy balance at the interface between soil surface and atmosphere by using the “force-restore” approach for land surface temperature. Theoretically, the control volume to which the balance is applied is the unit area bounded vertically by the surface of the soil and the top of the canopy, assuming the thermal capacity of this volume is negligible. The horizontal energy fluxes are neglected. In practice, the volume is extended to the unit cell of the numerical scheme used. This approximation is a fair trade-off between parsimony in parameterization and accuracy in the description of the processes.The conservation of energy at soil surface is given by:
where \(R_n\) is the net radiation, \(H\) the sensible heat flux, \(LE\) the latent heat flux and \(G\) the ground flux (all [E t−1 L−2 ]). This latter term closes the budget, and it represents the heat propagated by diffusion towards the deep layers of the soil. The shortwave component of \(R_n\) is derived from radiometer observations when the density of observations is appropriate. Otherwise, it is estimated by combining the extraterrestrial component of the radiation computed, attenuated using meteorological variables and cloud cover. The terrain parameter characterizations that influence both direct and diffuse components of the radiation are computed. The longwave components are rarely available from observations, and they are therefore estimated using the Stefan–Boltzmann law as a function of air temperature and humidity. The daily cycle of LST has the implicit signature of the energy balance. Maximum amplitudes of LST diurnal cycle are usually reached in the presence of bare and dry soil. The presence of moisture on the surface and in the subsurface soil greatly moderates the daily range of LST. The vegetation cover has a similar effect. The “force-restore” approach leads to the following equation for LST:
where \(\varphi\) φ [E L−2 T−1 t−(1/2) ] is the effective thermal inertia and \(T_{deep}\) [T] is a “restoring” deep ground temperature. \(T_{deep}\) is evaluated by filtering data for air temperature at ground level; \(\varphi\) is the thermal inertia, and it is a function of conductivity, density and specific heat capacity of soil, and it is eventually related to soil moisture. The fluxes are estimated using bulk formulations. The equation input variables are commonly observed by ground-based micrometeorological networks. The soil parameters used in the estimation of the thermal inertia, usually constant at basin scale, can be estimated by a data assimilation process, or related to soil type when reliable maps are available. In HMC the evapotranspiration \(LE\) [m s−1] is estimated as:
where \(\rho_{w}\) [m L−3 ] is the water density, and \(ET\) is deducted from the interception storage \(S_v\) if not empty, otherwise from the subsurface reservoir \(V(t)\) adding the following:
if \(S_v\) greater then 0 and:
if \(S_v\) greater equal 0.
Parameters
Six model parameters need calibration on the basis of input–output time series: \(c_{f}\), \(c_{t}\), \(h_{h}\), \(u_{c}\), \(R_{f}\), \(VW_{max}\). The first two parameters \(c_{f}\) and \(c_{t}\) mainly rule the generation of runoff and the movement of water in the different soil layers, while \(u_{h}\) and \(u_{c}\) control the surface water motion. \(VW_{max}\) represents the maximum storage capacity of the aquifer, and \(R_{f}\) summarizes the effect of soil porosity as well as of the ratio between vertical and horizontal saturated soil conductivity. The range of variation of the parameters has been defined based on prior knowledge of the parameter meaning, which defines their mathematical and physical range of validity.
Parameters |
Units |
Description |
uh |
[s-1] |
Flow motion coefficient in hillslopes |
uc |
[m-0.5s-1] |
Friction coefficient in channels |
cf |
[-] |
Defines the infiltration capacity at saturation |
ct |
[-] |
Defines the mean field capacity |
Rf |
[-] |
Related to anisotropy between the vertical and horizontal saturated conductivity, and to porosity |
VWmax |
[mm] |
Maximum water capacity of the aquifer on the whole investigated area |
The parameters \(u_{c}\) and \(u_{h}\) impact the water flow on the surface. High values of these two parameters lead to narrow and highly peaked hydrographs; \(u_{h}\) has influence on the general shape of the hydrograph while \(u_{c}\) has an increasing influence with the increasing length of the channelled paths (e.g. large/elongated basins). It modifies the peak flow value as well as the peak arrival time. The impact estimation of parameters \(u_{h}\) and \(u_{c}\) has been made considering a short period of simulation (16 to 18 August 2006) since they influence directly the overland and channel flow. The first subplot in the following figure shows that \(u_{h}\) has a considerable influence on both the Tiglieto and Casalcermelli outlets.
The peak values and the hydrograph shape have quite a large range of variation, while peak times are not significantly affected by this parameter. The second subplot shows the influence of \(u_{c}\). It mainly affects the shape and the peak times on the Casalcermelli outlet section, while hydrographs of the Tiglieto outlet show negligible differences. Note that Casalcermelli has a drainage area that is one order of magnitude larger in respect to Tiglieto. The parameter \(c_{t}\) is related to the soil field capacity and defines the fraction of water volume in the soil not available for percolation and subsurface flow. It has an impact on the dynamics of soil saturation between rain events: higher values of \(c_{t}\) reduce the soil drying time scale especially during the cold season, with consequently higher runoff coefficients for single rainfall events. However, the subsurface flow tends to vanish rapidly, because water level drops easily under the field capacity. The parameter \(c_{f}\) controls both the velocity of subsurface flow and the dynamics of saturation of the single cells. Low values of \(c_{f}\) (i.e. low values saturated hydraulic conductivity) tend to cause the rapid saturation during rainfall events associated with slow subsurface flow increasing runoff production. Higher values of \(c_{f}\) produce a rapid subsurface flow with saturated areas that quickly concentrate along the drainage network. The combination of the two soil parameters \(c_{t}\) and \(c_{f}\) controls the distribution of the volumes of soil and surface water in space and time, and it impacts soil humidity propagation; \(c_{t}\) and \(c_{f}\) influence the mass balance over long periods and regulate the exchanges between subsurface flow and runoff. The third and fourth panels in the previuous figure show how they affect the tails of the hydrographs and the values of peak flows in the period between 7 and 10 December 2006. The effect of the combination of these two parameters is quite complex, and it is only partially represented in the figures. They must be calibrated over long periods of time using, at best, external soil information when available.
The parameter \(R_{f}\) regulates the response of the deep flow and mainly influences low flow regimes, while for larger basins it also affects high discharges. In the figure below the period between 14 and 17 September 2006 is shown. The effects of \(R_{f}\) on the Tiglieto outlet are negligible during the flood while the influence on low flows is more relevant for both the outlet sections. Particular remarks need to be made about \(VW_{max}\) – a measure of the capacity of the basin for storing water in its aquifer and deep soil layer. All these simulations highlighted another important feature of the model.
Because of its internal structure, similarly to other complete hydrological models, it is possible to map different processes and therefore different parts of the hydrograph, onto the parameters, so that different parts of the hydrograph time series can be used separately to better identify model parameter values. Further analysis is needed to show sensitivity to spatial and temporal resolution.
It is not easy to define a value for \(VW_{max}\) that reproduces a correct or realistic distribution of the deep soil layer water storage throughout the basin due to the fact that this distribution is hard to observe. Tests made using different values of \(VW_{max}\) in a physically acceptable range and starting from the same initial condition show that the model has low sensitivity to this parameter when the period of simulation covers between 6–12 months. This is related to the slow temporal dynamic of the water table.
If data series for very long simulations (many years) are available, the parameter \(VW_{max}\) can be re-calibrated and adjusted. In the adopted scheme the initialization of the related state variable \(VW(t)\) is more important than its upper limit. In fact, practice demonstrates that the definition of the water table initial condition \(VW(t = 0)\) evidently influences simulated discharge. A reasonable initial condition produces a rapid stabilization of the water table with dynamics driven by the water input from upper soil layer. Two considerations are made in order to define these:
in correspondence with the drainage network, the water table is generally next to soil surface, because it is continuously recharged by the upstream portions of catchment;
the mountainous parts of the water table receive reduced contribution, because they drain small areas and are characterized by high gradients, and here the water table has lower levels.
Based on these considerations, water table initial conditions are set as follows; \(VW(t = 0)\), in correspondence with channels, is set close to \(VW_{max}\). In the hillslopes the level of \(VW(t = 0)\) is estimated supposing it is inversely proportional to the downslope index \(\alpha\).
In order to carry out a basic sensitivity analysis, we considered what appear to be the most sensitive parameters (\(u_c\), \(u_h\), \(c_t\), \(c_f\)) and a set of 2000 model runs has been generated using a Monte Carlo approach, sampling the parameters from a uniform distribution bounded by the parameter domain.
The runs have been carried out on a sub-period of the calibration period where two major flood events occurred (July – September 2006). To show the results a dot plot representation has been drawn using the Nash–Sutcliffe coefficient (NS) as skill estimator. From figure above, it is possible to clearly identify a behavioural domain for \(u_h\) and \(u_c\), while a lot of uncertainty in the streamflow simulation due to \(c_f\) and \(c_t\) is present. This could be due to the very nonlinear relationship that connects these last two parameters with streamflow.
