Journal of Hydrology 287 (2004) 214–236
www.elsevier.com/locate/jhydrol
Using a multiobjective approach to retrieve information
on surface properties used in a SVAT model
J. Demartya,b,*, C. Ottléa, I. Braudc,d, A. Oliosoe,
J.P. Frangib, L.A. Bastidasf,g, H.V. Guptaf
a
Centre d’étude des Environnements Terrestre et Planétaires-CETP/CNRS/UVSQ, 10-12, Avenue de l’ Europe, 78140 Vélizy, France
b
Laboratoire Environnement et Développement-LED/Université Paris 7, Case 7071, 2 place Jussieu, 75251 Paris Cedex 5, France
c
Laboratoire d’étude des Transferts en Hydrologie et Environnement, (LTHE UMR5564 CNRS, INPG, IRD, UJF),
BP 53 38041 Grenoble Cedex 9, France
d
CEMAGREF, UR Hydrologie-Hydraulique, 3bis Quai Chauveau, CP 220, 69336 Lyon Cedex 9, France
e
INRA Avignon, Unité CSE, Domaine Saint Paul, Site Agroparc, 84914 Avignon Cedex 9, France
f
SAHRA-Department of Hydrology and Water Resources, University of Arizona, Tucson, AZ 85721-0011, USA
g
Department of Civil and Environmental Engineering, Utah State University, Logan, UT 84322-8200, USA
Received 28 August 2002; revised 1 October 2003; accepted 17 October 2003
Abstract
The reliability of model predictions used in meteorology, agronomy or hydrology is partly linked to an adequate
representation of the water and energy balances which are described in so-called SVAT (Soil Vegetation Atmosphere Transfer)
models. These models require the specification of many surface properties which can generally be obtained from laboratory or
field experiments, using time consuming techniques, or can be derived from textural information. The required accuracy of the
surface properties depends on the model complexity and their misspecification can affect model performance. At various time
and spatial resolutions, remote sensing provides information related to surface parameters in SVAT models or state variables
simulated by SVAT models. In this context, the Simple Soil-Plant-Atmosphere Transfer-Remote Sensing (SiSPAT-RS)
model was developed for remote sensing data assimilation objectives. This new version of the physically based SiSPAT model
simulates the main surface processes (energy fluxes, soil water content profiles, temperatures) and remote sensing data in the
visible, infrared and thermal infrared spectral domains. As a preliminary step before data assimilation in the model, the
objectives of this study were (1) to apply a multiobjective approach for retrieving quantitative information about the surface
properties from different surface measurements and (2) to determine the potential of the SiSPAT-RS model to be applied with
‘little’ a priori information about input parameters. To reach these goals, the ability of the Multiobjective Generalized
Sensitivity Analysis (MOGSA) algorithm to determine and quantify the most influential input parameters of the SiSPAT-RS
model on several simulated output variables, was investigated. The results revealed the main influential input parameters
according to different contrasted environmental conditions, and contributed to the reduction of their a priori uncertainty range.
A procedure for specifying surface properties from MOGSA results was tested on the thermal and hydraulic soil parameters,
and evaluated through the SiSPAT-RS model performance. Although slightly lower than a reference simulation,
the performance were satisfactory and suggested that complex SVAT models can be driven with little a priori information
* Corresponding author. Address: INRA Avignon-Unité CSE, Domaine Saint Paul, Site Agroparc 84914 Avignon Cedex 9, France. Tel.:
þ 33-432-722427; fax: þ33-432-722362.
E-mail address: jdemarty@avignon.inra.fr (J. Demarty).
0022-1694/$ - see front matter q 2004 Elsevier B.V. All rights reserved.
doi:10.1016/j.jhydrol.2003.10.003
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
215
on soil properties, as in a future context of remote sensing data assimilation. Measurement acquired on a winter wheat field of
the ReSeDA (Remote Sensing Data Assimilation) experiment were used in this study.
q 2004 Elsevier B.V. All rights reserved.
Keywords: Soil vegetation atmosphere transfer models; ReSeDA experiment; Simple soil-plant-atmosphere transfer-remote sensing model;
Multiobjective sensitivity analysis; Multiobjective generalized sensitivity analysis algorithm; Remote sensing data assimilation
1. Introduction
Soil Vegetation Atmosphere Transfer (SVAT)
models describe energy and water exchanges between
the soil, the vegetation and the atmosphere. The
understanding and the quantification of these biophysical processes are important for meteorological,
agronomical and hydrological purposes. SVAT
models require a large set of input parameters and
initial state variables that are spatially and temporally
distributed. In general, model complexity and the
number of parameters increase simultaneously. For
example, the Interactions Soil Biosphere Atmosphere
(ISBA) model (Noilhan and Planton, 1989; Noilhan
and Mahfouf, 1996) used in the operational simulations of the French weather forecast model has
around 10 parameters and can easily be applied at
different spatial scales. On the other hand, the Simple
Soil-Plant-Atmosphere Transfers (SiSPAT) model
(Braud et al., 1995) has been developed as a
physically based research tool taking into account
many biophysical processes, such as detailed profile
of root extraction and coupled heat and water
exchanges in the soil. This model needs the
specification of around 60 parameters and initial
state variables (this number depends on the soil
description). Most of them vary in time and space and
are often assessed through in situ experiments. For
example, hydraulic properties in the soil are required
to establish the relationship between hydraulic
conductivity, soil matric potential, and soil water
content. Soil thermal properties are also required. The
specification of such soil properties affects significantly model behavior, implying that SVAT models
need to be calibrated or regularly corrected.
Generally, calibration in land surface modeling is
considered as the determination of a single optimum
parameter set, allowing the ‘best’ simulation of
several output variables. Due to errors in model structure and parameter uncertainties, the determination
of this set is often impossible. Such a
difficulty was the basis for the development of
multiobjective approaches (Yapo et al., 1998; Gupta
et al., 1998). These approaches are based on statistical
analysis of different objective functions in order to
determine parameter sets that cannot be distinguished
in terms of model performance. The choice and the
number of the objective functions depend on the case
studied, the model and the available data. For instance,
Madsen (2000) used a multiobjective approach to
calibrate a conceptual rainfall-runoff model. A single
output model variable (discharge) was studied but the
author tested and applied different mathematical
objective functions in order to simulate different
streamflow ranges (peak flow, low flow or overall
shape of hydrograph). Gupta et al. (1999) and
Bastidas et al. (1999) applied the multiobjective
approach to the Biosphere Atmosphere Transfer
Scheme (BATS, Dickinson, 1984) SVAT model
using only one mathematical objective function
defined as the Root Mean Square Error (RMSE),
but considering simultaneously four different surface
variables.
The multiobjective calibration has already been
applied in different modeling contexts. Several recent
studies have been devoted to the multiobjective
calibration of rainfall-runoff models (Madsen, 2000;
Boyle et al., 2000, 2001; Houser et al., 2001) and land
surface models (Gupta et al., 1999; Xia et al., 2002;
Leplastrier et al., 2002). Such a calibration requires
the use of a complex optimization algorithm for
retrieving the different parameter sets which lead to a
satisfactory simulation of the surface processes. To be
efficient, the optimization algorithm generally
requires a high number of model runs (Duan et al.,
1992; Yapo et al., 1996). In the case of complex
models, which have a significant simulation time,
this efficiency could be highly impacted.
To reduce the dimension of the parameter
optimization problem, multiobjective sensitivity
216
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
analysis approach has been also investigated in a few
recent studies. For example, the Multiobjective
Generalized Sensitivity Analysis (MOGSA; Bastidas
et al., 1999) algorithm was developed to determine the
main influential input parameters. It has already been
applied with a land surface model (Bastidas et al.,
1999) and a hydrochemical model (Meixner et al.,
1999). Although this algorithm was originally
developed to determine the model sensitivity to its
input parameter, it is also a useful tool for retrieving
quantitative information about influential parameters
(Demarty, 2001). Such a behavior is potentially
interesting for driving complex models which require
a significant simulation time.
Moreover, by providing frequent observations at
different spatial scales, remote sensing is an attractive
tool for retrieving information about surface
properties used in land surface models or to control
model trajectory. Such methods are referred to as
assimilation of remotely sensed data. Even if the
assimilation of remotely sensed data has been widely
used in atmospheric and oceanographic sciences, it
has been implemented only recently in modeling
surface processes. This can be done by coupling the
SVAT model with Radiative Transfer Models (RTM)
in order to simulate both surface processes and
remotely sensed data, such as directional reflectances,
thermal infrared brightness temperatures, passive
microwave brightness temperatures or radar
backscattering coefficients. Thus, Burke et al. (1997,
1998) calibrated soil hydraulic properties from L band
passive brightness temperature measurements in a soil
water and energy budget model, coupled with a
microwave emission model (MICROSWEAT).
Cayrol et al. (2000) assimilated series of AVHRR
surface reflectances and temperatures to calibrate a
coupled SVAT-Vegetation growth model. More
recently, Wigneron et al. (2002) assimilated L band
passive brightness temperature in the ISBA-Ags
model (Calvet et al., 1998), coupled with the
Tau-Omega microwave radiative transfer model, in
order to retrieve initial soil moisture and a parameter
controlling vegetation growth.
This article presents a multiobjective approach
which was performed on the Simple Soil-Plant-Atmosphere Transfer-Remote Sensing (SiSPAT-RS)
model. This new version of the SiSPAT model
(Braud et al., 1995) was developed to facilitate the
assimilation of remotely sensed data. It simulates
soil-vegetation-atmosphere energy and water
transfers as well as remote sensing measurements, at
the field scale, in the visible – near infrared and
thermal infrared spectral domains. As a preliminary
step before data assimilation, the objectives of this
study were (1) to apply a multiobjective approach for
retrieving substantial information about the surface
properties from surface measurements and (2) to
determine the potential of the SiSPAT-RS model to be
applied with ‘little’ a priori information about input
parameters. To reach these goals, this article focuses
on the potential of the MOGSA algorithm to
determine and quantify the most influential input
parameters of the SiSPAT-RS model. Its implementation was tested on a winter wheat field, from the data
base collected during the Alpilles-ReSeDA (Remote
Sensing Data Assimilation) experiment.
This article is divided into six sections. Section
2 provides an overview of the Alpilles-ReSeDA
experiment and the particular data set used in this
study. Section 3 describes the SiSPAT-RS model.
Section 4 presents the multiobjective approach
which was conducted with the MOGSA algorithm
and its implementation. Section 5 describes the
results. Finally, Section 6 opens the article for
discussion and presents some conclusions.
2. The Alpilles-ReSeDA experiment
The aim of the Alpilles-ReSeDA experiment
(Olioso et al., 2002a) was to provide a consistent
dataset for assessing crop and soil processes using
remote sensing data. This experiment was focused on
agricultural land and practices in order to develop
ReSeDA methods to better evaluate soil and vegetation processes (biomass production, crop yield,
energy and water budgets). Therefore, a small
agricultural area, characterized by a large diversity
of crops, was instrumented and monitored during one
year to document the whole crop cycle. The site was
located near Avignon, France, (N438470 , E48450 )
and the experiment lasted from October 1996 to
November 1997. A very flat area, with fields of size
200 m by 200 m characterized the site. Various crops
representative of the area were monitored (wheat,
sunflower, alfalfa).
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
Data acquisition was done at different levels:
(1) ground measurements of meteorological, soil and
vegetation variables on some specific fields in order to
calibrate, test and validate ReSeDA procedures and
(2) remote sensing measurements at different scales
with airborne and space-based instruments over
the whole experimental area in order to extrapolate
the procedures to this area.
In this article, results obtained from a winter
wheat field (field number 101) are investigated.
Available observations used in this study are
described in Table 1. The total simulation period
covers the whole growing and senescent period,
217
between Days Of Experiment (hereafter DOE) 387
(21st January 1997) and 542 (25th June 1997).
Atmospheric forcing was measured at the meteorological site, located in the middle of the
experimental area: data were acquired over a bare
soil surface with a 15 s time step and an averaging
period of 20 min. Soil properties including particle
size data, infiltration and dry bulk density were
measured at several depths along the soil profile.
Canopy height, green Leaf Area Index (LAI)
and Organ Area Index (OAI) for leaves, stems
and ears in green vegetative phase and yellow/
senescent vegetative phase were acquired regularly
Table 1
Information about wheat field 101 measurements over the ReSeDA experiment
Total simulation period
DOE 387 –542 (21st January to 25th June 1997)
Type of field measurements
Atmospheric forcing
Direct incoming solar radiation
Diffuse incoming solar radiation
Incoming atmospheric radiation
Wind speed
Air temperature/vapor pressure
Atmospheric pressure
Precipitation rate
Method/instrument/characteristics
Kipp pyranometer/0.3–3 mm spectral domain
Pyranometer with shadow ring/0.3–3 mm spectral domain
Eppley pyrgeometer/3 –100 mm spectral domain
Vector instruments A100L2 cup anemometer/2 m above the ground
Vaisala HMP35D/2 m above the ground
Electronic barometer
Automatic rain gauge
Soil properties
Soil granulometry profile
Dry bulk density
Retention curves
Saturated hydraulic conductivity
Thermal conductivity
Soil albedo
Laboratory
Laboratory
Pressure chamber method and wind method
Wind method
Laboratory
Vegetation properties
Leaf area index
Cover height
Emissivity
Planimeter
Average sampling
Initial and output variables
Soil moisture profile
Surface soil moisture
Soil water potential profile
Net radiation
Soil heat conduction
Turbulent heat fluxes (two methods)
Brightness temperature
d symbolized the estimated accuracy.
Neutron probes/Solo 40 and 25/0–140 cm
each 10 cm, weekly data
Capacitive probes/HMS 9000/0–5 cm, hourly data
Tensiometers/SKT850C/
Rebs Q7
Transducers Rebs HFT-1 and 3.1/Estimated accuracy 40 W m22
Mono-dimensional Eddy correlation Campbell CA27T=d ¼ 30 W m22
Bowen ratio system/Campbell/d ¼ 60 W m22
In situ radiometers/Heinman KT15/8–14 mm waveband, d ¼ 2 K
218
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
during the crop cycle and interpolated to daily
values. Several initial and modeled output variables
were observed, such as soil moisture down to a
depth of 140 cm, energy and mass fluxes and local
thermal infrared brightness temperatures (Table 1).
Eddy correlation (EC) and Bowen ratio (BR)
methods were implemented to assess latent and
sensible heat fluxes, over short periods and over the
whole crop cycle, respectively. Unfortunately,
failures in the BR method instrumentation, generating error in the system for measuring vapor and
temperature gradients, did not allow the estimation
of turbulent fluxes as accurately as expected.
The comparison with the measurements by EC
method showed an overestimation (bias of
15 Wm22) and a large scatter (root mean square
difference between 50 and 70 W m22) of the
sensible heat flux by the BR system (Olioso et al.,
2002a).
3. Modeling approach
In this study, the SiSPAT model (Braud et al., 1995;
Braud, 2000) with a Remote Sensing (SiSPAT-RS)
module (Demarty, 2001; Demarty et al., 2002) was
used. It was developed specifically to assess the value
of multi-spectral remotely sensed data within a
modeling framework. The ultimate objective of this
work is to improve the prediction of the SVAT
prognostic variables by assimilation of remote sensing
data. For this purpose, the SVAT model was coupled
with two RTM, in the visible – near infrared and
thermal infrared spectral domains, respectively.
3.1. The SiSPAT model
The SiSPAT model describes vertical heat and
water exchanges within the soil-plant-atmosphere
continuum. Following the classification proposed by
Huntingford et al. (1995), SiSPAT is a two-layers
(soil and vegetation) SVAT model (Fig. 1).
Three main modules can be distinguished.
Fig. 1. Schematic representation of the energy processes in SiSPAT model: net radiation (Rn), soil heat conduction (G), sensible and latent heat
fluxes (H and LE), bulk stomatal resistance (rs ) and aerodynamic resistances (rb ; ras ; ra ), temperatures (T) and specific humidity (q). Subscripts
‘s’, ‘v’ and ‘tot’ refer to the soil, vegetation and total contributions, respectively, ‘a’ to the reference height measurements above the canopy and
‘av’ to the fictive level inside the canopy.
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
(i)
The soil module describes vertical coupled heat
and water transfers in the soil. The soil column is
described as a juxtaposition of horizontal horizons, characterized by different thermal and
hydraulic properties. The module solves the
coupled heat and mass transfer equations including vapor phase and water extraction by roots,
following the formalism proposed by Milly
(1982). Numerically, each horizon is divided
into layers with an increased resolution at the
lower and upper boundaries of the horizon. The
dependant variables are the temperature T (K)
and the soil matric water potential h (m) which
are continuous at the interface of soil horizon,
unlike the volumetric water content (Vauclin,
1979). Specifications of retention curve and
hydraulic conductivity parameterizations for
each horizon in terms of volumetric water
content are necessary. Many parameterizations
of retention curve and hydraulic conductivity
have been proposed in the literature (Brooks and
Corey, 1964; Van Genuchten, 1980) and are
incorporated as options in the model.
(ii) The soil-plant module describes the water
extraction by the roots. Water extraction from
the soil is parameterized using a resistance model
developed by Federer (1979), considering that
plant transpiration equals the root extraction.
Two resistances per soil layer (one for the soil to
roots water path and the other one for the water
path through the plant) are computed.
This requires the specification of a total plant
resistance and a description of the root density
profile within the soil column.
(iii) The soil-plant-atmosphere interface module
computes energy fluxes between the soil, the
vegetation cover and the atmosphere.
Atmospheric forcing (Table 1) is imposed at a
reference level Za (2 m here). Net radiation (Rn) is
calculated through a shielding factor (Deardorff,
1978; Taconet et al., 1986) depending on the LAI
and controlling the partition between the
vegetation (Rnv) and the bare soil (Rn s).
The vegetation is assumed to be a semi-transparent ‘big leaf’ and multiple reflections of solar and
long-wave radiation between the soil and the
vegetation are considered. Sensible and latent
heat fluxes (respectively, H and LE) are expressed
219
using an electrical analogy. The computation of a
bulk stomatal resistance together with three
different aerodynamic resistances is necessary.
The bulk stomatal resistance (rs ) represents the
plant physiological response to climatic and
environmental conditions. Following Jarvis
(1976), it is modeled in terms of environmental
factors, incoming solar radiation, vapor pressure
deficit and leaf water potential. The three
aerodynamic resistances describing turbulent
transfers between the soil surface and the canopy
air space (ras ), between the leaf surface and the
canopy air space (rb ) and between the canopy air
space and the atmospheric reference level (ra ) are
determined using wind profile parameterizations
inside and above the canopy developed by
Shuttleworth and Gurney (1990). The atmospheric forcing (incoming radiation, wind
velocity, precipitation, air temperature, and
humidity at the reference level Za ) is required in
order to compute 5 prognostic variables; namely
the soil surface temperature Ts ; the soil surface air
humidity qs ; the air temperature inside the canopy
Tav ; the air humidity inside the canopy qav and the
effective canopy temperature Tv : The model time
step is automatically adjusted: it is prescribed to
2 min and decreases in case of rainfall or in case of
strong temperature or soil matrix potential
gradients; forcing data being linearly interpolated
to each model time step.
3.2. The SiSPAT-RS model
The previous description showed that the RTM
included in the original SiSPAT model were very
simple. This simplicity prevents the computation of
remotely sensed variables in the geometrical
conditions of observation (directional effect) and in
the particular bandwidth (spectral effect) of the
instrument. As a consequence, we decided to couple
SiSPAT with two new canopy RTM in order to make
it possible the simulation of remote sensing variables.
In the visible – near infrared spectral domain
(solar domain between 0.3 and 3 mm), the 2M version
(Multi-layer and Multi-element; Weiss et al., 2001)
of the SAIL model (Scattering by Arbitrary Inclined
Leaves; Verhoef, 1984, 1985) was used. The canopy
is described as several horizontal layers of vegetation
220
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
in a specific phenological phases (green, yellow or
senescent vegetation layer). Moreover, different
vegetative organs, such as leaf, stem or ear, compose
each layer of vegetation. Optical properties can be
either specified from measurements or modeled by the
PROSPECT model (Jacquemoud and Baret, 1990)
for vegetative organs and by the MRPV (Rahman
et al., 1993) or SOILSPECT (Jacquemoud et al.,
1992) models for the soil. Other main inputs of
2M-SAIL model are the different OAI and the spectral
densities of the solar and atmospheric radiations.
Finally, the scheme simulates the bi-directional
reflectances in the spectral and geometrical measurement conditions.
The 2M-SAIL model can also be used to
compute the radiative exchanges between the soil,
the vegetation and the atmosphere. For that, the
whole solar spectral domain is decomposed into 55
ranges on which the RTM is applied. Surface
albedo, absorbed radiation and diffuse radiation for
each vegetation layers are computed by convolution
over the whole spectral domain. For that, the
implementation of the RTM in the SVAT model
required to account for the different vegetative
layers into the SiSPAT structure. That was simply
done by the specification of an OAI to each
vegetative organ, instead of the LAI as in the
original version. However, we assumed that each
vegetative layers had a different contribution to the
surface processes according to their associated
phenological states. We considered that only
green vegetation layers contributed to plant
transpiration (through the bulk stomatal resistance)
and that total vegetation, including green layers and
yellow/senescent layers, contributed to radiative
(through the shielding factor) and aerodynamic
(through the aerodynamic resistances) processes.
In the thermal infrared domain, the original RTM
implemented in the SiSPAT model simulated the
canopy radiative temperature, which was representative of the whole long-wave spectral domain
(3 – 100 mm) and of the superior hemisphere. In fact,
thermal infrared remote sensing provide brightness
temperature in a short bandwidth (generally between
8 and 14 mm) and in a specific viewing configuration.
As a consequence, a directional parameterization
(François, 2001) was coupled with the SiSPAT model.
It simulated the brightness temperature of the canopy
accounting for the atmospheric incoming long-wave
radiation and the emitted radiations by the soil and by
the vegetation. The incoming long-wave
radiation coming from the atmosphere was prescribed
such as a model input. It was calculated in the
bandwidth of the instrument using the parameterization proposed by Olioso et al. (1995). The thermal
infrared radiations emitted, respectively, by the bare
soil and the vegetation were calculated by a numerical
integration of the Planck’s law on the bandwidth of
the instrument. These integrations required the
temperature of the soil surface and the temperature
of the vegetation which were both simulated by the
SiSPAT model.
4. Multiobjective methodology
A multiobjective framework was established in
order to progress toward the objectives of retrieving
information about the surface properties used in the
SiSPAT-RS model. The framework combined a
sensitivity analysis and a parameter retrieval
procedure. Unfortunately, the use of a complex
optimization algorithm was not consistent with the
large computer time required by the SiSPAT-RS
model. An alternative procedure involved evaluation
of the ability of the MOGSA algorithm (Bastidas et al.,
1999) to detect (model sensitivity) and quantify
(parameter calibration) the main input parameters
which lead to a satisfactory simulation of the surface
processes. This was done through two main steps.
First, in order to better understand model response to
the environment and to analyze the relative impact of
each model parameter and each initial condition,
MOGSA was used on three separate simulation periods
of the wheat cycle. These periods differed in terms of
climatic and phenological conditions (Table 2). To our
knowledge, the sensitivity of a complex SVAT
model coupled with RT models and the impact of
environmental conditions on the model sensitivity
have not previously been studied.
Second, we investigated the ability of the MOGSA
algorithm to retrieve quantitative information about
the surface properties from surface states variables.
In this context, we focused our analysis on the
hydraulic and thermal soil properties. These properties which highly affect the model behavior, are
221
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
Table 2
Specific information about the three periods retained for the SiSPAT-RS sensitivity analyses
Period
Climatic/environmental conditions
Measurements
DOE 440–460
Regular soil drying
Wheat-growing phase
Dry conditions
Wheat well-developed phase
Precipitation events (519– 522)
Regular soil drying (523– 530)
Wheat senescence phase
Eddy correlation method
DOE 505–517
DOE 518–530
usually estimated through in situ measurements.
Their retrieval via a modeling approach represents
an important research issue (Burke et al., 1997).
As MOGSA was able to provide an a posteriori
uncertainty range for each of the most sensitive
parameters (Section 4.3), this study suggests a
procedure to retrieve the hydraulic and thermal
properties of the soil. This procedure was evaluated
through (1) the comparison of retrieved soil
properties to experimental values obtained during
the Alpilles-ReSeDA experiment and (2) the analyze
of SiSPAT-RS model performance on the whole
period of simulation (Section 4.4).
4.1. Theoretical background of multiobjective
approach
Generally speaking, the multiobjective problem
can be stated as the following minimizing formulation
Min{F1 ðuj Þ; F2 ðuj Þ; …; Fm ðuj Þ}
ð1Þ
where Fi symbolizes a single objective function with
i ¼ 1; …; m and uj ¼ {uj1 ; uj2 ; …; ujp } a particular set of
p parameters included in the feasible parameter space.
For an objective function defined as the RMSE,
the following relation is used
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
u X
u1 n
j
ð2Þ
Fi ðu Þ ¼ t
ðZ ðuj Þ 2 Oi;t ðuj ÞÞ2
n t¼1 i;t
where Zi represents the modeled variable, Oi the
observations and n the total number of available
observations. Observations can be obtained from
measurements or from a reference simulation.
The solution of the multiobjective formulation
given in Eq. (1) does not lead to a unique solution,
Bowen ratio method
Bowen ratio method
but to a set of solution, generally named Pareto set or
‘behavioral’ set. Following Gupta et al. (1998),
the solutions of the Pareto set have the two following
properties:
1. For all non-members uj there exists at least
one member uk where Fi ðuk Þ , Fi ðuj Þ for all i ¼
1; 2; …; m:
2. It is not possible to find uj within the Pareto set such
that Fi ðuj Þ , Fi ðuk Þ for all i ¼ 1; 2; …; m:
The first property implies that the feasible parameter space can be partitioned into two parts:
‘behavioral’ solutions (Pareto set) and ‘ non-behavioral’ solutions (remaining set). The second
properties show that, as the individual objectives
functions are non-commensurable, it is not possible to
isolate objectively the best solution among the
‘behavioral’ solutions: each solution in the Pareto set
improves one or several criterion while causing
deterioration of another. Therefore, within the Pareto
set, none of the solutions is better than the other in
terms of all objective functions. A simple example for
only 2 objective functions is presented on Fig. 2.
Points within the criterion space represent the results
obtained from a sample of the feasible parameter
space. The thick line joins all the solutions in the Pareto
set. Points A and B represent particular solutions that
minimize the individual criteria F1 and F2 ;
respectively. By definition, these two particular points
are always members of the Pareto set. Point C is
included in the Pareto set too, because it is impossible
to find another element in the Pareto set providing
simultaneously better results for each individual
criteria F1 and F2 : Conversely, point D doesn’t belong
222
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
constrain the model calibration. The computation of
the Kolmogorov –Smirnov (KS) two samples test was
used as in Bastidas et al. (1999). The test estimates
whether the cumulative marginal distribution
functions of the two regions are different through
the following relation
dM;N ¼ sup ½PB ðxÞ 2 PB ðxÞ
ð3Þ
x
Fig. 2. Example of determination of the Pareto set and Pareto
ranking in a simple case of two criteria {F1 ; F2 }. Each point
represents a particular simulation realized. A maximum rank of 10
has been found. Partition of the sample into acceptable solutions
(linked by dark lines) and non-acceptable solutions (linked by dash
lines) has been realized through a Pareto rank of 2 as threshold.
to the Pareto set because point C provides better results
for each individual criterion F1 and F2 :
4.2. Model sensitivity to input parameter
The MOGSA methodology proposed by Bastidas
et al. (1999) was applied. It is a robust and efficient
method accounting for the joint multiparameter and
multiresponse interactions. Using a Monte Carlo
search of the feasible space, it is based on the notion
of Pareto ranking (Goldberg, 1989). First, it
determines the Pareto set of the entire sample and
assigns to all its members a rank 1 (thick solid line on
Fig. 2). Setting aside these solutions, a new Pareto set
is determined for the remaining solutions, and the new
points determined are assigned to rank 2 (solid line on
Fig. 2). Then, the same procedure is applied until all
the points of the sample have been processed.
Lower ranks provide better results in a multiobjective
sense. Fig. 2 shows results of Pareto ranking for the
preceding example. In this case, a maximum rank of
10 has been found. The choice of a rank as threshold
allows the partition of the sample into the so-called
‘acceptable’ and ‘non-acceptable’ regions. The statistical analysis of each parameter distributions
between these two regions allows both an assessment
of parameter sensitivity and a reduction in the
associated uncertainty ranges which are required to
where dM;N represents the maximum distance between
the two cumulative distribution functions PB and PB ;
which are defined for the N points of the acceptable
region and the M points of the non-acceptable region,
respectively. Thereafter, the KS test associates a
probability value to dM;N : Significance levels for
probability value have been chosen to define ‘high’,
‘medium’ and ‘low’ parameter sensitivities as in
Bastidas et al. (1999).
Fig. 3a shows an example of cumulative
distribution functions obtained for a high
sensitive parameter named X: The feasible space
of X is [0.3; 2.0]. The dark and dotted lines indicate
the cumulative distribution functions obtained for
acceptable solutions and non-acceptable solutions,
respectively. The vertical dark line represents the
maximum distance dM;N :
4.3. Towards calibrated parameters
In addition to the determination of the main
influential parameters on the simulated variables,
the MOGSA algorithm is potentially useful for
retrieving quantitative information about these parameters (Demarty, 2001). In the case where the
cumulative distribution function of acceptable solutions increases faster than cumulative distribution
function of non-acceptable solutions (as on Fig. 3a),
lower values of X provide, in general, better results on
the considered objective functions than higher values.
Fig. 3b shows this behavior in terms of distribution
functions where a X value of 1.05 appears to be
statistically the upper limit for the acceptable
solutions. As a consequence, a parameter estimation
can be carried out using a value in this preferential
range ([0.3; 1.05] for X). It is important to note,
however, that the higher values of X do not provide
systematically worst results. It is a consequence of
parameter interactions in model structure.
223
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
Fig. 3. Comparison of cumulative distribution functions (a) and distribution functions (b) of parameter X obtained for ‘acceptable’ solutions
(solid line) and ‘non-acceptable’ solutions (dashed line) in multiobjective sensitivity analysis. Vertical bar represents the maximum distance
between the two cumulative distribution functions on which Kolmogorov–Smirnov test is applied.
4.4. Modeling implementation
The sensitivity of the SiSPAT-RS model was first
investigated with the MOGSA algorithm.
Three specific time periods over the crop cycle of
the winter wheat field 101 were considered,
according to the available ground observations and
to the various environmental conditions and
vegetation stages (Table 2). The first period (DOE
440 –460; March 15 to April 4, 1997), corresponded
to a regular soil drying, without precipitation, during
the wheat-growing period. A very dry soil and a
well-developed canopy characterized the second
period (DOE 505 –517; May 19– 31). The third
period (DOE 518– 530; June 1– 13) corresponded to
the senescent phase of wheat, with two precipitation
events (DOE 519 and 523) preceding a new drying
phase. For each period, MOGSA results were
analyzed in terms of relative parameter sensitivity.
Five variables were considered in MOGSA as single
criteria: soil heat conduction flux {G}, sensible heat
flux {H}, latent heat flux {LE}, water content of the
upper five soil centimeters {u0 – 5 }, and local
directional brightness temperature {Tb }. Objective
functions associated with each individual criteria
were defined as the RMSE between observed and
model-simulated variables, from 7 a.m. to 4 p.m.
The second objective of this article was to
determine the ability of the MOGSA algorithm for
retrieving quantitative information about thermal and
hydraulic soil properties. Thus, following the methodology presented in Section 4.3, the preferential
range of the most influential parameters were first
deduced for each of the three previous periods.
Results were analyzed in comparison with the
observed values during the Alpilles-ReSeDA
experiment. To better understand the impact of soil
parameters within the model structure, a very simple
procedure for specifying soil properties from MOGSA
results was tested through a SiSPAT-RS model
simulation on the whole crop period, between DOE
387 and DOE 542. It consisted of attributing a middle
value within the preferential range. A simple rule was
established to prescribe parameters for which no
model sensitivity was obtained. To emphasize the
performance of the simulation, a comparison simulation based on the soil parameters determined during
the Alpilles-ReSeDA experiment was also performed.
The two simulations were finally compared using the
three following statistical criteria: the Bias (Eq. (4)),
the RMSE (Eq. (2)) and the Efficiency (Eq. (5))
Bi ðuÞ ¼ Zi ðuÞ 2 Oi ðuÞ
n
X
ð4Þ
ðZi;t ðuÞ 2 Oi;t ðuÞÞ2
Ei ðuÞ ¼ 1 2 t¼1
n
X
t¼1
ð5Þ
2
ðZi;t ðuÞ 2 Zi ðuÞÞ
224
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
i are, respectively, the simulated and
and O
where Zi
observed average values of the output i: All the other
parameters were prescribed to values calibrated or
determined during the Alpilles-ReSeDA experiment
(Olioso et al., 2002b).
For all SiSPAT-RS simulations, a soil column of
200 cm was layered in 3 horizons, corresponding to
the 0 – 10 cm soil depth (denoted H1 hereafter),
10 –30 cm (H2) and 30 –200 cm (H3) soil layers.
This soil description appeared relatively realistic as a
first approximation, for crops, according to agricultural practices (Chanzy, personal communication).
The retention and hydraulic conductivity curves were
parameterized using the Van Genuchten (1980) model
and the Brooks and Corey (1964) model, respectively.
Fuentes et al. (1992) showed that this combination
was better fulfilling mathematical constraints imposed
by static and dynamical considerations. The vegetation was modeled with two layers: a green
vegetation layer and a yellow-senescent vegetation
layer. The sensitivity of the SiSPAT-RS model to 50
parameters (23 for soil, 27 for vegetation) and 10
initial variables were investigated in the MOGSA
algorithm. An uncertainty range was attributed to each
parameter and initial variable (Table 3). For hydrodynamic soil parameters and initial variables, variability of the entire ReSeDA area was retained. For
optical properties, a simple shifting of the entire
spectra was considered for wet soil, dry soil and
leaves. For OAI and roots parameters, a function
described by two parameters (controlling, respectively, the amplitude and the spreading of the function)
was used to account for temporal variations during the
crop cycle. Associated uncertainty ranges are quite
representative of errors occurred in the measurement.
Finally, a uniform distribution was associated to each
parameter uncertainty range. The exception is the
saturated hydraulic conductivity for which a uniform
distribution was associated to the logarithm of the
uncertainty range.
5. Results and discussion
5.1. Multiobjective sensitivity analysis
Samples of maximum 2500 simulations were
realized for the three studied periods. Tests were
performed to verify that the size of the samples was
large enough to obtain robust results. They showed
that a number of 1500 simulations was the
minimum sample size necessary both to stabilize the
number of sensitive parameters and to eliminate
the impact of initial parameter sampling on the
sensitivity analysis results. Furthermore, a Pareto
rank 2 was chosen as threshold for the first specific
period, while a Pareto rank 1 was considered for the
two other ones. This choice ensures a compromise
between the number of sensitive parameters and a
reasonable ratio between the number of acceptable
and non-acceptable simulations.
Figs. 4 – 6 present sensitivity analysis results
obtained on each period. Each figure is composed
by two subplots showing the relative sensitivity of 30
particular parameters or initial variables. A vertical
bar indicates the relative sensitivity in terms of
probability result of the KS test. The horizontal
dashed lines indicate the transition levels between
‘high’, ‘medium’ and ‘low’ sensitivity. If the vertical
bar crosses the level 0.01, the parameter is considered
to have a high sensitivity. If the vertical bar crosses
only the level 0.05, the parameter is considered to
have a medium sensitivity. Under the level 0.05,
the parameter is insensitive.
Concerning the first period (Fig. 4), results show 18
sensitive parameters, 13 of which with ‘high’
sensitivity and five with ‘medium’ sensitivity.
Moreover, eight parameters were related to the soil
properties (Wsa1, hg2, hg3, Ks2, Ks3, La1, La2, Emis),
five to the vegetation (Rsm, Rp, PFC, Emif, Alaiv)
and five to the root profile (Azrpm, Azrt, Apmr, Lzrt,
Lfdr). As might be expected from the knowledge of
the crop development and the drying of the soil, we
found a significant impact of parameters controlling
stomatal regulation (Rsm, Rp, PFC, Alaiv) and water
exchanges in the deep soil layers (hg2, hg3, Ks2, Ks3,
Azrpm, Azrt, Apmr, Lzrt, Lfdr). Parameters La1 and
La2 controlled the thermal conductivity in soil layers
between 0 and 30 cm (H1 and H2 ). They played a role
on the computation of the soil heat conduction flux
G. The soil and vegetation emissivities (Emis, Emif)
were high sensitive parameters. This was due to their
high impact on the brightness temperature Tb :
Results obtained for the second period indicate 18
sensitive parameters (Fig. 5), 10 of which had ‘high’
sensitivity and eight ‘medium’ sensitivity. Of these,
225
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
Table 3
List of parameters and initial variables used in the MOGSA algorithm and their considered uncertainty ranges
Name
Van Genuchten retention curve
1
Wsa1
2
Wsa2
3
Wsa3
4
hg1
Description (units)
Saturated water content for H1 (m3 m23)
Saturated water content for H2 (m3 m23)
Saturated water content for H3 (m3 m23)
Scale factor in the VG retention
curve model for H1 (m)
5
hg2
Scale factor in the VG retention
curve model for H2 (m)
6
hg3
Scale factor in the VG retention
curve model for H3 (m)
7
n1
Shape parameter in the VG retention
curve model for H1 (–)
8
n2
Shape parameter in the VG retention
curve model for H2 (–)
9
n3
Shape parameter in the VG retention
curve model for H3 (–)
Hydrodynamic and thermal soil conductivities
10
Ks1
Saturated liquid hydraulic conductivity
for H1 (m s21)
11
Ks2
Saturated liquid hydraulic conductivity
for H2 (m s21)
12
Ks3
Saturated liquid hydraulic conductivity
for H3 (m s21)
13
La1
Multiplicative coefficient of thermal
conductivity for H1 (– )
14
La2
Multiplicative coefficient of thermal
conductivity for H2 (– )
15
La3
Multiplicative coefficient of thermal
conductivity for H3 (– )
16
Cd1
Volumetric heat capacity of minerals
for H1 (J m23 K21)
17
Cd2
Volumetric heat capacity of minerals
for H2 (J m23 K21)
18
Cd3
Volumetric heat capacity of minerals
for H3 (J m23 K21)
Stomatal conductance
19
Rsm
Minimal stomatal resistance (s m21)
20
RsM
Maximal stomatal resistance (s m21)
21
Rp
Total plant resistance (s m21)
22
PFC
Critical leaf water potential (m)
23
PST
VDP stress function parameter in the
stomatal resistance (Pa21)
Optical properties
24
Decv
Shifting of green leaves spectrum ( –)
25
Decj
Shifting of yellow leaves spectrum ( –)
26
Emif
Leaf emissivity ( –)
27
Ass
Shifting for dry albedo spectrum ( –)
28
Ash
Shifting of wet dry albedo spectrum ( –)
29
Wd
Maximum 0–5 cm water content
for dry soil albedo (m3 m23)
A priori uncertainty range,
All periods
0.37–0.53
0.37–0.43
0.37–0.40
22.0 –20.02
24.0 –20.6
25.0 –22.5
2.116–2.151
2.115–2.151
2.110–2.144
5.0 £ 10210 –5.0 £ 1026
2.0 £ 10210 –2.0 £ 1026
2.0 £ 10210 –2.0 £ 1026
0.3333–4.0
0.3333–4.0
0.3333–4.0
0.75 £ 106 –1.25 £ 106
0.75 £ 106 –1.25 £ 106
0.90 £ 106 –1.50 £ 106
25 –160
3500–7000
5 £ 1011 –5 £ 1012
2110 –2170
1.0 £ 1024 –5.0 £ 1024
20.05–0.05
20.05–0.05
0.96–0.99
20.04 –0.04
20.04 –0.04
0.15–0.24
(continued on next page)
226
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
Table 3 (continued)
Name
Description (units)
A priori uncertainty range,
All periods
30
Ww
0.25–0.33
31
32
Emis
Gv
33
Gth
Minimum 0–5 cm layer water content
for wet soil albedo (m3 m23)
Soil emissivity ( –)
Visible–infrared parameter
of the shielding factor ( –)
Thermal infrared parameter
of the shielding factor ( –)
Initial variables
34
35
36
37
38
39
40
41
42
43
Ts
T2.5
T10
T50
T110
Ws
W2:5
W10
W50
W110
Canopy structure
44
Alaiv
45
Alaij
46
47
Aht
Liaiv
48
Liaij
49
Lht
Root profile
50
51
52
53
Zri
Zrmu
Zrmd
Azrpm
54
Azrt
55
Apmr
56
Afdr
57
Lzrpm
58
Lzrt
59
Lpmr
60
Lfdr
0.94–0.98
0.40–0.60
0.7 –0.95
Initial surface temperature (K)
Initial soil temperature at 2.5 cm (K)
Initial soil temperature at 10 cm (K)
Initial soil temperature at 50 cm (K)
Initial soil temperature at 110 cm (K)
Initial surface water content (m3 m23)
Soil water content at 2.5 cm (m3 m23)
Soil water content at 10 cm (m3 m23)
Soil water content at 50 cm (m3 m23)
Soil water content at 110 cm (m3 m23)
282.1–285.1/287.0–291.2/287.0–291.2
282.9–285.9/289.7–293.3/289.7–293.3
283.2–286.1/290.6–294.1/290.6–294.1
283.4–285.0/288.7–291.7/288.7–291.7
283.2–284.5/286.9–289.4/286.9–289.4
0.20–0.30/0.09–0.14/0.09–0.14
0.20–0.24/0.13–0.17/0.13–0.17
0.22–0.25/0.19–0.23/0.19–0.23
0.31–0.33/0.23–0.26/0.23–0.26
0.37–0.39/0.275–0.29/0.275– 0.29
Amplitude coefficient of green organ
area index (m2 m22)
Amplitude coefficient of yellow organ
area index (m2 m22)
Amplitude coefficient of cover height (m)
Spreading coefficient of green organ
area index ( –)
Spreading coefficient of yellow organ
area index ( –)
Spreading coefficient of cover height
width parameter ( –)
0.6 –1.4/1.2–2.5/1.2–2.5
1st length of the root density profile (m)
2d length of the root density profile (m)
3rd length of the root density profile (m)
Amplitude coefficient of 4th length
of the root density profile (m)
Amplitude coefficient of 5th length
of the root density profile (m)
Amplitude coefficient of parameter pmr
of the density profile ( –)
Amplitude coefficient of maximum root
length density (m m23)
Spreading coefficient of 4th length
of the root density profile ( –)
Spreading coefficient of 5th length
of the root density profile ( –)
Spreading coefficient of parameter pmr
of the density profile ( –)
Spreading coefficient of maximum
root length density ( –)
0.2 –0.5/0.9–1.8/0.9–1.8
0.55–0.82
2 £ 1024 –5 £ 1024/2.2 £ 1023
– 4.5 £ 1023/2.2 £ 1023 –4.5 £ 1023
1 £ 1023 –3 £ 1022/4.6 £ 1023
– 2 £ 1022/4.6 £ 1023 – 2 £ 1022
2 £ 1024 –3 £ 1024
0.0225–0.0375
0.0375–0.0625
0.1125–0.1875
0.51–0.85
1.26–2.1
0.36–0.6
5666–8500
1.02 £ 1024 –2.5 £ 1024
1.25 £ 1024 –2.5 £ 1024
1.5 £ 1024 –4.5 £ 1024
1 £ 1024 –2 £ 1024
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
227
Fig. 4. Multiobjective sensitivity analysis results on period 440 –460 for the 60 parameters and initial variables of SiSPAT-RS model. Vertical
bars indicate relative sensitivity of parameters in terms of probability result of the KS test. Horizontal dashed lines indicate transition levels
between ‘high’ (above 0.01), ‘medium’ (between 0.01 and 0.05) and ‘low’ (under 0.05) sensitivities.
seven parameters concerned the soil properties,
notably four for the surface horizon (Wsa1, hg1, n1 ;
La1) and three for the deepest horizon (Wsa3, hg3, n3 ).
The predominance of these parameters was consistent
with a well-developed canopy and the drying
conditions. In such dry conditions for instance,
parameter n became the predominant parameter of
the retention curve and of the hydraulic conductivity
function. This explained the high model sensitivity
which was observed on this parameter, instead on the
saturated hydraulic conductivity as in the first period.
Moreover, the model was sensitive to the parameters
controlling water exchanges near soil surface as well
as in deeper soil layers, according to the high
evaporative demand of the atmosphere. This was
emphasized by the five sensitive vegetation
parameters: three of them controlled the stomatal
regulation (Rsm, Rp, PFC), while the two others
concerned the canopy structure (Alaiv, Llaiv).
For roots, the key parameter controlling the amplitude
of the maximum root depth (Azrt) was highly
sensitive. Correlated to the dry conditions, initial
water contents in the soil had a significant impact on
the simulated processes.
In the third case, sensitivity analysis results for the
senescent period (Fig. 6) show 11 sensitive
parameters, nine of which having ‘high’ sensitivity
and two ‘medium’ sensitivity. For this case, only soil
parameters of the surface horizon were concerned for
the soil (Wsa1, hg1, n1 ; La1). No sensitivity was
observed for stomatal regulation. These results were
explained by the precipitation and the low
228
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
Fig. 5. Multiobjective sensitivity analysis results on period 505 –517 for the 60 parameters and initial variables of SiSPAT-RS model. See
legend on Fig. 6.
transpiration. Moreover, the high model sensitivity to
green and yellow OAI (Alaiv, Alaij, Llaiv) and to
cover height (Aht) were explained by the strong
impact of canopy structure on the simulation.
Leaf emissivity (Emif) had again a great impact
because of its effect on the simulation of the
brightness temperature Tb : Finally, except for the
Lpmr parameter for which a ‘medium’ sensitivity was
found, no significant sensitivity was observed for root
parameters during the senescent period.
Results obtained above clearly revealed model
sensitivity to different parameters all over the
vegetation cycle. We found that only three parameters
were associated with a high sensitivity during the
three specific periods: the saturated water content of
H1 (Wsa1), the multiplicative coefficient of the thermal
conductivity of H1 (La1) and the green OAI (via
Alaiv). Wsa 1 and La 1 had a high impact on
the simulations of the surface water content u0 – 5
and {G}, respectively. On the other hand, the
parameter Alaiv had an impact on the computation
of many processes, such as partition of incoming
radiative energy between soil and vegetation, turbulent transfers and evapotranspiration. All of the other
model sensitivities which were observed, were
correlated with the specific climatic and environmental conditions. Moreover, an high sensitivity was
observed for the hydraulic and thermal properties of
the soil. These soil properties are generally estimated
through laboratory methods and field measurements,
for a result which can be highly affected by modeler
subjectivity and measurement errors. Results obtained
on the three elementary short periods revealed the
ability to specify few of them without a priori
information. The two following sections describe
results obtained in this context.
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
229
Fig. 6. Multiobjective sensitivity analysis results on period 518 –530 for the 60 parameters and initial variables of SiSPAT-RS model. See
legend on Fig. 6.
5.2. Analyze of the preferential uncertainty ranges
of the soil properties
Following the methodology described in Section
4.3, the preferential uncertainty ranges of the sensitive
soil parameters were retrieved for the three specific
periods of simulation (Table 4). It is important to note
that each period provided complementary information
according to the climatic and environmental
conditions, and that no contradiction in the preferential ranges were observed. Such results showed that
the choice of three separate simulation periods was
pertinent. More specifically, the MOGSA algorithm
returned substantial information about several soil
properties, especially the Van Genuchten retention
curves. For example, the a priori uncertainty range of
the saturated water content of the first horizon Wsa1 ;
which was between 0.37 and 0.53 m3 m23, was highly
reduced to between 0.37 and 0.44 m 3 m 23.
For hydraulic conductivities, however, the preferential uncertainty ranges can still be relatively large, as
for the third horizon which ended up between
5 £ 1029 and 2 £ 1026 m s21. An other point concerned the comparison of the preferential uncertainty
ranges with the derived values from the AlpillesReSeDA experiment. Except for the n1 parameter,
all of the observed values were included within the
preferential range. So, to better analyze the impact of
the soil properties in the SiSPAT-RS performance, a
simple specification procedure were applied.
5.3. Parameter specification of the soil properties
The specification of soil properties simply
consisted in prescribing a near middle value of the
preferential range (Table 4). As no model sensitivity
230
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
Table 4
Prescribed thermal and hydraulic soil properties in the SiSPAT-RS model. Values were deduced from preferential parameter ranges and
observations made during ReSeDA experiment for Rspe and Rref model simulations, respectively
Sensitive parameter
Preferential ranges
Name (unit)
Period 1
Specified value (Sspe )
Prescribed value (Sref )
Period 2
Period 3
Whole period
Whole period
Van Genuchten retention curve
0.37–0.45
Wsa1 (m3 m23)
/
Wsa2 (m3 m23)
Wsa3 (m3 m23)
/
/
hg1 (m)
hg2 (m)
22.2–20.6
hg3 (m)
24.7–23.0
n1 ( –)
/
n2 ( –)
/
n3 ( –)
/
0.37–0.45
/
0.37–0.39
20.5 –20.02
/
23.5 –23.0
2.140–2.151
/
2.130–2.140
0.37–0.44
/
/
20.5 –20.02
/
/
2.140–2.151
/
/
0.42
0.40
0.38
20.3
21.0
23.3
2.145
2.140
2.135
0.43
0.41
0.38
20.4
20.8
23.0
2.129
2.133
2.136
Hydraulic and thermal soil conductivities
/
Ks1 (m s21)
Ks2 (m s21)
2 £ 1028 –2 £ 1026
21
Ks3 (m s )
5 £ 1029 –2 £ 1026
La1 ( –)
0.333–0.8
La2 ( –)
0.7–4.0
La3 ( –)
/
Cd1 (J m23 K21)
/
Cd2 (J m23 K21)
/
Cd3 (J m23 K21)
/
/
/
/
0.333–0.7
/
/
/
/
/
/
/
/
0.333–0.6
/
/
/
/
/
5 £ 1027
2 £ 1027
5 £ 1028
0.5
2.0
1.0
1.0 £ 106
1.0 £ 106
1.20 £ 106
5 £ 1027
2 £ 1027
5 £ 1028
0.5
1.0
1.0
1.0 £ 1026
1.0 £ 1026
1.0 £ 1026
was observed for three parameters, a priori rules were
deduced from physically based considerations.
Thus, the saturated water content and the scale factor
for the second horizon (Wsa2 and n2 ; respectively)
were deduced from the average of corresponding
values for the first and the third horizons. The case of
the saturated hydraulic conductivity of the first
horizon (Ks1) was more problematic and a greater
value than Ks2 was a priori prescribed in order to
account for a lower surface density.
To evaluate the impact of this soil properties
specification, a simulation (hereafter Sspe ) was
performed with the SiSPAT-RS model on the
whole crop cycle (Table 1). A comparison
simulation (hereafter Sref ) was performed using
the thermal and hydraulic soil parameters observed
during the ReSeDA experiment (Braud and Chanzy,
2000; Olioso et al., 2002b). For the two simulations, other commune sensitive parameters which
are listed in Table 5 were derived from Olioso et al.
(2002a,b). Root profile and OAI were also derived
from Alpilles-ReSeDA experiment.
Fig. 7 shows examples of simulated soil water
content results, obtained for 0 –5 (surface), 10 –20,
50– 80 and 0 – 140 cm layers. In complement, Table 6
presents results in terms of the performance obtained
for soil water content, brightness temperatures and
surface fluxes. The Sspe simulation provides better
results than Sref on surface layer water content
(RMSE ¼ 0:040 m3 m23 ; E ¼ 0:44 against RMSE ¼
0:044 m3 m23 and E ¼ 0:33). As u0 – 5 was considered
as one of the single criteria in the multiobjective
sensitivity analysis, the hydraulic and thermal
Table 5
Other prescribed main parameters used in the SiSPAT-RS model
Parameter
Prescribed value
Rsm (s m21)
Rp (s m21)
PFC (m)
Emif (–)
Emis ( –)
Lfdr ( –)
70
10212
2130
0.98
0.965
1.25 £ 1024
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
231
Fig. 7. SiSPAT-RS model performance on the soil water content results for the Sref simulation (solid line) and the Sspe simulation (dashed line).
Available measurements collected with capacitive probes (þ) and neutron probes ( p ) are also indicated.
parameters of the soil seemed to be correctly
specified. Conversely, better results were obtained
for the second horizon (10 – 30 cm) for Sref : The lack
of sensitivity observed for H2 retention curve
parameters did not make it possible to specify
correctly these parameters. After DOE 480, the
soil water content of deeper layers was generally
underestimated for the two simulations (negative
biases), implying underestimation of u0 – 140 : Moreover, u0 – 140 was better simulated for the simulation
Sref ðRMSE ¼ 0:006 m3 m23 ; E ¼ 0:98Þ than for Sspe
ðRMSE ¼ 0:011 m3 m23 ; E ¼ 0:95Þ:
Fig. 8 presents results for the surface fluxes and
brightness temperature for the Sspe simulation. For the
brightness temperature, good agreement was observed
for the wheat-growing period, but the scatterplot over
the whole experiment shows a model underestimation
for the higher values of brightness temperature. It was
observed that the underestimation was mainly due to
the senescent period (not presented here). For sensible
and latent heat fluxes, comparisons with the
measurements acquired by EC and BR systems were
done. The growing period is always well simulated.
The RMSE on HEC and LEEC were close to the errors
232
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
Table 6
SiSPAT-RS model performance for the Rspe and Rref model simulations
Variable
Sspe model simulation
Name (unit)
Biasa
Soil water content
u0 – 5 (m3 m23)
u0 – 10 (m3 m23)
u10 – 20 (m3 m23)
u20 – 30 (m3 m23)
u30 – 50 (m3 m23)
u50 – 80 (m3 m23)
u80 – 120 (m3 m23)
u0 – 140 (m3 m23)
Brightness temperature
Tb (K)
Surface fluxes
Rn (W m22)
G (W m22)
HEddy (W m22)
LEEddy (W m22)
HBowen (W m22)
LEBowen (W m22)
a
0.011
20.005
0.023
0.015
0.021
20.008
20.020
20.005
0.6
214.8
0.62
220.2
11.4
240.7
24.8
Sref model simulation
RMSEa
Efficiency
Biasa
0.040
0.027
0.034
0.023
0.025
0.018
0.025
0.011
0.44
0.87
0.77
0.89
0.81
0.84
0.66
0.95
0.026
0.011
0.009
0.002
0.025
20.003
20.016
20.003
1.9
0.95
38.7
23.6
31.4
31.2
67.1
56.3
0.96
0.72
0.79
0.91
0.65
0.60
RMSEa
Efficiency
0.044
0.025
0.021
0.014
0.028
0.012
0.019
0.006
0.33
0.88
0.91
0.96
0.75
0.93
0.80
0.98
0.6
1.9
0.95
214.9
0.6
218.5
9.9
237.9
21.7
38.8
24.6
29.4
29.8
67.4
57.1
0.96
0.70
0.81
0.91
0.65
0.59
In the unit of the considered variable.
estimated during the Alpilles-ReSeDA experiment
(around 30 W m22). Using observations acquired
with the BR system, poor performance was obtained.
For the HBR flux, a bias of 2 40 W m22 and a RMSE
around of 60 –70 W m22 were observed, for instance.
However, such values were consistent with the
differences between the BR and EC methods which
were observed during the Alpilles-ReSeDA
experiment:, BR system overestimated H by
15 W m22 (bias) compared to EC system. Moreover,
the root mean square difference between the two
systems was between 50 and 70 W m22 (Olioso et al.,
2002a). Finally, as for the brightness temperature, the
main differences between simulation and observations
were during the senescent period. We cannot clearly
establish whether the origins of these discrepancies
are in the model structure. For instance, the assumption that the yellow vegetation layer did not contribute
to the evapotranspiration fluxes still remained quite
simple. However, Demarty (2001) showed that
accounting for the yellow vegetation layer in
the SiSPAT model allowed for an improvement in
the simulation of soil heat flux G and soil water
content u0 – 140 : Inclusion of the 2M-SAIL model in
the SiSPAT structure had a positive impact in this
context.
6. Conclusion
This study presents a multiobjective framework
which was applied to the SiSPAT-RS. This model was
developed to simulate the main surface processes
(energy fluxes, soil water content, temperatures) and
remote sensing observations in the visible– infrared
and thermal infrared spectral domains. For the
visible – infrared wavebands, the Multi-layer and
Multi-element (2M) version (Weiss et al., 2001) of
the Scattering by Arbitrary Inclined Leaves (SAIL;
Verhoef, 1984, 1985) canopy radiative transfer model
was implemented within the SiSPAT model structure.
It allowed (1) the distinction between different
vegetation organs, and possible inclusion of a yellow
vegetation layer, in the radiative transfer processes
and (2) the simulation of the bi-directional
reflectances. For thermal infrared wavebands,
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
233
Fig. 8. SiSPAT-RS model performance on the brightness temperatures and surface fluxes for the Sspe model simulation. Right column shows few
days of simulated time series (solid line) and available measurements ((). Left column presents scatterplots obtained on the entire period of
simulation. Measurement of turbulent fluxes acquired by Eddy correlation and Bowen ratio methods are separately used.
234
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
the simple parameterization proposed by François
(2001) was chosen in order to simulate directional
brightness temperature, using the soil surface temperature and the vegetation temperature which were
simulated by the SVAT model.
The multiobjective framework consisted of applying the MOGSA algorithm proposed by Bastidas et al.
(1999) in order to determine and quantify the main
influential parameters on the simulation of many of
the output variables. The MOGSA algorithm was first
applied for the ReSeDA winter wheat field 101 in
order to determine the model sensitivity for three
separate periods. Differences in model sensitivity to
input parameters between the three periods were
consistent with the climatic and environmental
conditions which were observed for each period. For
instance during the wheat-growing period, high model
sensitivity was obtained for parameters controlling
stomatal regulation and water exchanges of deep soil
layers. For the senescent period, high model sensitivity were found for parameters describing the
canopy structure, but there was no model sensitivity
for stomatal regulation and root parameters. More
generally, analysis of model sensitivity is an important issue for a better understanding of the model
structure. Results found in this study emphasize the
high impact of climatic and environmental conditions.
Moreover, they revealed that several independent
allowed the constraining a complex model with a high
number of parameters. On the other hand, results
indicated that on a short period of simulation, a
maximum number of 18 parameters out of the 60 had
a great impact on the simulation of the five output
variables were considered. Thus, the sensitivity
analysis was useful to reduce the number of
parameters that must be considered for
the calibration of a complex model, such as
SiSPAT-RS. It may be noticed that this method,
based on a Monte Carlo approach, is potentially
applicable to the analysis of the correlation between
input parameters, such as the correlation between the
soil hydraulic properties.
Finally, the ability of the MOGSA algorithm to
retrieve quantitative information about sensitive
parameters was revealed. It was shown that the
retrieved hydraulic and thermal parameters were
relatively close to ‘observed’ values during the
Alpilles-ReSeDA experiment. This analyze allows
the definition of a very simple procedure for
specifying the soil properties from MOGSA results.
This procedure was evaluated by comparing model
performance over the whole experimental season with
a reference simulation. The results obtained show
potential for driving complex SVAT models with
little a priori information on soil properties within an
agronomic context. However, the procedure for
specifying hydraulic and thermal soil properties tested
in this study was relatively simple and partially
subjective. Further approaches are already envisaged,
as for example, an iterative procedure based on the
successive contraction of the parameter uncertainty
ranges. Moreover, in this study, surface fluxes were
used as simple objective functions. They provided
substantial information in the retrieving of the input
parameters. Further works will examine model
calibration in the context of ReSeDA, given that
only brightness temperatures and surface soil water
content will be known.
Acknowledgements
The authors would like specially to thank all the
partners of the ReSeDA program for providing ground
truth data. This work was supported by the
INSU/CNRS French national programs (PNRH and
PNTS). We are pleased to thank Doctor
Pierre Guillevic and anonymous reviewers for very
constructive suggestions.
References
Bastidas, L.A., Gupta, H.V., Sorooshian, S., Shuttleworth, W.J.,
Yang, Z.L., 1999. Sensitivity analysis of a land surface scheme
using multicriteria methods. Journal of Geophysical Research
104 (D16), 19481–19490.
Boyle, D.P., Gupta, H.V., Sorooshian, S., 2000. Toward improved
calibration of hydrologic models: combining the strengths of
manual and automatic methods. Water Resources Research 36
(12), 3663–3674.
Boyle, D.P., Gupta, H.V., Sorooshian, S., Koren, V., Zhang, Z.,
Smith, M., 2001. Toward improved streamflow forecasts: value
of semidistributed modeling. Water Resources Research 37
(11), 2749–2759.
Braud I. SiSPAT, a numerical model of water and energy fluxes in
the soil-plant-atmosphere continuum. Version 3.0, SiSPAT
user’s manual, 107 pp., 2000 (http://www.lthe.hmg.inpg.fr).
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
Braud, I., Chanzy, A., 2000. Soil properties, initial and boundary
conditions for use within SVAT models in the framework of the
intercomparison of SVAT models used in the Alpilles-ReSeDA
project, Alpilles-ReSeDA data base, 43 pp.
Braud, I., Dantas-Antonino, A.C., Vauclin, M., Thony, J.L., Ruelle,
P., 1995. A simple soil-plant-atmosphere transfer model
(SiSPAT) development and field verification. Journal of
Hydrology 166, 213–250.
Brooks, R.H., Corey, A.T., 1964. Corey, Hydraulic Properties of
Porous Media. Colorado University, Fort Collins.
Burke, E.J., Gurney, R.J., Simmonds, L.P., Jackson, T.J., 1997.
Calibrating a soil water and energy budget model with remotely
sensed data to obtain quantitative information about the soil.
Water Resources Research 33, 1689–1697.
Burke, E.J., Gurney, R.J., Simmonds, L.P., O’Neill, P.E., 1998.
Using a modeling approach to predict soil hydraulic properties
from passive microwave measurements. IEEE Transactions on
Geoscience and Remote Sensing 36 (2), 454–462.
Calvet, J.C., Noilhan, J., Roujean, J.L., Bessemoulin, P.,
Cabelguenne, M., Olioso, A., Wigneron, J.P., 1998. An
interactive vegetation SVAT model tested against data from
six constrasting site. Agricultural and Forest Meteorology 92,
73–95.
Cayrol, P., Kergoat, L., Moulin, S., Dedieu, G., Chehbouni, A.,
2000. Calibrating a coupled SVAT-Vegetation growth model
with remotely sensed reflectance and surface temperature—a
cas study for the HAPEX-Sahel grassland sites. Journal of
Applied Meteorology 39, 2452– 2472.
Deardorff, J.W., 1978. Efficient prediction of ground surface
temperature and moisture, with inclusion of a layer of
vegetation. Journal Geophysics Research 83, 1889–1903.
Demarty J., Développement et application du modèle SiSPAT-RS à
l’échelle de la parcelle dans le cadre de l’expérience Alpilles
ReSeDA, University Paris 7 thesis, Méthodes physiques en
Télédétection, December 2001, 220 pp. (in French).
Demarty, J., Ottlé, C., François, C., Braud, I., Frangi, J.P., 2002.
Effect of aerodynamic resistance modeling on SiSPAT-RS
simulated surface fluxes. Agronomie 22, 641–650.
Duan, Q., Gupta, H.V., Sorooshian, S., 1992. Effective and efficient
global optimization for conceptual rainfall-runoff models.
Water Resources Research 28, 1015–1031.
Federer, C.A., 1979. A soil-plant-atmosphere model for transpiration and availability of soil water. Water Resources Research 15
(3), 555–562.
François, C., 2001. The potential of directional radiometric
temperatures for monitoring soil and leaf temperature and
soil moisture status. Remote Sensing of Environment 79,
1–12.
Fuentes, C., Haverkamp, R., Parlange, J.Y., 1992. Parameter
constraints on closed-form soil water relationships. Journal of
Hydrology 134, 117–142.
Goldberg, D.E., 1989. Genetic Algorithms in Search. Optimization
and Machine Learning. Addison-Wesley, Reading, MA.
Gupta, H.V., Sorooshian, S., Yapo, P.O., 1998. Toward improved
calibration of hydrologic models: multiple and non-commensurable measures of information. Water Resources Research 34 (4),
751–763.
235
Gupta, H.V., Bastidas, L.A., Sorooshian, S., Shuttleworth, W.J.,
Yang, Z.L., 1999. Parameter estimation of a land surface
scheme using multicriteria methods. Journal of Geophysical
Research 104 (D16), 19491–19503.
Houser, P.A., Gupta, H.V., Shuttleworth, W.J., 2001. Multiobjective
calibration and sensitivity of a distributed land surface water and
energy balance model. Journal of Geophysical Research 106,
(D24), 33421–33433.
Huntingford, C., Allen, S.J., Harding, R.J., 1995. An intercomparison of single and dual-source vegetation-atmosphere transfer
models applied to transpiration from Sahelian Savannah.
Boundary Layer Meteorology 74, 397–418.
Jacquemoud, S., Baret, F., 1990. PROSPECT: A Model of Leaf
Optical Properties Spectra. Remote Sens. Environ. 34, 75– 91.
Jacquemoud, S., Baret, S., Hanocq, J.F., 1992. Modeling spectral
and bidirectional soil reflectance. Remote Sens. Environ. 41,
123 –132.
Jarvis, P.G., 1976. The interpretation of the variations in leaf water
potential and stomatal conductance found in canopies in the
field. Philosophical Transactions of the Royal Society of
London, Series B 273, 593–602.
Leplastrier, M., Pitman, A.J., Gupta, H.V., Xia, Y., 2002. Exploring
the relationship between complexity and performance in a land
surface model using the multi-criteria method. Journal of
Geophysical Research 107 (D20).
Madsen, H., 2000. Automatic calibration of a conceptual rainfallrunoff model using multiple objectives. Journal of Hydrology
235, 276–288.
Meixner, T., Gupta, H.V., Bastidas, L.A., Bales, R.C., 1999.
Sensitivity analysis using mass flux and concentration. Hydrological processes 13, 2233–2244.
Milly, P.C.D., 1982. Moisture and heat transport in hysteretic
inhomogeneous porous media: a matric head-based formulation
and a numerical model. Water Resources Research 18, 489–498.
Noilhan, J., Mahfouf, J.-F., 1996. The ISBA land surface parameterisation scheme. Global and Planetary Change 13, 145– 159.
Noilhan, J., Planton, S., 1989. A simple parameterization of land
surface processes for meteorological models. Monthly Weather
Review 117, 536 –549.
Olioso, A., 1995. Estimating the difference between brightness and
surface temperatures for a vegetal canopy. Agricultural and
Forest Meteorology 72, 237 –242.
Olioso,A., Braud,I.,Chanzy, A.,Demarty,J., Ducros,Y., Gaudu, J.C.,
Gonzales-Sosa, E., Jacob, F., Lewan, L., Marloie, O., Ottlé, C.,
Prévot, L., Thony, J.L., Autret, H., Bethenot, O., Bonnefond, J.M.,
Bruguier, N., Buis, J.P., Calvet, J.C., Caselles, V., Chauki, H.,
Coll, C., Gouget, R., Jongschaap, R., Kerr, Y., King, C.,
Lagouarde, J.P., Laurent, J.P., Mc Aneney, J., Moulin, S.,
Rubio, E., Weiss, M., Wigneron, J.P., 2002a. SVAT modeling
over the Alpilles-ReSeDA experiment: experimental setup for
monitoring energy and mass transfers. Agronomie 22, 651 –668.
Olioso, A., Braud, I., Chanzy, A., Courault, D., Demarty, J.,
Kergoat, L., Lewan, E., Ottlé, C., Prévot, L., Zhao, W., Calvet,
J.C., Cayrol, P., Jongschaap, R., Moulin, S., Noilhan, J.,
Wigneron, J.P., 2002b. SVAT modeling over the AlpillesReSeDA experiment: comparison of SVAT models, first results
on wheat. Agronomie 22, 597–610.
236
J. Demarty et al. / Journal of Hydrology 287 (2004) 214–236
Rahman, H., Pinty, B., Verstraete, M.M., 1993. Coupled surfaceatmosphere reflectance (CSAR) model 2. Semiempirical surface
model usable with NOAA advanced very high resolution
Radiometer Data. Journal of Geophysical Research 98 (D11),
20791– 20801.
Shuttleworth, W.J., Gurney, R.J., 1990. The theoretical relationship
between foliage temperature and canopy resistance in sparse
crops. Quarterly Royal Journal of Meteorological Society 116,
497–519.
Taconet, O., Carlson, T.N., Bernard, R., Vidal Madjar, D., 1986.
Evaluating of a surface vegetation/model using satellite infrared
surface temperatures. Journal of Climate and Applied Meteorology 50, 239 –243.
Van Genuchten, M.T., 1980. A closed equation for predicting the
hydraulic conductivity of unsaturated soils. Soil Science Society
of American Journal 44, 892–898.
Vauclin, M., Haverkamp, R., Vachaud, G., 1979. Résolution
numérique d’une équation de diffusion non-linéaire.
Application à l’infiltration de l’eau dans les sols non saturés.
Presse Universitaire de Grenoble, Grenoble, France, 183 pp.
Verhoef, W., 1984. Light scattering by leaf layers with application
to canopy reflectance modelling: the SAIL model. Remote
Sensors Environment 16, 125 –141.
Verhoef, W., 1985. Earth observation modeling based on layer
scattering matrices. Remote Sensing of Environment 17,
165 –178.
Weiss, M., Troufleau, D., Baret, F., Chauki, H., Prévot, L., Olioso,
A., Bruguier, N., Brisson, N., 2001. Coupling canopy functioning and radiative transfer models for remote sensing data
assimilation. Agricultural and Forest Meteorology 108,
113 –128.
Wigneron, J.P., Chanzy, A., Calvet, J.C., Olioso, A., Kerr, Y., 2002.
Modeling approaches to assimilating L band passive microwave
observations over land surfaces. Journal of Geophysical
Research, 107.
Xia, Y., Pitman, A.J., Gupta, H.V., Lepastrier, M., HendersonSellers, A., Bastidas, L.A., 2002. Calibrating a land surface
model of varying complexity using multi-criteria methods and
the Cabauw data set. Journal of Hydrometeorology V3 (2),
181 –194.
Yapo, P.O., Gupta, H.V., Sorooshian, S., 1996. Automatic
calibration of conceptual rainfall-runoff models: sensitivity to
calibration data. Journal of Hydrology 181, 23 –48.
Yapo, P.O., Gupta, H.V., Sorooshian, S., 1998. Multiobjective global optimization for hydrologic models. Journal
of Hydrology 204, 83 –97.