Multiobjective calibration and sensitivity of a distributed land surface water and energy balance model

. The feasibility of using spatially distributed information to improve the predictive ability of a spatially distributed land surface water and energy balance model (LSM) was explored at the U.S. Department of Agriculture Agricultural Research Service (USDA-ARS) Walnut Gulch Experimental Watershed in southeastern Arizona. The inclusion of spatially variable soil and vegetation information produced unrealistic simulations that were inconsistent with observations, which was likely an artifact of both discretely assigning a single set of parameters to a given area and inadequate knowledge of spatially varying parameter values. Because some of the model parameters were not measured or are abstract quantities a multiobjective least squares strategy was used to find catchment averaged parameter values that minimize the prediction error of latent heat flux, soil heat flux, and surface soil moisture. This resulted in a substantial improvement in the model's spatially distributed performance and yielded valuable insights into the interaction and optimal selection of model parameters. estimates of saturated areas estimated from remote sensing and topography using the weighting of estimates in a fuzzy set framework.

center of the watershed. Soil types range from clays and silts to well-cemented boulder conglomerates [Kustas et al., 1992], with the surface soil textures being gravelly and sandy loams containing, on average, 30% rock and little organic matter. The mixed grass-brush rangeland vegetation ranges from 20 to 60% in coverage. This rangeland region receives 250 to 500 mm of precipitation annually, typically two thirds of it as convective precipitation during the summer monsoon season [Renard et al., 1993].  [Kustas and Goodrich, 1994].

Eighty
Precipitation is the most important spatial forcing variable in semiarid regions because of its highly variable, convective nature; thus much effort was devoted to deriving spatially distributed precipitation data sets for the Monsoon '90 experiment. A multiquadric-biharmonic interpolation algorithm [Syed, 1994] was used to produce spatially distributed precipitation values for the entire model domain from the available rain gauge data. All other meteorological forcing was assumed to be spatially constant and derived from averaging observations at the eight Metflux stations in place during the experiment. A spatial averaged soil moisture profile was derived from several in situ profile observations ] and used to initialize the LSM on July 22 (Table 1).

Land Surface Model
The land surface water and energy balance model (LSM) [Famiglietti and Wood, 1994] simulates land surface runoff, energy fluxes, and soil moisture dynamics in three layers. This LSM is designed to predict diurnal dynamics of the water and This section examines the hypothesis that the use of spatially distributed soils and vegetation information can help to improve the predictive ability of the LSM. The LSM was used to simulate the land surface water and energy dynamics in a spatially distributed manner for the entire WGEW at a 40 m resolution. Predictions were made hourly from July 22 to August 15, or day of year (DOY) 204 to 228, of 1990. LSM parameter estimates were first specified by using catchment-averaged parameter values based largely on observations. Next, spatially distributed values for the parameters were derived from geographic information system's (GIS) maps.

Spatially Constant Parameters
Spatially constant parameters for the LSM were specified primarily on the basis of observations made during Monsoon '90 [Daughtry et al., 1991;Kustas et al., 1992Kustas et al., , 1994aKustas et al., , 1994bKustas et al., , 1994cStannard et al., 1994;Moran et al., 1994aMoran et al., , 1994bHumes et al., 1994aHumes et al., , 1994bWeltz et al., 1994;Menenti and Ritchie, 1994;Hipps et al., 1994;Spangler, 1969]. The albedo of dry soil was specified by adopting the parameters suggested by Dicla'nson et al. [1993]. Soil profile temperatures measured in three trenches at the Lucky Hills and Kendall sites were analyzed and modeled to determine the temperature diurnal damping depth. The saturation soil moisture was assumed to be close to the porosity calculated from bulk density measurements. The surface saturated hydraulic conductivity was based on porosity using the Kozeny-Carmen equation with a coarse fraction correction, following RawIs et al. [1993]. The depth of the surface zone was specified equal to the maximum reported L band microwave penetration depth [Jackson, 1993] to aid compatibility with remote sensing data. Because of the deep water table at Walnut Gulch, surface processes show no sensitivity to the LSM water table parameters, so these were assigned to arbitrary values. The vegetation stress soil moisture parameter was assumed to be halfway between the saturated and residual soil moisture values, while the wilting point was assumed to be equal to the value of residual soil moisture because desert vegetation rarely wilts. The soil moisture at which evapotranspiration reaches the potential rate was assumed to be 10% below saturation, this being the approximate ratio suggested by Sellers et al. [ 1986]. Values for albedo of dry vegetation for semiarid areas were obtained from Dicla'nson et al. [1993], and the albedo of wet vegetation was assumed to be 5% lower than that for dry vegetation. The stomatal resistance, root activity factor, root density, root resistivity, and critical leaf water potential were not measured; indeed, because many of these parameters were not observable quantities, estimates were used. Values in parentheses produce the "best" Monsoon '90 simulation. Read l e9 as 1 fro.

Spatially Variable Parameters
Spatial distributions of the soil and vegetation parameters were estimated using GIS maps of Walnut Gulch vegetation and soils. The potential rooting depths for each soil series were taken from a recent Walnut Gulch soil survey (D. J. Breckenfeld, unpublished document, 1993). Values of the effective porosity, •b e, residual soil moisture, 0r, and the saturated hydraulic conductivity, Ks, were determined by Mayeux [1995] using values suggested by Rawls et al. [1993] and Bouwer [1966]. J. Berglund, unpublished report, (1995) calculated saturated soil moisture, 0e, for each soil class as overlaps were averaged, then a smoothing algorithm was applied to reduce noise, and isolated runoff sinks were filled. The watershed was delineated and the topographic index was derived using the standard methodology detailed by Beven [ 1995]. A selection of representative spatially variable Walnut Gulch parameters is shown in Figure 1 a-1 d, and the watershed averages of these parameters are reported in Table 1.

Results and Discussion
1 e' 1 fi It is evident that the spatially variable parameters have a large impact on the spatial patterns of the simulation. Further, the soil moisture patterns derived using the spatially variable parameters appear to be unrealistic and do not compare well with observed push broom microwave radiometer (PBMR) surface soil moisture (Figures 1 g-1 h). A series of simulations, where only one spatially variable parameter set was used at a time (the other parameters were set to their estimated catchment-averaged values), were performed to determine which subset of spatial parameters contribute The simulated spatial patterns of surface soil moisture for most to these patterns (the topographic index and the August 7, 1990 (DOY 219), using both the spatially constant precipitation fields were spatially variable in all these and the spatially variable parameters, are shown in Figures simulations). The simulations with spatially variable vegetation parameters perform similarly to the control run, indicating that spatial variation in these parameters has little effect on predictions. All of the simulations using spatially variable soil parameters show distinct polygon patterns. The parameter specifying saturated soil moisture has the most influence on simulated spatial patterns, while those which specify the percentages of sand and clay, the saturated hydraulic conductivity, and the residual soil moisture have a more moderate influence.
The soil characteristics are extremely important, and the vegetation characteristics are less important in regulating soil moisture variability and other surface water and energy components in this watershed. The results of this sensitivity also show that the spatially variable topographic index has very little influence on the simulations; this result is a consequence of the large depth to the water table at Walnut Gulch.
The pattern of enhanced spatial soil and vegetation polygons apparent in the simulations is an artifact of both assigning lumped values of the parameters to discretely partitioned areas as well as the misspecification of parameter values. A more appropriate specification of spatial parameters would reflect their continuous variability in space, as obtained with remote sensing. It is also clear that the specification of the spatial parameter values using scattered observations, table lookup, and physical relationships (as described in section 4.2) result in inadequate parameter specification. These results indicate that the available WGEW parameter observations simply have insufficient accuracy and spatial frequency to realistically constrain the LSM.
Because the simulations using spatially constant vegetation and soil parameters compare well to the PBMR patterns, the subsequent parameter calibration studies reported in section 5 assume spatially constant soil and vegetation parameters across the catchment, leaving only topographic index and precipitation as spatially varying entities.
The use of catchment-averaged parameters in this study is supported by White et al. [ 1997], who found that the difference between the calculated area-averaged surface energy fluxes given for the WGEW by a single point land surface model with aggregate parameters and that given by a distributed array of land surface models was small. This conclusion was tested here; although there was more variability than reported by White et al. [1997] because of the finer spatial resolution, there was overall agreement. The relevance of this finding to the present study is that it supports the use of catchment-averaged parameters and forcing for the WGEW if catchment-averaged surface fluxes are required. This conclusion is only valid when catchment-averaged fluxes are of interest; for cases where downslope flows or water availability in valley bottoms is of interest, then spatially variable parameters and forcing are required. In this study, the spatially distributed LSM is required to allow prediction of spatial patterns of soil moisture at the resolution of the available remotely sensed observations.

Model Parameters
The purpose of this study was to obtain estimates of the nonobserved model parameters that would result in improved simulations of the overall behavioral responses of the watershed. In particular, the desire was to select values for the parameters that minimized the error in predicted soil moisture and sensible, latent, and ground heat fluxes. As discussed by Gupta et al. [1999], this situation involves the problem of finding parameter estimates which can simultaneously minimize several noncommensurable criteria. This section examines the hypothesis that a multiobjective parameter calibration method can be used to improve the estimates of the nonobserved model parameters of the LSM model.

Data Available for Calibration
The primary observations used to calibrate the LSM parameters in this study were hourly time series of catchment-averaged, water balance corrected, latent heat flux during unstable atmospheric conditions, soil heat flux, and gravimetric calibrated near-surface soil moisture measured with resistance sensors (Table 1)

Parameters Calibrated
The 13 parameters listed in Table 2 were identified for calibration. Table 2 also lists the typical range of their values; these ranges were specified larger than might normally be expected in order to allow the calibration procedure some leeway in accounting for model and data error. Moreover, because some of the parameters are nonphysical, the specification of their ranges is somewhat arbitrary. Note that the soil heat flux parameters (temperature of the deep soil layer and penetration of diumal heating) were included in this list even though they were measured with confidence. This is because the soil heat flux was poorly simulated by the LSM owing to the simplicity of its soil heat flux submodel; this suggests that nonphysical values for these parameters may be necessary for the model to perform adequately. Note also that the five soil water retention parameters, which can be estimated from observations, were also included in the calibration, because this greatly improved the simulation of soil moisture and yielded insight into the use of observed soil parameters.  Figure 2b shows the problem plotted in the function space with f• (19) and f2(0) on the axes, and the curve represents the trajectory obtained by varying 19 over its feasible range. Clearly, the objective is to find 0 such that we come as close to the origin • ,f2}={0,0} as possible. Note that the segment of the curve between B and C represents the trade-off region, while other portions of the curve represent inferior solutions. Any one of the trade-off solutions can be selected by the user as "best," depending on which trade-off in matching of model performance is considered most acceptable by the user.

Objective Function
Three objective functions were used to calibrate the model parameters: the hourly mean square error in matching of latent heat flux during unstable atmospheric conditions, soil heat flux, and gravimetric calibrated near-surface soil moisture. In each case the mean square error objective function, F, was computed using the normalized sum of the squared difference between the catchment-averaged, water balance corrected observations, Oi, and the predictions, Pi, whenever the acceptable solution can be selected by the user from this region observation was available; thus of trade-off solutions using subjective judgment.

Inspeclion of these results shows that the five solutions giving
The procedure adopted for the multiobjective calibration of best matching to soil moisture also give relatively good the LSM parameters at Walnut Gulch was to conduct a matching of the other four fluxes. In addition, it was decided Monte-Carlo search of the feasible parameter space indicated that the soil moisture objective was most important for the in Table 2. A total of 200,000 parameter set realizations were purposes of this study. Therefore the parameter set that randomly generated, and the three objective function values produced the lowest soil moisture objective was chosen were computed at each parameter set. From these, the (bolded line in Figure 4). The selected parameter values are "noninferior" trade-off parameter region was estimated by finding those points for which there did not exist at least one other member of the original 200,000 realizations for which all three objective function values were better. Only 125 parameter set realizations were found to belong to this estimated noninferior region, indicating a rather substantial reduction in overall parameter uncertainty. In addition, at each noninferior point, the mean square error in matching of hourly sensible heat flux was also computed for diagnostic purposes.
The estimated noninferior region is displayed in Figure 3 .g., deep soil 5). The range of prediction shows that latent heat flux can be temperature and percent clay), showing that these parameters calculated well by the LSM but at the expense of soil moisture. are identified well. The noninferior parameter values also are By choosing to predict soil moisture most accurately, rather generally confined to moderate ranges of the parameter space. than latent heat flux, we obtain "flat tops" or "hats" in the However, quite a lot of parameter interaction is observed, with two clusters of parameters evident (i.e., see minimum stomatal resistance, potential evaporation soil moisture, diurnal heating penetration, residual soil moisture, etc.).
Despite the relatively low objective function values produced by the noninferior parameter range the corresponding total range of LSM flux prediction remains high. The eight Metflux site average ranges of prediction obtained from the Pareto parameter sets are shown in Figure 5. There is still a large amount of prediction variability even within the Pareto set, and hence the final parameter set must be chosen carefully.
All of the noninferior parameter sets are good but with different compromises as to which objective is more fully minimized. To select a single "best" parameter set from the 125 noninferior solutions, a necessarily subjective procedure was followed. For each of the three objectives, the five noninferior parameter sets giving the best value of that objective function were selected (a total of 15 solutions); in addition, the five members of the noninferior region giving the best value of hourly mean square error in matching of sensible heat flux also were selected as shown in

Conclusions
High-quality land surface water and energy balance simulations are essential for the prediction of drought, floods, agricultural production, land surface inputs to the atmosphere, and the response of land surfaces to climate change. However, the land surface water and energy balance models that make these predictions generally utilize large numbers of parameters, many of which are not regularly measured and some of which  Table 2. performance. Therefore, this study has explored the feasibility of using spatially distributed soils and vegetation information and a multiobjective parameter calibration method to improve the predictive ability of a spatially distributed LSM at the WGEW.
Although it seems reasonable that the inclusion of spatially variable soil and vegetation information in a spatially distributed LSM should improve local predictions, this study shows that the polygon nature of these data sets results in unrealistic simulations which were inconsistent with observations. A parameter sensitivity study revealed that spatial variations in vegetation parameters have little effect on predictions, but that soil parameters have a large effect. This problem is both an artifact of discretely assigning a single set such manipulations, and we seriously doubt that such manipulation would change our conclusions given the very strong spatial patterns observed in the subsequent predictions.
There is doubt that the process of deriving the spatially variable parameter fields from a few scattered observations using established physical relationships and table lookup is sufficient for realistically constraining a complex distributed land surface model. The spatially variable parameter fields were derived from a variety of sources, which leads to inconsistencies in the resulting parameter set and degraded LSM performance. Ideally, the model should be calibrated for each point in the domain to produce consistent parameter sets that account for the spatial variability of surface characteristics. For the WGEW this is a poorly posed problem, being that there are 90,000 model grid points and only eight surface flux observation points. Given that the available spatial information is insufficient to accurately specify the LSM spatial parameter distributions, we chose to calibrate the catchment average, parameters to provide the best possible simulation. Because precipitation is the dominant driver of spatial variability in soil moisture at the WGEW, this approach results in reasonable LSM spatial predictions. It is acknowledged that the simulation should improve if there were sufficient information to truly calibrate the LSM spatially. It may have been possible to use the spatially distributed PBMR soil moisture remote-sensing estimates to calibrate the spatially distributed parameters. However, this was not done because our intention in this study was to perform a multiobjective calibration that would improve the prediction of a number of quantities. If the LSM had been Soil Heat Flux calibrated spatially using only remotely sensed soil moisture, were derived by optimization. A multiobjective LSM then the surface energy flux prediction skill may have calibration technique was used to calibrate these parameters by degraded. estimating the noninferior parameter region based on a large The use of catchment-averaged parameters in this study is number of Monte-Carlo parameter set realizations, then also supported by the small difference between the calculated selecting the "best" parameter set in a semisubjective way. surface energy fluxes given by a single LSM with aggregate The chosen parameter set contained physically reasonable parameters and that given by a distributed array of LSMs. values which, in the case of the soil water retention parameters, This finding also suggests that a spatially distributed LSM applied at the WGEW with the resolution used in this study is probably not needed if catchment-averaged surface fluxes are required. However, a spatially distributed LSM is required for prediction of spatial patterns of soil moisture at the resolution of the available observations. LSMs have many parameters which must be carefully specified for their optimal implementation. Several parameters were directly observable at Walnut Gulch, and their values could be readily specified, but others were nonphysical and were superior to measured values (because of the removal of the coarse fraction prior to measurement).
Several additional insights into multiobjective calibration were gained as a result of additional experiments not reported in this paper and for which additional investigation is required.
These include the following: (1) Some parameters tend to have similar values for minimum objectives from several sites, indicating that watershed average parameters may be acceptable. However, the spatial patterns present in the PBMR data suggest that improved distributed representations may be possible, and therefore a lumped representation may not be acceptable for all the parameters. (2) Pareto optimal parameter sets often produce objective functions with large error indicating that some model components may require improvement (the soil heat flux component for example).
(3) Some objectives that are linked through model physics, such as