Exploiting remote sensing land surface temperature in distributed hydrological modelling: the example of the Continuum model

Full process description and distributed hydrological models are very useful tools in hydrology as they can be applied in different contexts and for a wide range of aims such as flood and drought forecasting, water management, and prediction of impact on the hydrologic cycle due to natural and human-induced changes. Since they must mimic a variety of physical processes, they can be very complex and with a high degree of parameterization. This complexity can be increased by necessity of augmenting the number of observable state variables in order to improve model validation or to allow data assimilation. In this work a model, aiming at balancing the need to reproduce the physical processes with the practical goal of avoiding over-parameterization, is presented. The model is designed to be implemented in different contexts with a special focus on data-scarce environments, e.g. with no streamflow data. All the main hydrological phenomena are modelled in a distributed way. Mass and energy balance are solved explicitly. Land surface temperature (LST), which is particularly suited to being extensively observed and assimilated, is an explicit state variable. A performance evaluation, based on both traditional and satellite derived data, is presented with a specific reference to the application in an Italian catchment. The model has been firstly calibrated and validated following a standard approach based on streamflow data. The capability of the model in reproducing both the streamflow measurements and the land surface temperature from satellites has been investigated. The model has been then calibrated using satellite data and geomorphologic characteristics of the basin in order to test its application on a basin where standard hydrologic observations (e.g. streamflow data) are not available. The results have been compared with those obtained by the standard calibration strategy based on streamflow data.


Introduction
Continuous streamflow simulation is of fundamental importance in the support of water management decisions (e.g.best use of water resources) and civil protection actions (e.g.flood and drought mitigation actions) (Schlosser et al., 1997;Middelkoop et al., 2001;Karsten et al., 2002;Bartholomes and Todini, 2005).These are only some examples of the huge variety of cases where continuous hydrological models have been applied.The application of models to different problems resulted in the development of a number of hydrological models, which sometimes showed very different characteristics (e.g.Beven, 1997;Todini and Ciarapica, 2001;Rigon et al., 2006).
Over recent decades computation capacity has developed exponentially.Meanwhile, due to the progress of Earth observation techniques, a large amount of territorial information (digital elevation models, land use, soil and vegetation parameters) has become readily available.As result of this, full process description and distributed hydrological modelling, assisted by detailed catchment descriptions, has become feasible, leading to the improvement in the understanding and representation of both runoff formation and propagation dynamics (Winchell et al., 1998;Giannoni et al., 2000).Distributed modelling allows us to understand the role played by space and time rainfall distribution (Giannoni et al., 2003), by soil and vegetation heterogeneity, and by the Published by Copernicus Publications on behalf of the European Geosciences Union.
Several hydrological models are based on a detailed description of all the main hydrological processes so that they can provide a continuous simulation in support of specific tasks such as long-term water balances and flood forecasting.Some examples starting from the early 1970s can be cited: Sacramento model (Burnash et al., 1973;Burnash, 1995), SHE (Abbott et al., 1986), RIBS (Garrote and Brass, 1995), CASC2D (Julien et al., 1995), TOPMODEL (Beven, 1997;Wang et al., 2006), TOPKAPI (Todini and Ciarapica, 2001), GEOtop (Rigon et al., 2006), and MOBIDIC (Campo et al., 2006).They can usually work in different regimes and simulate continuously the spatiotemporal evolution of the state of the catchment.The advantages of these types of models have been highlighted in various works (e.g.Bartholomes and Todini, 2005;Liu et al., 2005;Castelli et al., 2009).
The use of models capable of exploiting distributed information has consolidated, with the penalty of increasing the number of model parameters with the associated problem of reliably estimating them (Abbott et al., 1986;Beven, 1993;Madsen, 2000;Anderson et al., 2006).
Because of the detailed and complete description of the hydrological cycle, these advanced models are characterized by an abundant number of parameters, which need a detailed knowledge of the catchment to which they are applied.Some of them, such as those describing soil physics and land use, have a direct physical meaning.They can therefore be estimated on the basis of maps of catchment information and geophysical cartography exploiting also the continuous improvement of satellite-derived information; nevertheless, the scale of representation, the detail and the temporal updating of information often do not match requirements for the reliable estimation of these parameters.Furthermore, relationships between land information and model parameters are sometimes indirect and affected by great uncertainty, nullifying, in practical applications, the benefit of an accurate process description.
The presence of many parameters, in general, allows modellers to obtain good results in the calibration phase because of the adaptation skills of the model; however, this leads to an increased probability of obtaining similar results with different parameter sets (Beven and Binley, 1992;Beven and Freer, 2001;Savenije, 2001).This limits the possibility of reliable parameter estimation, hampering the prediction of the abilities of the model.There are two ways to face this issue: develop a parsimonious model or increase the means of identifying parameters from available data.The design of the "Continuum" model is based on a combination of the abovementioned strategies.On the one hand, special attention is paid to reducing, as much as possible, the parameterization of the physical processes so that land information can be extensively used as a constraint to parameter calibration.On the other hand, the model parsimony has to be compatible with the necessity for a detailed description of all the terms of the hydrological cycle.This results in an increased number of observables in the model state variables, in order to exploit the potential offered by today's remote sensing (R. S.).
In the last decade many advances in satellite data development and analysis for Earth observation have been reached.R. S. provides inputs to the model (e.g.meteorological input, catchment description, vegetation characterization) or observation of the state variables (e.g.soil moisture, land surface temperature).It is therefore a concern for new models to include prognostic equations of satellite observables so that they can be used in the calibration/validation phase, or as constraints in a data assimilation framework (Winsemius et al., 2006;Kumar et al., 2008).
The Continuum model aims at an equilibrium between simplicity and rigorous physical modelling while maintaining comparable performances to existing models; this is demonstrated by a calibration and validation exercise using standard methods and data.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; along these lines this work proposes a novel calibration methodology based on remotely sensed data only.
The article is organized as follows.In chapter two the components of the model are presented; chapter three describes the model's sensitivity analysis, the case study, as well as the parameter calibration and model validation; the chapter includes the calibration experiment carried out using only R. S.-derived information (i.e.LST -land surface temperature -data and morphological information); chapter four contains discussion and conclusions.

Model description
Continuum is a continuous, distributed hydrological model that strongly relies on a morphological approach, based on a novel way for the drainage network component identification (Giannoni et al., 2005).Continuum is designed to be applied on a wide range of catchment classes, from small (O(Area) = 10 1 : 10 2 km 2 ) to medium (O(Area) = 10 3 km 2 ) size basins and with time resolutions lower than 24 h.Such scales are normally modelled using data-intense models (e.g.DHI, 2003), and application of such models to basins with O(Area) > 10 3 km 2 is possible (e.g.Bartholomes and Todini, 2005), but the data demand in order to have a reliable implementation is usually very high.
The basin is represented using a regular square mesh based on digital elevation model (DEM); the flow directions are identified using a D-8 approach (O'Callaghan and Mark, 1984).The drainage network is represented distinguishing between hillslope and channelled flow with a filter defined by the expression A S k = C where A is the contributing area upslope of each cell [L 2 ] and S the modified local slope [−] (Giannoni et al., 2000).The exponent k [−] is a function of the morphologic environment in which the network developed and weighs the importance between flow accumulation and slope in determining active channelled flow (Giannoni et al., 2000).In mature morphological environments, k can be set equal to 1.7 (Giannoni et al., 2003) while the threshold C can be calibrated following fractal theory (Giannoni et al., 2005) and reproducing topographic independent information like the so-called "blue lines".
Infiltration and subsurface flow is described using a semiempirical, but quite detailed, methodology based on a modification of Horton algorithm (Bauer, 1974;Disikin and Nazimov, 1997) and focuses especially on exploiting land use and climatology information, which are easily available, to set the infiltration parameters.The energy balance is based on the so-called "force-restore equation" (Dickinson, 1988), which balances forcing and restoring terms, with explicit soil surface temperature prognostic computation.The overland runoff is distributed with differentiation between hillslope and channel flow.Vegetation interception and water table flow have been also schematized.The different approaches are detailed in the following paragraphs.

Overland and 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 (Wooding, 1965;Todini and Ciarapica, 2001) with a nonlinear dependence between discharge and water velocity (see Appendix A for details on the equations).
The water depth for the i-th channel cell is evaluated combining the momentum equation and the mass balance equation (Liu et al., 2005): where I i represents the input per unit of area (the sum of runoff, saturation excess and inflow discharge from upstream) to the grid cell [L T −1 ].
On the hillslopes the overland flow has a linear equation for the motion: where u h parameterizes the main morphologic characteristics of the hillslopes [T −1 ] (slope, roughness, etc.).The final 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 similar to the one used by Rey (1999), Kozak et al. (2007), Zhao (2003), among others.A maximum retention capacity S v is introduced, and it is estimated as a function of the leaf area index (LAI) by the relationship (Kozak et al., 2007): The water in the reservoir with capacity S v is evaporated at the evaporation rate derived by the latent heat flux estimation (see Sect. 2.5) without affecting the infiltration computation; the input is the precipitation (see Fig. 1).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 (Diskin and Nazimov, 1994;Gabellani et al., 2008) 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).

F. Silvestro et al.: Exploiting remote sensing LST in distributed hydrological modelling
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 (Mishra and Singh, 2003): The method proposed by Gabellani et al. (2008) 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 new configuration (see Fig. 1), the dynamic mass-balance equation for the soil can be written for each cell: where The details of dynamic mass balance equations are described in Appendix B.
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 (Risse et al., 1995).Following Gabellani et al. (2008) 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 α is such that tg(α) is the downslope index as in Hjerdt et al. (2004); sin(α) 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 (see Hjerdt et al., 2004 for details on this issue).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 (Dunne and Black, 1970).

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 (Todini and Ciarapica, 2001;Rigon et al., 2006;Campo et al., 2006).However, it is often difficult to have the data necessary for the correct implementation and parameterization of water table dynamics (Castelli et al., 2009).
In Continuum 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 by following Saulnier et al. (1997) using the surface slope as a proxy.The maximum water content in every cell (i) of the basin is given by where V W 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 (see Fig. 2).The reservoir is fed with r d (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, f 1i the final infiltration rate estimated as in Eq. ( 5), h W the water table level (for details on the equation defining h W , see Appendix C) 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 (V Wi (t) = V W mi ), the deep percolation term in Eq. ( 9) is inhibited, while the condition V W (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 (Deardorff, 1978) 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.Bhumralkar (1975) and Blackadar (1976) independently 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.Deardorff (1978) referred to this approach 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 (e.g.Lin, 1980;Dickinson, 1988;Dickinson et al., 1993;Noihan and Planton, 1989;Caparrini et al., 2004).Hu and Islam (1995) 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 Continuum model 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 (Dickinson, 1988).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 (Boni et al., 2001;Sini et al., 2008).
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 on the basis of Yang et al. (2001Yang et al. ( , 2006)), attenuated using meteorological variables and cloud cover (including MSG cloud cover).
The terrain parameter characterizations that influence both direct and diffuse components of the radiation are computed following the approach in Dozier and Frew (1990).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  is evaluated by filtering data for air temperature at ground level (Caparrini et al., 2003;Caparrini and Castelli, 2004).ϕ 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.Details about the component of the energy balance and its parameterization are reported in Appendix D. 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 (Peters et al., 1997) when reliable maps are available.
In Continuum the evapotranspiration ET [m s −1 ] is estimated as where ρ 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 terms to Eq. (B1) in Appendix:

Model summary
In summary, Continuum is a distributed model based on a space-filling representation of the network, directly derived from a DEM.The DEM resolution coincides with the model resolution.The mass and energy balances are solved at cell scale referring to the schematizations of subsurface flow, deep flow and vegetation interception.The overland and channel flow are described by a linear and a nonlinear tank schematization respectively.Figure 3 offers a sketch of the Continuum structure.Two consecutive cells are considered to highlight the interaction between the single modules.
Continuum can work with different time steps ( t m ), generally smaller than 1 day.Time steps, in the range between 10 min to 1 h, are appropriate when the objective is a flood simulation in small to medium catchments, and they are consistent with commonly available meteorological data.The overland and channel flow module works with a time step ( t s ) that is smaller than t m , dependent on the DEM resolution to avoid problems of numerical stability.The mass-balance equations are integrated using the Heun time stepping scheme (Clark and Kavetski, 2010), which is a semi-implicit method (predictor-corrector scheme), and it is second-order accurate.Such a choice allows a good balance between accuracy of the solution and performance when an appropriate time step is used (this would change in applications where the integration step is large, e.g.daily).A different scheme is used to solve the "force-restore" equation where more accuracy is needed and a Runge-Kutta forthorder method is used.Since small time steps are used, minor effects due to the adopted time stepping scheme can be expected (Clark and Kavetski, 2010).
Six model parameters need calibration on the basis of input-output time series: c f , c t , u h , u c , R f , V Wmax (see Table 1).The first two parameters c f , and c t mainly rule the generation of runoff and the movement of water in the different Table 1.Summary of the model parameters that need calibration with their brief description.
Parameter Description Defines the infiltration capacity at saturation Defines the mean field capacity Related to anisotropy between the vertical and horizontal saturated conductivity, and to soil porosity Maximum water capacity of the aquifer on the whole investigated area soil layers, while u h and u c control the surface water motion.V Wmax 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.

Watershed description
The model has been tested on the Orba basin located in the southern part of the Piemonte region in the Apennines (see Fig. 4).The Orba River originates from Mt. Reixa (1183 m a.s.l.) in the Beigua Massif, and it flows into the Bormida River, a tributary of the Tanaro River, before it reaches the town of Alessandria.Three main morphological areas can be identified: a mountain part characterized by very steep sub-catchments with a very deep river bed, a mild part with an average slope of 1 %, and finally the alluvial part characterized by very small slope values.The Orba River has mainly a torrential regime with recurrent flash floods during the autumn and spring rainfall seasons and very low flows during summer.The Orba mean annual flow in correspondence of the confluence with Bormida River is around 20 m 3 s −1 (here the drainage area is about 810 km 2 ).

Dataset
The micrometeorological networks of Liguria and Piemonte Italian regions provide meteorological inputs.In the Orba basin there are 31 rain gauges, 27 thermometers, 6 hygrometers, 4 radiometers (shortwave) and 4 anemometers.The temporal resolution of the observations is 1 h.This latter is also the temporal resolution used in the model application ( t m = 1 h).The overland and channel flow module uses t s = 30 s.
Observations from two nested stage discharge gauges, with reliable and constantly updated discharge rating curves, are used for model calibration and validation (Fig. 4).The Tiglieto station (75 km 2 ) is located in the upper part of the basin, characterized by a mountainous morphology with high

Influence of the calibration parameters on simulations
A first set of 50 simulations has been carried out to evaluate the effects on the simulated hydrographs of the calibration parameters shown in Table 1.The most significant are used in the following in order to show how the parameters affect the simulations.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 (Beven and Binley, 1992;Maidment, 1992;Dingman, 2002).The ranges are reported in Table 2.
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 of Fig. 5 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 of Fig. 5 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 Fig. 5 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 Fig. 6 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 V Wmax -a measure of the capacity of the basin for storing water in its aquifer and deep soil layer.
It is not easy to define a value for V Wmax 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 V Wmax 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 V Wmax can be re-calibrated and adjusted.
In the adopted scheme the initialization of the related state variable V W (t) is more important than its upper limit.In fact, practice demonstrates that the definition of the water table initial condition V W (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: (i) 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; (ii) 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.V W (t = 0), in correspondence with channels, is set close to V Wmax .In the hillslopes the level of V W (t = 0) is estimated supposing it is inversely proportional to the downslope index α (Hjerdt et al., 2004).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 (Beven and Binley, 1992;Shen et al., 2012) has been drawn using the Nash-Sutcliffe coefficient (NS) (Nash and Sutcliffe, 1970) as skill estimator (Fig. 7).From Fig. 7 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.
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 (Reusser et al., 2009) 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.

Standard calibration based on streamflow data
The model has been firstly calibrated with a standard approach based on streamflow data.The dataset collected for the period 1 June to 31 December 2006 is particularly suitable for calibration under different hydrological conditions as it contains a range of significant floods and a preceding drought period.Runs of the model have been carried out for the whole year 2006, but the first five months are considered as model startup to ensure that the initial conditions do not affect the results.
The first part of the calibration period is quite dry with an absence of precipitation for more than two months.A very intense and persistent event started on 16 August 2006 and lasted about 10 h with peak rainfall intensities larger than 70 mm h −1 .On 14 September 2006 another intense event with a 24 h duration and considerable rainfall accumulation took place.
For the calibration period we referred to empirical formulations that use wind speed and air temperature data to estimate the bulk transfer coefficient for heat and moisture (Deardoff, 1968).They do not account for vegetation variability.
Several skill estimators are considered in order to evaluate the performance of Continuum: Nash and Sutcliffe (1970), Chiew and McMahon (1994), correlation coefficient (CORR), root-mean-square error (RMSE) and percentage error on six flood events in the period.The first three estimators assume the value 1 if observation and simulation match 100 %, while the other estimators tend to be 0 for a perfect simulation.
V WMax has not been considered for calibration for the reasons already explained, and it has been set to the value V WMax = 1500 mm.
The calibration has been performed in two steps.Initially, to reduce their range of variation, the two overland flow parameters have been calibrated on a short time period (16 to 18 August 2006).The errors on peak flow and on peak time on both the outlet sections have been minimized.Finally, a global calibration has been made including the three soil parameters c t , c f and R f maximizing the Nash-Sutcliffe coefficient, which is one of the most commonly used skill estimators in hydrology to compare observed and simulated hydrographs (Legates and McCabe, 1999).This has been carried out by varying parameter values in their range of validity following an iterative brute-force approach.Two iterations have been done; the minimum gradient varies for each parameter and it is reported in Table 2. Table 3 reports the set of parameters obtained by the calibration procedure.
Table 4 shows the skill estimators calculated for the entire simulation; Table 5 reports the percentage errors on peak flow (PPE) for the two outlet sections, Tiglieto and Casalcermelli.Figures 8 and 9 report the simulated hydrographs compared with observations; we plot the entire simulation     The bulk transfer coefficient for heat C H used in the energy balance equations is in general estimable by using empirical expressions that depend on wind speed, pressure and on vegetation characteristics that define the neutral transfer coefficient C HN .The vegetation parameterization is often derived by literature values, and it is considered constant in space and time because of the difficulty of evaluation in a real context.By using proper models that describe the vegetation and soil interaction with the atmosphere, it is possible to produce a more detailed C HN estimation that takes into account the vegetation spatial distribution and its seasonal variability.In this work the C HN estimation is carried out using a variational assimilation scheme as described in Caparrini and  Castelli (2004); maps of C HN are produced every 15 days for the entire validation period.The application of this model has been possible for the validation period only.

Streamflow validation
Table 6 reports the skill estimators for the whole validation period, while Table 7 reports the percentage errors of the peaks of the main flood events.Skill estimators assume similar or better values during validation period with respect to the calibration period; this indicates the robustness of the model.Figures 10 and 11 report the simulated hydrographs compared with observations for the two outlet sections.As can be seen, the model reproduces well the observed discharges.In particular, during the most important flood events, Continuum simulates the shape of the hydrograph, the peak flows and the peak times surprisingly well.
In both calibration and validation periods, in few periods the differences between observed and simulated streamflow are not negligible; this can be due to the flash flood regime of the catchment, where the precipitation distribution could seriously hamper the performance of the model in reproducing the hydrograph due to input volume mismatch: this problem was highlighted in previous works facing the same challenging goal (Liu and Todini, 2002;Rigon et al., 2006, Ruiz-Villanueva et al., 2012).

Land surface temperature analysis
The LSTs estimated by Continuum have been compared with LSTs retrieved by satellite measurements.The database of LST estimations provided by LAND SAT application facility (SAF) of EUMETSAT has been used (EUMET-SAT, 2009).Land surface temperature estimations are available every 15 min with a spatial resolution of about 0.04 • .Data are available for the validation period only.The retrieval of LST (Freitas et al., 2010) is based on clear-sky measurements from the MSG system in the thermal infrared window (MSG/SEVIRI channels IR108 and IR120).The analysis has been carried out for the period 1 June 2009 to 31 December 2009.
Due to the complex topography of the Orba basin, the LST satellite estimates cannot be directly compared to model outputs because of the following problems: (i) the georeferencing of model pixels and satellite pixels, (ii) the shadowing due to the presence of the mountains, (iii) the variation of the satellite viewing angle for the different detected areas, and (iv) the different spatial resolution (satellite-estimated spatial resolution is about 4.5 km while model output is about 0.1 km).A land application model (F.Castelli, personal communication, 2009) that projects the radiance obtained from the model (obtained from the simulated LST) to the same geometry of the satellite observations solves these four cited issues.The land application model produces a correlation matrix that weights the model radiance to estimate the portion of energy of each model pixel that contributes to the energy of the satellite pixel.The application of the land surface model can be formalized as where M is a matrix operator that weights the model output and maps it on the same grid of satellite data, ε m and ε o the model and the satellite thermal emissivity, LST m and LST o the modelled and the satellite-estimated land surface temperatures.The model assumes a constant ε m , and ε o is estimated  as the mean thermal emissivity of the two sensor channels used for LST-SAF retrieval.
We compared the basin mean LST and the LST of three selected pixels.The three pixels have been chosen from different areas of the basin to investigate the model behaviour in different environments.One pixel is in the mountainous part of the basin (we name it Mountain Pixel), one in the middle part (we name it Hill Pixel) and one near the outlet section in the flood plain (we name it Flat Pixel).Four skill estimators have been used to evaluate the performances of the model: the mean bias (BIAS), the root-meansquare error (RMSE), the mean absolute error (MAE) and the correlation coefficient (CORR).
Figures 12 to 15 report the scatter plot of satellite estimates and simulated LST for the mean basin comparison and for the pixel analysis.
The model slightly overestimates the LST during some periods of the warm season (Figs.16 and 17), while it slightly underestimates LST during the cold season, in particular during the central hours of the day.This behaviour could be related to a slight overestimation of the model soil thermal inertia.
However, the general trends and diurnal cycle amplitude are well reproduced and the skill estimators show quite a good performance (see Table 8); BIAS, RMSE, and MAE are quite small, while CORR is next to 1. Model errors are in most of the cases smaller than the root-mean-square error of satellite LST estimation with respect to in situ measurements, which is about 2-3 • C (Freitas et al., 2010).The model reproduces well the trend and the daily periodicity of LST (Figs. 16 and 17).
Figure 18 shows satellite estimates and simulated LST plotted with the mean incoming solar radiation at basin scale for a period of a few days.The LST is strongly related to the radiation, and the model is able to correctly reproduce the rapid changes in LST due to solar radiation variations.
The soil humidity is a factor that influences the thermal inertia and, as a consequence, the LST.We noted that, for example, the simulated LST overestimates the satellite estimates during a dry period at the beginning of summer (June-July), where Continuum produces very low values of soil humidity at basin scale.This can be related to the fact that AE and CORR are similar for all the target areas (Table 8), while BIAS is higher and positive on mountain pixel where Continuum tends to dry the soil quickly because of the percolation and the soil humidity propagation.This behaviour becomes evident in Fig. 19 that depicts the RMSE map of LST distributed in the catchment layered over the DEM relief.RMSE high values tend to concentrate in mountainous areas.Another interpretation of this result is that in such areas the variance of LST is higher due to altitude variation, and this influences the statistics when aggregated at the MSG pixel size.However, we can deduce that a better LST simulation could be obtained by varying model parameters so as to obtain a different soil humidity distribution.Further analysis  is needed to verify this hypothesis, and it was not carried out in this study.It is, however, interesting to note that LST could be actually a useful constraint in the calibration phase as discussed more in detail in the next section.

Parameter calibration using land surface temperature and morphological data
To demonstrate the possibility of reliably implementing the model in data-scarce environments, we performed a calibration supposing no streamflow data are available for the catchment under study.While it is always possible to find usable surrogates of meteorological variables using satellites or meteorological model analyses on one side and information about morphology and land cover is globally available at the needed resolution, on the other, it is often the case of not having reliable streamflow data in order to perform a standard hydrologic calibration of the model.In this exercise we used information that can be extracted by land cover maps, DEMs and the LST satellite estimations.The analysis is focused on the most sensitive parameters (u c , u h , c t , c f ); the other two parameters have been set equal to the values V Wmax = 2000 mm, R f = 1 assuming that we have no information in support of their identification.V Wmax is set to an average value for the environment considered, while R f is set to one making the hypothesis of isotropy of the soil characteristics.

Surface parameters (u h and u c )
The estimation of the overland and channel flow parameters has been carried out by exploiting geomorphologic information derived by the DEM.The idea is to calibrate the parameters against some hydrological time scales that can be derived by geomorphological information using wellestablished morphological relationships.One of the most influential parameters in that sense is the lag time (t l ) of the basin, defined as the temporal distance between the centre of mass of the hydrograph and the centre of mass of the mean hyetograph at basin scale.
The lag time (t l ) can be derived by geomorphologic features of the basin; we chose the formulation proposed by the National Resources Conservation Service (Ward and Trimble, 2003): where L is the hydraulic length of the watershed, i f the mean slope of the basin, S a function of the curve number representing an index of the soil water storage capacity.
We generated a set of synthetic rainfall events (Eq.7) with constant intensity in space and time to feed the model.In this way we generate hydrographs that mainly depend on the geomorphologic characteristics of the basin.The features of the synthetic events are reported in Table 9.They are characterized by increasing accumulated rainfall while their duration is a function of the basin area.We calculated the simulated lag time of these synthetic events for each combination of the parameters and minimized the objective function: where t lo is the t l derived by the geomorphologic characteristics of the basin and t lm is the t l obtained by the model simulations.
The soil has been considered impermeable in order to eliminate the nonlinearities related to the wetness condition and to the infiltration-runoff process.The absence of infiltration and of dependency on humidity conditions at the beginning of the event (in the t lo formula S is set to 0) makes the model parameters that regulate subsurface and deep flow  (c t , c f , R f and V Wmax ) irrelevant.L and i f are derived by the DEM.
On the basis of the analysis carried out in Sect.3.3, it is clear that, in upstream sections, with reduced paths in channelled network, the influence of the parameter u c is scarce; therefore a first guess value for u c can be fixed and u h calibrated.Once u h is calibrated, u c can be calibrated considering a downstream section.This procedure is iterated until  the two parameters reach stability.In the case study we referred to the two gauged stations of Orba basin: Tiglieto (A = 75 km 2 ) and Casalcermelli (A = 800 km 2 ).
For each P cum the procedure rapidly converged (N iterations = 3).
The optimal values of u c or u h are different for different values of P cum .This depends on the intrinsic characteristics of the model, and it is coherent with reality.The average of the optimal values and the standard deviation are calculated and reported in Table 10.
The mean values of the parameters have been used to proceed with the calibration.

Subsurface parameters (c t and c f ) -exploiting LST satellite estimate
Once the surface flow parameters are estimated, it is necessary to calibrate the subsurface soil parameters.To do this the LST retrieved by satellite observations can be compared with that modelled by Continuum with the objective of minimizing the differences.LST is influenced by the soil characteristics and by the humidity content; in Continuum the two parameters c t and c f regulate the subsurface water propagations as well as the field capacity in each cell and, as a consequence, the humidity patterns through the basin.
The mean LST at basin scale (LST Av ) has been considered, and the following score based on the BIAS between simulated (LST Avm ) and satellite (LST Avo ) estimates has been considered: where T max is the number of the considered temporal steps.The entire validation period has been used.The score should guarantee that, on long periods, the temporal average of model LST is similar to that estimated by satellite data.
The BIAS AV minimization returns the following values for the parameters: c t = 0.56 [−] and c f = 0.03 [−] with BIAS AV = 0.01.
The value of c t is weakly increased in respect to streamflow-based calibration as well as c f .This means that Table 11.Values of skill estimators for the calibration and the validation periods obtained using c t and c f calibrated basing on the LST and u c and u h calibrated using geomorphologic data.to reduce the total bias of LST, larger values of field capacity and infiltration capacity at saturation are needed.

Streamflow validation
The final set of parameters has been used for simulating both calibration and validation periods, and the results are reported in  4 and 6; this proves that such calibration approach produces a reliable implementation of the model that can be replicated in data-scarce environments where real problems in finding reliable input-output series can be encountered.In Figs.20 and 21 the hydrographs obtained through the two calibration approaches are reported for the Casalcermelli outlet section.The flood events are really similar, while there are more evident differences in the baseflow and in the recession periods of the hydrographs.

Discussion and conclusions
The article describes a distributed and continuous hydrological model that balances the necessity for a complete description of hydrological cycle with a simple and versatile structure resulting in a limited number of parameters.It has been designed for a variety of purposes: flood forecast and simulation, water resource management and droughts.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 (Giannoni et al., 2005;Gabellani et al., 2008) 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 Continuum 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 with reference to two different outlet sections of the Orba basin.These sections represent different characteristics in terms of soil use, slope and response time.The model produces good results in terms of discharge for both outlet sections.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 (Hjerdt et al., 2004).
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 Continuum 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 Continuum 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.

Fig. 1 .
Fig. 1.Sketches of vegetation retention and subsurface flow at cell scale.S v is the capacity of the vegetation reservoir, V max the capacity of the soil reservoir, V the actual water volume in the soil and c t V max the field capacity of the soil.

Fig. 3 .
Fig. 3. Sketch of flux partition in the Continuum model with the integration of the single modules.Two consecutive cells are illustrated.

Fig. 4 .
Fig. 4. Orba river location in north-west Italy.In the lower left corner is a zoom on the basin with the micro-meteorological stations.White squares represent the level gauges with rating curve, black squares the level gauges without rating curve.

Fig. 5 .
Fig. 5. Influence of the calibration parameters on hydrographs.c t and c f are the subsurface flow parameters while u h and u c are the overland and channel flow equation parameters.

Fig. 6 .
Fig. 6.Influence of the parameter R f on hydrographs; both the outlet sections are considered.The same period is shown with different y-axis scales to highlight the differences for both high and low flows.

Fig. 7 .
Fig. 7. Dot plot representation of the sensitivity analysis for the parameters u c , u h , c t , and c f .

Fig. 8 .
Fig. 8.Comparison between observed and simulated hydrographs for the calibration period for Casalcermelli outlet section.

Fig. 9 .
Fig. 9. Comparison between observed and simulated hydrographs for the calibration period for Tiglieto outlet section.

Fig. 10 .
Fig. 10.Comparison between observed and simulated hydrographs for the validation period for Casalcermelli outlet section.

Fig. 11 .
Fig. 11.Comparison between observed and simulated hydrographs for the validation period for Tiglieto outlet section.

Fig. 12 .
Fig. 12.Comparison between satellite estimates and simulated LST.Average LST at basin scale.The graph refers to the period from 1 June to 31 December 2009.

Fig. 13 .
Fig. 13.Comparison between satellite estimates and simulated LST, mountain pixel.Graph refers to the period from 1 June to 31 December 2009.

Fig. 14 .
Fig. 14.Comparison between satellite estimates and simulated LST, hill pixel.Graph refers to the period from 1 June to 31 December 2009.

Fig. 15 .
Fig. 15.Comparison between satellite estimates and simulated LST, flat pixel.Graph refers to the period from 1 June to 31 December 2009.

Fig. 16 .
Fig. 16.Comparison between air temperature, satellite estimates and simulated LST for three periods belonging to the validation period: from 1 June to 31 December 2009.Average LST at basin scale.

Fig. 18 .
Fig. 18.Simulated (continuous grey line) and satellite-estimated (grey circles) LST plotted with observed net radiation (continuous grey black line).The graphs represent the basin-scale average values.A brief period of 6 days of June 2009 is shown.

Fig. 19 .
Fig. 19.RMSE map layered over the DEM relief.The value for each pixel is calculated as the mean of the RMSE on the whole validation period.

Fig. 20 .
Fig. 20.Comparison of the streamflow simulations obtained with calibration using remote sensing data (Cal.R. S.) and calibration using streamflow data (Cal.Streamflow) for Casalcermelli outlet section for calibration period.

Fig. 21 .
Fig. 21.Comparison of the streamflow simulations obtained with calibration using remote sensing data (Cal.R. S.) and calibration using streamflow data (Cal.Streamflow) for Casalcermelli outlet section for validation period.

Table 2 .
Range of variation of the parameters used for the calibration process and minimum parameter gradient ( Par) used for calibration with streamflow data.

Table 3 .
Set of parameters obtained by the calibration.

Table 4 .
Values of skill estimators for the whole calibration period.
a few precipitation events occurred with negligible ground effects.Consistent rainfall events started in October and the most relevant occurred in November and December.Similarly to what was done for the model calibration, a warm-up

Table 5 .
Values of percentage error for the peak flows (PPE) of the main events that occurred during the calibration.

Table 6 .
Values of skill estimators for the whole validation period.

Table 7 .
Values of percentage error for the peak flows (PPE) of the main events that occurred during the validation period.

Table 8 .
Comparison of satellite estimates and simulated LST.Values of the skill estimators for the period: mean bias (BIAS), rootmean-square error (RMSE), mean absolute error (AE) and correlation coefficient (CORR).

Table 9 .
Characteristics of the synthetic rainfall events related with the two surface parameters and the considered sections.

Table 10 .
Values of parameters u c and u h calibrated using geomorphologic data.
Table 11 in terms of skill estimator values.The performance of the model in terms of streamflow reproduction is good.The skill estimators have overall values similar in respect to those reported in Tables