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.