Global parameters sensitivity analysis of modeling water, energy and carbon exchange of an arid agricultural ecosystem

Agricultural and Forest Meteorology 271 (2019) 295–306

Contents lists available at ScienceDirect

Agricultural and Forest Meteorology journal homepage: www.elsevier.com/locate/agrformet

Global parameters sensitivity analysis of modeling water, energy and carbon exchange of an arid agricultural ecosystem ⁎

Mousong Wua,b, Youhua Ranc, , Per-Erik Janssond, Peng Chene,f,g, Xiao Tanh, Wenxin Zhangb,

T

⁎

a

International Institute for Earth System Science (ESSI), School of Geography and Ocean Science, Nanjing University, Nanjing, 210023, China Department of Physical Geography and Ecosystem Science, Lund University, Lund, SE-22362, Sweden c Key Laboratory of Remote Sensing of Gansu Province, Heihe Remote Sensing Experimental Research Station, Cold and Arid Regions Environmental and Engineering Research Institute, Chinese Academy of Sciences, Lanzhou, China d Department of Sustainable Development, Environmental Science and Engineering, KTH Royal Institute of Technology, Stockholm, Sweden e State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, China f Joint International Research Laboratory of Global Change and Water Cycle, Nanjing, China g College of Hydrology and Water Resources, Hohai University, Nanjing, China h State Key Laboratory of Hydraulics and Mountain River Engineering, College of Water Resource & Hydropower, Sichuan University, 610065 Chengdu, Sichuan, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Arid agricultural ecosystem Parameter sensitivity index Equifinality Water-carbon coupling Remote sensing

Agricultural ecosystems are important for regulating terrestrial hydrological and carbon cycles. Hydrological and carbon processes in agricultural ecosystem models are complex due to interactions between parameters. It is therefore crucial to identify parameter sensitivity before a process-based model is applied for simulations and predictions of water, energy and carbon fluxes in agricultural ecosystems. In this study, we investigated the sensitivity and equifinality of the CoupModel parameters in modeling an arid agricultural ecosystem in northwestern China. In total, 27 model parameters were analyzed using a global parameters sensitivity analysis approach and a combination of multiple in situ and remotely sensed data sets. Among the five major model processes, we found that the energy balance process account for much of the importance in the model, followed by soil hydrology, plant growth, soil heat, and soil carbon processes. Meanwhile, parameters from the plant growth process exhibited higher equifinalities than other processes. We found that net ecosystem exchange (NEE) is controlled by soil heat, soil hydrology and energy balance processes, which is mainly due to a high equifinality (0.91) between the parameters gmax (maximal stomatal conductance) and Vcmax (maximal carboxylation rate). The equifinalities between different parameters result in a trade-off in model performance metrics (i.e. determination coefficient R2 and mean error ME) in the water, energy and carbon balance simulations. We revealed that daytime and yearly accumulated eddy fluxes (sensible heat Hs, latent heat LE and NEE) can constrain the model parameters better. Remotely sensed data were also promising as additional constraints on soil water contents and energy fluxes. This study introduced a systematic global parameter sensitivity analysis approach together with the equifinality identification in an ecosystem model. The approach proposed here is applicable to other studies and the equifinalities detected in this study can be important implications for modelling arid agricultural ecosystems. Additional exploration on remotely sensed data in constraining the model from different aspects are highly recommended in modeling agricultural ecosystems.

1. Introduction Agricultural ecosystems play a vital role in controlling the hydrological and biophysical cycles on Earth. The agricultural ecosystem occupies approximated 12% of the global land surface, and most of them is distributed at middle to high latitudes (e.g. Europe and US) in areas with an arid climate (Wood et al., 2000; Ramankutty et al., 2008). These arid agricultural ecosystems have a stronger carbon uptake than

⁎

natural vegetation because of their high productivity (MacBean and Peylin, 2014). In last several decades, climate change has shown a threat to agricultural productivity in China due to changing water resources (Piao et al., 2010). Therefore, a process-based understanding in arid agricultural ecosystems is essential in addressing the issues such as food security and groundwater depletion in a changing climate. A better understanding of the hydrological and biophysical cycles in agricultural ecosystems resort to not only field experiments but also

Corresponding authors. E-mail addresses: [email protected] (Y. Ran), [email protected] (W. Zhang).

https://doi.org/10.1016/j.agrformet.2019.03.007 Received 12 August 2018; Received in revised form 2 February 2019; Accepted 11 March 2019 0168-1923/ © 2019 Elsevier B.V. All rights reserved.

Agricultural and Forest Meteorology 271 (2019) 295–306

M. Wu, et al.

Fig. 1. Study site location and experimental field setup. The bottom left panel is the location of the study region in China. The top left panel is a map of the Heihe River Basin and the location of the Yingke Irrigation District. The top right panel is the land use types in the Yingke Irrigation District (with a scale) and the irrigation canals are depicted. The bottom right panel depicts the pictures of maize and an EC tower.

uncertainty in modeling. Usually, a well-designed model calibration approach is necessary and it is achieved by using the Monte Carlo sampling method combining multiple model performance metrics or likelihood functions. The behavioral simulations are then selected using acceptance criteria with respect to different observations. During model calibration, a more efficient way to reduce equifinality of parameters requires measurements from different aspects of ecosystems (Beven and Freer, 2001). Wu et al. (2011a) used soil moisture and soil temperature to calibrate CoupModel in bare agricultural soil with seasonal frost conditions; and found that there was a trade-off in modeling soil temperature and soil water content in frozen soils. Senapati et al. (2016) used data from both soil profile observations and eddy covariance measurements to investigate the modeling uncertainty in a grassland area. They also noticed the influences of different observations on model performance in a multi-objective calibration. Metzger et al. (2016) used eddy covariance data and soil profile observations to study uncertainty and equifinality when modeling a peatland ecosystem; and addressed the importance of multiple observations in constraining model performance. Recently, remotely sensed data from different sources have also been simultaneously used in ecosystem modeling at multiple sites (Wu et al., 2019). In present study, we performed a global parameter sensitivity analysis using multiple sources of observations and multiple criteria to constrain water, energy and carbon fluxes in an arid agricultural ecosystem of northwestern China. Our objectives are: 1) to identify the parameter sensitivity with respect to various constraints on model outputs; 2) to detect interactions between parameters and how they result in equifinality in modeling the water, energy and carbon fluxes; and 3) to discuss the usefulness of various sources of observations in constraining the model processes from different aspects.

process-based models. On one hand, different sources of observations representing various components of agricultural ecosystems can be used to improve understanding of mechanisms underlying various processes; on the other hand, process-based models are needed to transfer the knowledge gained from measurements to various sites with similar conditions. Therefore, a combination of field experiments with processbased models can provide insight into the processes underlying the observations and lead to finding a better management strategy under a changing climate in the future. Over the last few decades, many process-based soil-plant-atmosphere continuum models have been developed, which combine observations to understand the ecosystems from different perspectives. Examples of these models include CENTURY (Parton, 1996), LPJGUESS(Smith et al., 2001), ORCHIDEE (Krinner et al., 2005), CoupModel (Jansson, 2012), Community Land Model (Bilionis et al., 2015), and MIKE SHE (Larsen et al., 2016). The CoupModel has been used to study different ecosystems including grassland (Senapati et al., 2016), forest (Wu et al., 2011b), and cropland (Nylinder et al., 2011). This model has shown the flexibility in modeling agricultural ecosystems: 1) as a process-based model, CoupModel can simulate agricultural ecosystems while considering detailed management, and 2) CoupModel allows for calibration based on the Monte Carlo sampling method, which is important before applying a model to an un-investigated site. Many parameter estimation methods have been developed, such as the simulated annealing algorithm, genetic algorithm, and Bayesian approaches (Li et al., 2004; Khu and Madsen, 2005; Su et al., 2009; Doherty and Hunt, 2010; Zhu et al., 2011). For process-based models with a large number of parameters, some methods may not be suitable due to extensive computational needs. The Monte Carlo-based approach was thus developed and is suitable for the estimation and calibration of a nonlinear model (Beven and Binley, 1992). Instead of searching for a single optimal parameter set, the Monte Carlo-based approach investigates an ensemble of different parameter sets that achieve equally good simulations, called ‘equifinality’ (Beven, 2006). The equifinality results from the interactions between different model parameters causes 296

Agricultural and Forest Meteorology 271 (2019) 295–306

M. Wu, et al.

2. Materials and methods

Table 2 Irrigation and fertilization management for maize at the study site.

2.1. Experimental site

Irrigation ʹ

ʹ

The Yingke Irrigation District (100°25 E, 38°51 N) is located in the middle reaches of the Heihe River Basin in northwest China. This site is a typical irrigated oasis cropland region with maize as the major crop (Fig. 1). The climate at the study site is temperate, with annual mean air temperature of 7.6 °C, annual mean precipitation of 117 mm, and annual mean potential evaporation of 2390 mm from 2008 to 2015. This site is one experimental area for the Watershed Allied Telemetry Experimental Research (WATER) project (Li et al., 2009). Biogeophysical and biogeochemical parameters, as well as meteorological and hydrological data, were collected via the Cold and Arid Regions Science Data Center at Lanzhou (http://westdc.westgis.ac.cn), which is a regular member of the ICSU World Data System (https:// www.icsu-wds.org/).

Date 18-May 15-Jun 16-Jul 15-Aug 8-Sep

Fertilization Irrigation amount 150 mm 150 mm 180 mm 180 mm 225 mm

Date 5-Apr 16-May 15-Jun 14-Aug /

Nitrogen amount 6.46 g/m2 7.02 g/m2 12.6 g/m2 18.4 g/m2 /

(HMP45AC, Vaisala), global radiation (CNR 4, Kipp&Zonen), precipitation (TE525 MM, Campbell), and wind speed at 2 m height (010C, MetOne). These data were aggregated to hourly values to drive the model. 2.5. Soil profile water and heat measurements Soil temperature and soil water content at six depths (10, 20, 40, 80, 120, 160 cm) were obtained by temperature (107, Campbell) and water content (CS616, Campbell) sensors from 2008 to 2015. Two heat flux plates (HFP01SC, Campbell) were placed at 5 cm and 15 cm depths to measure soil heat fluxes (2008–2015).

2.2. Soil data Soil texture at the study site was investigated in seven layers from the surface to 100 cm depth, as shown in Table 1. Soil texture information, including the clay-silt-sand fractions, was determined by a laser particle size analyzer (Microtrac S3500). Porosity and bulk density were determined by the cutting ring method, which takes soil samples from the field by using the standard cutting ring with a known volume (diameter of 5 cm and height of 5 cm) and dries them. Saturated hydraulic conductivity was determined in the laboratory by conducting fixed pressured head infiltration experiments. Soil total organic carbon (TOC), total carbon (TC), and inorganic carbon (IC) were then determined using an organic matter analyzer (TOC-VCPH total organic carbon analyzer).

2.6. Eddy covariance measurements An eddy covariance system (EC) was also installed at the study site for long-term use at a height of 2.85 m. Canopy CO2 and H2O fluxes were continuously measured from 2008 to 2015. Additional processing of data included quality check and the calculations of hourly sensible heat (Hs), latent heat (LE) and net ecosystem exchange (NEE). Meanwhile, to investigate the usefulness of data sampled from different periods of the day, flux data were divided into daytime (09:30–15:30 CET) and nighttime (22:30–02:30 CET) values representing gross primary productivity and ecosystem respiration, respectively (Metzger et al., 2016). Hourly Hs, LE and NEE fluxes were also aggregated to yearly accumulation values.

2.3. Crop management Maize was hole seeded on April 20th and harvested on September 22nd every year. The row spacing for maize was 55 cm, and the plant spacing was 22 cm. The height of maize measured at different growth periods was used as an input for the model. The maximum height reached 2.5 m. During the growing season (spanning from May to September), five irrigation events occurred (Table 2). The irrigation amount was obtained from the local water service and an agrometeorological experimental station. Four fertilization events were also applied during the growing season, according to data provided by the local water service (Wang et al., 2013). The nitrogen amount in the fertilization is shown in Table 2.

2.7. Remotely sensed data Remotely sensed surface soil moisture data extracted from the reprocessed Soil Moisture and Ocean Salinity (SMOS) L3 (ver. 300) daily soil moisture product (L-band, 1.4 GHz) with a spatial resolution of ˜25 km (2010–2015) from CATDS (www.catds.fr) was also used for calibration. The SMOS soil moisture product represents moisture within a thin surface layer (4 cm) at ˜25 km spatial resolution. The data were filtered by considering radio frequency interference (RFI) according to Scholze et al. (2016). Similar filtering was applied at grid points with a high probability of snow, frozen soil, water bodies, strong topography or urban areas. We also used the fraction of absorbed photosynthetically active radiation (FAPAR) data (16-day resolution, ˜1 km spatial resolution) derived from the Joint Research Center with a two-stream inversion

2.4. Meteorological data Meteorological data were collected from an automated weather station located at the Yingke experimental station from 2008 to 2015. The station provided half-hourly measurements including air temperature at 2 m height (HMP45AC, Vaisala), relative humidity Table 1 Soil texture, hydraulic properties and soil carbon. Soil depth (cm)

Clay (%)

Silt (%)

Sand (%)

Porosity

Bulk density (g/cm3)

Ks (m/s)

TOC (%)

TC (%)

IC (%)

0–5 5–10 10–20 20–40 40–60 60–80 80–100

5.27 4.78 4.46 5.40 5.41 9.93 8.48

66.01 65.03 67.81 71.94 63.57 74.90 73.88

28.72 30.19 27.73 22.66 31.02 15.17 17.64

0.51 0.50 0.49 0.42 0.40 0.42 0.41

1.31 1.37 1.47 1.46 1.53 1.57 1.52

5.00E-06 1.46E-06 1.88E-06 4.86E-06 2.25E-04 3.57E-06 9.26E-07

1.12 1.41 / / 0.58 / /

2.21 2.53 / / 1.75 / /

1.09 1.12 / / 1.17 / /

TOC: total organic carbon; TC: total carbon; IC: inorganic carbon. 297

Agricultural and Forest Meteorology 271 (2019) 295–306

M. Wu, et al.

Fig. 2. CoupModel structure in modeling water, energy and carbon processes. The left panel is a scheme for soil hydrological processes considered in this study. The middle panel is the soil heat and energy balance within the soil-plant-atmosphere continuum. The right panel is the carbon cycle in the soil and plant.

package (JRC-TIP; Pinty et al., 2007). The JRC-TIP FAPAR product (2010–2015) was derived by the JRC-TIP algorithm from MODIS 0.01degree albedo data (Pinty et al., 2001). Uncertainties in the FAPAR data were derived from the MODIS input albedo uncertainty by error propagation (Voßbeck et al., 2010). For more detailed information on the uncertainty propagation, we refer to the related papers (Pinty et al., 2011). 2.8. CoupModel This study used CoupModel v5.4, the executable of which can be downloaded from (http://coupmodel.com). A detailed description of CoupModel can be found in Jansson (2012). CoupModel can be used for applications with the options for modules and formulations with different levels of complexity. In Appendix 1 in Supplementary material we describe the processes applied in this study. The model represented the soil-plant-atmosphere system with a detailed description of C and N, as well as water and heat transport in soil and plants (Fig. 2). The soil profile was divided into 28 layers with an incremental thickness from 5 cm to 1 m, resulting in a total soil depth of 10 m. The initial conditions for the soil carbon pools as well as the soil water and heat status were obtained with 100 repeated runs of the model with forcing data from 2008. This process made the system more insensitive to the initial values. The calibrated parameters and most relevant equations can be found in Table A1. The major model assumptions with respect to the application of the CoupModel to arid agricultural systems are presented in the following sections.

Fig. 3. The processes of global parameter sensitivity analysis.

(LGM) method (Lindeman et al., 1980) to identify important parameters and to investigate the equifinality between parameters based on 30,000 model realizations; 3) tests on the usefulness of additional observations except in situ data in reducing equifinality. We acknowledge that 30,000 parameter sets will not fully explore the entire parameter space for a 27-parameter model. However, our approach here serves as a demonstrative exercise to provide an example of how a modeler could perform a global sensitivity analysis via a process-based model and multiple sources of observations. As such, a greater number of parameter sets should not influence our general conclusions. We have compared the differences in model performance with 20,000 runs and

3. Model parameterization and sensitivity analysis 3.1. Model sensitivity analysis procedure The sensitivity analysis consisted of three steps (Fig. 3): 1) a precalibration procedure with the selection of potentially important parameters using Morris screening method (Morris, 1991); 2) parameter sensitivity identification using Lindeman, Gold and Merenda 298

Agricultural and Forest Meteorology 271 (2019) 295–306

M. Wu, et al.

30,000 runs and did not find obvious changes.

3.5. Selection of behavioral simulations For the basic and extra selections, different criteria were used for the selection of behavioral simulations. In the basic selection, two performance metrics were applied: coefficient of determination (R2) to assess the fit of the model to the measurements with respect to dynamics and mean error (ME) to minimize the bias between model and measurements. For soil temperature observations (B1 in Table 3, hereafter the same), R2 ≥ 0.9 and |ME |≤ 1.2 °C were applied to soil temperature at 10, 20 and 40 cm depths (T1–T3); for soil water observations (B2), R2 ≥ 0.5 were used to constrain soil water content at 10, 20, and 40 cm depths (θ1–θ3), but |ME |≤ 4.5%, |ME| ≤ 1.4%, and |ME| ≤ 6.0% cm3 cm−3 were applied to each depth; for soil heat flux observations (B3), soil heat flux at 5 cm depth (q1) was constrained by selecting simulations with R2 ≥ 0.33 and |ME| ≤ 80,000 J m−2 d−1, and soil heat flux at 15 cm depth (q2) with R2 ≥ 0.18 and |ME| ≤ 40,000 J m−2 d−1; for eddy flux observations (B4), sensible heat (Hs) was selected using R2 ≥ 0.28 J m−2 d−1, soil latent heat (LE) was selected by R2 ≥ 0.13 J m−2 d−1, and net ecosystem exchange (NEE) was selected by R2 ≥ 0.18 gC m−2 d−1. For each group of observations, the rejection efficiency, i.e. the ratio of rejected runs to total runs, was calculated by constraining each observation variable separately or all variables together. The basic selection was performed step-by-step, from constraining soil temperature as the first step, then soil water, soil heat fluxes and eddy fluxes. For the extra selection, we focused on only the performance metrics R2. We constrained R2 for achieving 2% of the best as behavioral simulations (i.e. rejection efficiency was controlled to 98%) based on 30,000 runs. For daytime and nighttime net ecosystem exchange (NEE), in addition to fitting the simulations to NEE, we also fitted them to total photosynthesis and total respiration. Net radiation Rn was compared with modeled net radiation. We fitted SMOS soil moisture to water storage at the surface layer from 0 to 4 cm depth. For FAPAR, there was no corresponding model variable in the output, but we knew that it is related to radiation absorbed by photosynthesis (i.e. PAR); therefore, we fitted FAPAR to photosynthetically active radiation to check the performance in dynamics.

3.2. Morris parameter screening The Morris method is an efficient screening approach for identifying the few potentially important parameters among a large number of parameters from a physically-based model. In this study, 131 parameters were thought to have effects on five processes (soil carbon, plant growth, soil hydrology, energy balance and soil heat) in the CoupModel. It would be computationally expensive to calibrate all parameters, as most of them do not have large impacts on model output variables of interest. We thus used the Morris screening method to rank the parameter importance for the water and energy balance and selected the most important parameters for calibration. This was performed by varying the value of each parameter randomly within its range and calculating the elementary effect as

Ei =

f (p1 , p2 , …, pi + Δ, …, pk ) − f (p1 , p2 , …, pi , …, pk ) Δ

(1)

where Δ is the step size of a change in parameter pi , and function f represents the model output for which the sensitivity is studied, such as multiyear water balance, energy balance and carbon balance. For estimation of μi and σi for the elementary effect, we repeated the calculation of the elementary effect several times with random starting points for the parameters (20 times here). Then we separated the important parameters from the unimportant parameters using

βi =

μi2 + σi2 σ

(2)

where βi is a measure of importance for parameter pi , and σ is the standard deviation of μi2 + σi2 for all parameters. A detailed description of how to perform parameter screening using the Morris method is presented in Appendix 2 in Supplementary material. Finally, 27 parameters related to different processes in CoupModel were selected as potentially important parameters, as shown in Table A1. For the rest of the parameters, fixed values were assigned according to model default values or literature-suggested values.

3.6. LGM relative importance

3.3. Potentially important parameters

As a way of quantifying the sensitivity of parameters to certain model outputs in the calibration, the relative importance of each calibrated parameter was calculated based on the standardized covariance matrix of parameters and related model performance metrics using the LGM method (Lindeman et al., 1980). The LGM method averages the sequential sums of squares over all orders of regressors and can obtain the relative importance of each parameter in impacting model results (see Appendix 3 in Supplementary material). The LGM method is a compensation for the Morris method, because the Morris method accounts for the influences of only a single parameter. The LGM relative importance can be formulated as

The potentially important parameters are related to formulas describing five processes: 1) soil carbon ( pθSatact , tmax , tmin , fe, l ), 2) plant growth (tpQ10 , θAmin , Emerge Tth, pmin , gmax , Vcmax ), 3) soil hydrology (SurfPoolInit, ψa (4) , ψa (5) , ψa (7) , ψa (8) , InitialPressuredHead), 4) energy balance (nc , a1, a2 , rk1, rk2 , TSnowL , sk ), and 5) soil heat (Cmd , b1, qh, low ). Prior ranges shown were selected according to the literature or model default setting. Uniform (or logarithm uniform) prior distributions were assumed for the parameters. 3.4. Measurement variables The variables of interest (Table 3) were separated into two groups: 1) basic selection and 2) extra selection. For the basic selection, soil temperature (T) at 10, 20 and 40 cm, soil water content (θ) at 10, 20 and 40 cm, soil heat fluxes (q) at 5 and 15 cm, sensible heat (Hs), latent heat (LE) and net ecosystem exchange (NEE) from 2008 to 2015 were chosen. For extra selection, daytime and nighttime eddy fluxes (Hs, LE and NEE) as well as the yearly accumulated values were used to evaluate the usefulness of different types of datasets in constraining daytime, nighttime, and yearly accumulated fluxes. Additionally, net radiation Rn, SMOS soil moisture and FAPAR from 2010 to 2015 were used in extra selection to constrain modeled net radiation, surface soil moisture and PAR (photosynthetically active radiation), respectively.

LGM (xk ) =

⎞ ⎛ p−1 seqR2 ({xk }|S ) ⎟ 1 ⎜ ∑ ∑ p−1 ⎟ p j = 0 ⎜ S ⊆ {x1, … , xp } ∖ {xk } i ⎟ ⎜ = n S j ( ) ⎠ ⎝

( )

(3)

where seqR2 ({xk }|S ) is the additional R2 when adding the regressor {xk } to a regression model with regressors S . We calculated the LGM relative importance based on a multivariate regression between the parameters and coefficient of determination R2 for different model output variables (11 variables shown in Table 3 for basic selection). 299

Agricultural and Forest Meteorology 271 (2019) 295–306

M. Wu, et al.

Table 3 Calibration variables and selected criteria. Category Basic selection Soil temperature B1

T1+T2+T3 Soil water B2

θ1+θ2+θ3 B1+B2 Soil heat flux B3

Measured variable

Unit

Observations

R2

ME

Rejection efficiency (%)

T1 T2 T3

o

C C o C 95

62153 62154 62154

≥0.9 ≥0.9 ≥0.9

± 1.2 ± 1.2 ± 1.2

95 95 87

θ1 θ2 θ3

% % % 80

62060 62154 62154

≥0.5 ≥0.5 ≥0.5

± 4.5 ± 1.4 ± 6.0

73 44 67

J m−2 d−1 J m−2 d−1 32

62155 33282

≥0.33 ≥0.18

± 80000 ± 40000

18 32

J m−2 d−1 J m−2 d−1 gC m−2 d−1 57

44530 43699 46453

≥0.28 ≥0.13 ≥0.18

– – –

47 45 54

J m−2 d−1 J m−2 d−1 gC m−2 d−1 gC m−2 d−1

15997 15708 15845 15845

≥0.57 ≥0.87 ≥0.29 ≥0.29

– – –

98 98 98 98

J m−2 d−1 J m−2 d−1 gC m−2 d−1 gC m−2 d−1

7325 7198 8078 8078

≥0.02 ≥0.09 ≥0.01 ≥0.01

– – – –

98 98 98 98

J m−2 d−1 J m−2 d−1 gC m−2 d−1

44530 43699 46453

≥0.63 ≥0.44 ≥0.31

– – –

98 98 98

J m−2 d−1

70888

≥0.91

–

98

mm –

66 61

≥0.12 ≥0.32

– –

98 98

o

99 q1 q2

q1+q2 B1+B2+B3 Eddy flux B4

99 Hs LE NEE

Hs+LE + NEE B1+B2+B3+B4 Extra selection Daytime eddy flux E1

99.7

Hs (daytime) LE (daytime) NEE (daytime) NEE (daytime fitting photosynthesis)

Nighttime eddy flux E2

Hs (nighttime) LE (nighttime) NEE (nighttime) NEE (nighttime fitting respiration) Yearly accumulated eddy flux E3 Hs (yearly accumulation) LE (yearly accumulation) NEE (yearly accumulation) Meteorological data E4 Rn Remote sensing E5 SMOS_sm FAPAR

−

T1–T3: Soil temperature at 10, 20, 40 cm; θ1–θ3: Soil water content at 10, 20, 40 cm; q1–q2: Soil heat flux at 5, 15 cm; Hs: Sensible heat flux; LE: Latent heat flux; NEE: Net ecosystem exchange; Rn: Net radiation; SMOS_sm: Soil moisture at surface 4 cm depth derived from SMOS (Soil Moisture and Ocean Salinity) satellite; FAPAR: Fraction of Absorbed Photosynthetically Active Radiation.

of the parameters from the soil carbon process was 0.6%, from the plant growth process was 22%, from the soil hydrology process was 27%, from the energy balance process was 34%, and from the soil heat process was 17%. For soil temperatures at 10, 20 and 40 cm, the parameters controlling the energy balance process had the greatest importance (approximately 66% to 71%) to soil temperature R2, followed by the parameters related to the soil heat process (relative importance approximately 24% to 29%). For soil moisture at 10 and 20 cm depths, the parameters from the soil hydrology process group showed the greatest importance (68% to 73%). For soil moisture R2 at 40 cm depth, the parameters related to the energy balance had the greatest importance (52%), followed by the parameters related to soil heat (42%); soil hydrology related parameters only showed 6% relative importance. The soil heat flux (q1, q2) R2 was highly influenced by soil hydrology, with parameters in the soil hydrology process showing 54% to 60% relative importance. For eddy fluxes R2, namely, Hs, LE and NEE, parameters belonging to the plant growth group had the greatest importance (approximately 74% to 94%).

3.7. Equifinality estimates To quantify the interactions between parameters, the equifinality for each parameter was defined based on the correlation coefficients of the interacting parameter pairs in the accepted runs (Metzger et al., 2016): 2

Eqi =

∑2 × j

210 × Rij, i ≠ j 10

(4)

where Ri, j is the correlation coefficient between parameters i and j, calculated from the behavioral parameter sets. 4. Results 4.1. Model process relative importance The relative importance of the parameters to the model performance R2 for different variables was categorized into five groups (i.e., soil carbon, plant growth, soil hydrology, energy balance and soil heat) (Fig. 4). For each variable, the relative importance of all parameters was normalized to total 100%. In general, the mean relative importance 300

Agricultural and Forest Meteorology 271 (2019) 295–306

M. Wu, et al.

Fig. 4. LGM relative importance of parameters from different groups of model processes on model performance metrics R2 from different model output variables. T1-T3 denotes soil temperature at 10, 20, 40 cm depths, θ1-θ3 denotes soil moisture at 10, 20, 40 cm depths, q1-q2 denotes soil heat flux at 5, 15 cm depths, Hs is sensible heat fluxe, LE is latent heat flux, NEE is net ecosystem exchange.

Fig. 5. LGM relative importance of most sensitive parameters in each parameter group. Parameters are defined in Table A1.

301

Agricultural and Forest Meteorology 271 (2019) 295–306

M. Wu, et al.

coefficient for snow), with a PSI of 0.70. Regarding the soil heat parameter group, two parameters that showed quite high sensitivity were Cmd and qh, low , with PSI values of −0.37 and 0.19, respectively.

4.2. Most important parameters The most important parameters in each parameter group are shown in Fig. 5. Detailed relative importance for each parameter is shown in Fig. A1. In the energy balance calculation, rk1 showed relative importance of 37.98%–40.83% to the soil temperature (Fig. 5a–c). Another parameter Cmd had around 18.31%–26.03% relative importance to soil temperatures at 10 cm and 20 cm depths. For the deeper soil layer e.g., 40 cm depth, soil temperature was more sensitive to the lower boundary heat flux qh, low . At 10 cm and 20 cm depths (Fig. 5d–e), the air entry at the upper layer (ψa (4) , air entry in soil water retention curve) was the most sensitive parameter to soil water content, with relative importance of 20.73% and 54.94%, respectively. At 40 cm depth (Fig. 5f), the most sensitive parameters were Cmd and rk1, showing relative importance of 39.77% and 31.25% to soil water content, respectively. For soil heat fluxes (q1 and q2) at 5 cm and 15 cm depths (Fig. 5g–h), the most sensitive parameter was ψa (4) , with 33.52% and 55.89% of relative importance, respectively. Cmd had 18.30% and 21.74% of relative importance to soil heat fluxes at 5 cm and 15 cm depths, respectively; followed by rk1, with 16.29% and 7.03% of relative importance, respectively. For eddy covariance fluxes (Fig. 5i–k), the most sensitive parameters were shown from the plant growth group. For sensible heat (Hs) and latent heat (LE) fluxes, the most sensitive parameter was gmax (maximal conductance of fully open stomata), with relative importance of 89.60% and 56.74%, respectively. Our result is consistent with Lu et al. (2013), who found that the leaf stomatal conductance is most important to mean annual heat flux, as also noticed by Milly and Shmakin (2002). For net ecosystem exchange (NEE), the most important parameter was Vcmax (maximum Rubisco capacity per leaf area at top of the canopy), with 41.62% of relative importance. This is similar to previous studies on land surface models (Dai et al., 2004; Alton, 2011; Kattge et al., 2009; Wang and Jarvis, 1990), which have shown that Vcmax is very important for photosynthesis since most of these models use the Farquhar model (Farquhar et al., 1980) for leaf photosynthesis.

4.4. Equifinality The total equifinality for each parameter was categorized into different parameter groups (Fig. 7). The parameters from the soil carbon process group showed equifinality of 5.53 to 5.64, with fe, l having the highest equifinality of 5.64. High equifinality ranging from 5.46 to 6.43 existed in parameters from the plant growth process group, with the parameters gmax and Vcmax having the highest equifinalities of 6.27 and 6.43, respectively. In soil hydrology process group, high equifinality of 5.41 to 5.82 was detected between the parameters, and kmat (7) had the highest equifinality value of 5.82. Equifinality for the parameters inside the energy balance group had the value of 5.51 to 6.05. For the parameters from soil heat process group, high equifinality was also detected, with value of 5.61 to 5.92, and qh, low showed the highest value of 5.92. 5. Discussion 5.1. Parameter sensitivity The sensitivity of different model variables to various parameters indicated that the tight interactions between different model processes. In the arid agricultural ecosystem, we found that the energy balance process is tightly coupled to the soil profile water and heat transport, suggesting that the calibration of one process would influence the other processes. The sensitivity of the parameters also varied with variables of interest. As shown in Fig. 6, when we used the same criteria (0.5% acceptance) for selecting different variables, we got different posterior parameter ranges, and even non-overlapping ranges (such as rk2 and Cmd ). This is an important finding since a lot of previous work only investigated the parameter sensitivity of certain processes of interest or only processes from the same group of model output variable (e.g. Wania et al., 2010; Zhu and Zhuang, 2014). This would lead to results that are not robust. In this study, we proposed a systematic global sensitivity analysis approach which considered the interactions of parameters and their influences on model performance. This approach together with the PSI calculated in Section 4.3 would help modelers to develop or select an appropriate parameter set (including parameter ranges) in a model with proper processes considered. Besides, the global parameter sensitivity analysis approach is an indispensable step for model calibration and prediction (Saltelli et al., 2000), by selecting the most sensitive parameters with well-constrained ranges. It should be taken into accounts carefully in modeling ecosystems with interacted processes. Although we do not expect to transfer the same parameter set from one site to another based on our analysis, the interactions between parameters from different processes and their influences on model performance are supposed to be site independent. For example, in Section 4.2 we found the parameters gmax and Vcmax have large influences on model performance, these two parameters are also found in other studies as key parameters for NEE modeling, e.g., MAESTRO (Wang and Jarvis, 1990), CLM2L (Dai et al., 2004), JSBACH (Kattge et al., 2009), BEPS (Ju et al., 2010), CABLE (Lu et al., 2013). This indicated that the sensitivities of these parameters are quite independent of model structure, as well as study sites and scales. We suggested the ranges of these parameters from one site should be compared to other ecosystems and models to identify site-specific parameters before they used for prediction in a larger scale.

4.3. Parameter sensitivity index (PSI) To quantify the sensitivity of each parameter for calibration, a general parameter sensitivity index (PSI) was thus defined here. We selected behavioral simulations based on the constraint on 11 model output variables. We chose the best 0.5% of R2 values for each variable separately. Then we compared the overlapped ranges of each parameter from the 11 parameter sets and calculated the PSI as the ratio of overlapping parameter ranges to the prior parameter ranges. The overlapping parameter range was defined as the range between the minimum upper ranges and the maximum lower ranges of parameter values from all 11 parameter sets. A smaller PSI value indicates less overlap and thus a more sensitive parameter. If the PSI < 0, i.e., parameter ranges from 11 selections are not overlapping, this parameter is seen as a highly sensitive parameter. Parameters belonging to the soil carbon group showed a high overlapping range, with PSI values from 0.88 to 0.91 (Fig. 6). In plant growth group, two sensitive parameters, namely, gmax (maximum stomatal conductance), and Vcmax (maximum carboxylation rate), had PSIs of 0.07 and 0.29, respectively. In addition, the parameter θAmin (the minimum amount of air that is necessary to prevent any reduced uptake of water from the soil) had a PSI of 0. In the soil hydrology parameter group, the parameters related to soil water retention showed high sensitivity (e.g., air entry ψa at different soil depths). In the energy balance parameter group, two parameters controlling atmospheric emissivity calculations showed quite high sensitivity, namely, rk1 and rk2 , with PSIs of 0.05 and −0.09, respectively. Another sensitive parameter was sk (thermal conductivity

5.2. Equifinality In calibration of a complex model, the interactions between parameters are more important in influencing model performance than a particular parameter (e.g. Beven and Freer, 2001). This results in the 302

Agricultural and Forest Meteorology 271 (2019) 295–306

M. Wu, et al.

Fig. 6. Parameter range and parameter sensitivity index (PSI). All parameters are normalized to 0 to 1, and different filled patterns indicate the overlapping parameter ranges from 11 parameter sets for different model processes. The PSI for each parameter is shown with a red star (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Fig. 7. Total equifinalities of parameters categorized into different parameter groups. For each parameter, the total equifinality is the sum of its equifinality with parameters from different groups, as indicated by different color bars.

Fig. 8. Counts of parameters from different process groups with strong correlation (|r| > 0.2) to R2 and ME for different variables.

We also noticed that the equifinalties existed not only between two parameters, but also among several parameters (e.g. ψa (5) − kmat (7) − nc − a2 as shown in Fig. A2). These equifinalities needed further investigation since the current simple regression method is not sufficient to describe interactions between several parameters. In the basic selection, we used 89 behavioral runs to obtain the equifinalities, and the definition of constraint criteria was subjective. This might not be enough to do more analysis on the equifinalities between more than two parameters in such a complex model as CoupModel. Nevertheless, the threshold applied in basic selection was high enough (more than 98% rejection efficiency) in identifying strong relationships between two parameters from various processes. High equifinalities with values larger than 0.26 (corresponded to correlation value of 0.2) were still useful in future calibration design for better understanding

trade-off in constraining certain parameters. Some previous studies have shown that in carbon modeling, the parameter sensitivities are influenced by other parameters (e.g. Verbeeck et al., 2006; Quillet et al., 2013). Therefore, when only parts of the parameters and processes are calibrated, an un-comparable or un-transferable parameter set would be obtained. As we have shown in Section 4.4, equifinalities existed not only between parameters inside the same model process, but also between different processes. This indicated the calibration of only a few parameters or certain process would result in problems when transferring parameter sets to another model, because models can differ in structures and parameterizations (Wania et al., 2010; Sándor et al., 2016). We thus recommend that better knowledge on equifinalities is necessary for better constraining parameters in model calibration in the future. 303

Agricultural and Forest Meteorology 271 (2019) 295–306

M. Wu, et al.

Fig. 9. Ratio of selected posterior parameter range to prior range for extra selection a) E3 and b) E4 + E5. The extra selection is illustrated in Table 3.

data from soil water, energy and carbon are all necessary, and they cannot be fully replaced by another (Metzger et al., 2016). For example, in B3 we did not get obvious changes in rejection efficiency when soil heat flux data were used as constraints, however, we still found they are useful because they showed different sensitivities to parameters as compared to soil temperatures (Fig. 5). When daytime eddy fluxes were used in calibration (E1), most parameter ranges did not show large differences. The performance metrics R2 was much higher when using the daytime eddy fluxes than using both the daytime and the nighttime data in calibration. The use of the nighttime fluxes (E2) did not change parameter ranges too much, and the R2 was also very low. This might be due to high uncertainty in deriving the nighttime fluxes by instruments. When the yearly accumulated fluxes were used for calibration in E3, some parameter ranges showed obvious differences, and the model performance was improved significantly (Fig. 9a). Use of the yearly accumulated Hs in calibration reduced the range of gmax to 0.57 of its prior range. The yearly accumulated LE also reduced ψa (5) to 0.90 of its prior range, and reduced ψa (7) to 0.64 of its prior range. Meanwhile, use of the yearly accumulated NEE in calibration reduced ψa (5) to 0.82 of the prior range, and reduced ψa (7) to 0.80 of its prior range. This implied the potential of using the accumulated flux data for better constraining model in a long term, as also suggested by Metzger et al. (2016), who combined the daytime and the accumulated eddy fluxes in constraining the model. When Rn, SMOS soil moisture and FAPAR were used to constrain model variables, three parameters ( gmax , a1, qh, low ) showed a slight reduction in posterior ranges compared to their prior ranges (Fig. 9b). Since Rn already showed quite good performance based on the energy balance scheme, the use of Rn in calibration did not significantly improve model performance. In another modeling work with ORCHIDEE, the use of Rn was found redundant if NEE and LE data were available (Santaren et al., 2007). However, Metzger et al. (2016) demonstrated the usefulness of Rn in constraining the modeling of peatland processes. They explained that variables which might have been identified as redundant in one study could be of high importance on another study when a different ecosystem or a different model is applied. As we have already addressed the importance of energy balance on model performance in both soil heat and hydrology, we thus recommend using other observations such as albedo in constraining energy balance process. SMOS soil moisture had only limited records due to the filtering out of data affected by snow cover and soil freezing/thawing. For the study site, with limited SMOS data we showed slight improvement on surface soil moisture. The finding was already promising because we have demonstrated that the soil hydrological process is important for modeling ecosystem carbon fluxes, and the use of additional observations for constraining soil water and energy balance will be beneficial to the ecosystem carbon cycle. Remotely sensed surface soil moisture has shown to be effective in constraining the global carbon cycle (Scholze et al., 2016). Therefore, obtaining more remote sensing data in cold regions would be necessary in future study of hydrological and carbon cycles.

and description of equifinalities. 5.3. Parameter interactions and their influences on model performance As equifinalties exist between parameters from different processes, they will inevitably have influences on model performance. In Fig. 8, the counts of parameters from different process groups that show strong correlations with R2 and ME with respect to different variables of interest are depicted. If a parameter shows correlation coefficient to model performance metrics with absolute value larger than 0.2, it is assumed sensitive to the model performance. Generally, there are around 2 to 7 parameters that are sensitive for two model performance metrics among the 11 selected variables. For different variables, model performance metrics are not only sensitive to parameters from the most relevant processes, but also sensitive to parameters from other processes. This results in strong correlations between parameters and model performance metrics of various variables. High equifinalities between gmax , Vcmax , nc , a2 , rk1, rk2 resulted in strong correlations between rk2 and soil water contents R2 (0.28 to 0.49), and between ψa (4) and soil heat fluxes R2 (0.54 to 0.60). This was mainly due to the tight coupling in water-carbon-energy balance in this agricultural ecosystem. We had the highest equifinality between gmax and Vcmax (0.91). This indicated a strong coupling in water-carbon cycle in the agricultural ecosystem. Lu et al. (2013) found that leaf stomatal conductance gmax and Vcmax are more important for annual latent heat flux and annual gross annual primary productivity based on sensitivity analysis within a land surface model (CABLE, Krinner et al., 2005). Due to the fact that plant transpiration and photosynthesis are tightly coupled, parameters related to these two processes should be calibrated carefully. We also found the differences in two performance metrics (R2 and ME) with respect to their correlations to different parameters. These two metrics represent model performance from different aspects, i.e. R2 for dynamics and ME for mean of error between model and measurements. The different links between these two metrics and parameters indicated that in evaluating the model performance, we need to choose a generally proper model performance metrics/likelihood or combine some of them. We have done this here by combining R2 and ME to constrain the model, and this combination was shown efficient to find more reliable parameter ranges (Wu et al., 2011a; Metzger et al., 2016; Senapati et al., 2016; Wu et al., 2016). 5.4. Usefulness of measurement variables To investigate the usefulness of measurement variables in reducing modeling uncertainty, both results from basic selection and extra selection were analyzed. The rejection efficiency increased from 95% in B1 to 99% in B2, until 89 behavioral runs were obtained in B4, by rejecting 99.7% of the total runs (Table 3). This implied that for a complex model with strong interactions between different processes, 304

Agricultural and Forest Meteorology 271 (2019) 295–306

M. Wu, et al.

ecosystem models.

FAPAR is a variable controlling plant phenology. The inclusion of FAPAR as an additional constraint resulted in a high correlation with absorbed radiation. FAPAR has been shown to be pronounced in constraining ecosystem carbon cycles at different scales in data assimilation schemes (Kaminski et al., 2012; Wu et al., 2019). In the next step, it would be necessary to implement a scheme in CoupModel for calculating FAPAR and include FAPAR as a model output for better calibration of model phenological process.

Acknowledgements This study was jointly supported by National Key Research and Development Program of China (2016YFA0600204, 2017YFC0403304), National Natural Science Foundation of China projects (41471359), and the Youth Innovation Promotion Association of Chinese Academy of Sciences (Grant number: 2016375). Study was also funded by the Opening Research Foundation of the State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University (2017NSG02). We thank Michael Voßbeck for providing SMOS soil moisture original data and the JRC-TIP FAPAR data for the study site.

6. Conclusions The investigation of parameter sensitivity is important for improved modeling of the water, energy and carbon balance in a representative arid agricultural ecosystem. Therefore, we performed the global parameter sensitivity analysis in modeling an arid agricultural ecosystem in northwestern China using multiple observations and the process-based CoupModel. The major findings are as follows: We found that the surface energy balance is most sensitive in calibration and can affect soil hydrology and soil heat transport, as well as the ecosystem carbon cycle. The surface energy balance and the snow energy balance are sensitive to estimates of atmospheric emissivity. In soil profile, calculations of soil heat and soil water fluxes are very sensitive to parameters controlling soil thermal and hydraulic properties. While above the ground, plant transpiration and photosynthesis are tightly linked. Parameters showing high sensitivity to these processes should be taken into accounts carefully by assigning proper ranges in future modeling. In addition to the influence of a single parameter, equifinalities between parameters have strong impacts on model performance and result in trade-offs in different performance metrics. When modeling the surface energy balance, the latent heat (LE) and sensible heat (Hs) performance exhibit a trade-off due to high equifinalities between the parameters in calculating energy balance, which indicates a more accurate constraint on surface energy balance is necessary. When modeling the net ecosystem exchange (NEE), a trade-off was also detected between model performance in ecosystem energy balance and carbon cycle. This trade-off is mainly due to the strong coupling of watercarbon cycles in ecosystems, especially the high equifinality of 0.91 between g max (maximal stomatal conductance) and Vcmax (maximal carboxylation rate). Therefore, when modeling agricultural ecosystems, we need to address these couplings by using multiple observations and different model performance metrics to constrain ecosystem processes from various aspects. By constraining the model with different types of observations, we found that the daytime and yearly accumulated eddy fluxes can achieve higher model performance than nighttime data due to the high level of noise in nighttime flux data. These two types of eddy fluxes have the potential in constraining the agricultural ecosystem in a combined way. Remotely sensed data have the advantage of being easily accessed at various temporal and spatial resolutions. As discussed above, the coupling of the water, energy, and carbon cycles in agricultural ecosystems is strong, therefore, we suggest that remotely sensed data with high temporal and spatial resolutions used as additional observations in constraining ecosystem processes from more aspects in the future. Generally, we can conclude that in modeling of arid agricultural ecosystems, we need to identify parameter sensitivity and equifinality using a systematic global parameter sensitivity analysis approach as provided in this study. Since the ecosystem has shown strong couplings between water, energy and carbon cycles, the use of multiple observations is necessary in constraining parameters and processes from different aspects, together with a proper selection of the model performance metrics. Moreover, except for in situ observations, the remotely sensed data with high temporal and spatial resolutions are promising in reducing uncertainty in modeling agricultural ecosystem. We thus recommend more remotely sensed data such as albedo, soil moisture, and FAPAR for better parameter estimation in calibration of

Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.agrformet.2019.03. 007. References Alton, P.B., 2011. How useful are plant functional types in global simulations of the carbon, water, and energy cycles? J. Geophys. Res. Biogeosci. (G1), 116. Beven, K., 2006. A manifesto for the equifinality thesis. J. Hydrol. 320 (1-2), 18–36. Beven, K., Binley, A., 1992. The future of distributed models: model calibration and uncertainty prediction. Hydrol. Process. 6 (3), 279–298. Beven, K., Freer, J., 2001. Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the GLUE methodology. Journal of Hydrology 249 (1-4), 11–29. Bilionis, I., Drewniak, B.A., Constantinescu, E.M., 2015. Crop physiology calibration in the CLM. Geosci. Model. Dev. 8 (4), 1071–1083. Dai, Y., Dickinson, R.E., Wang, Y.-P., 2004. A two-big-leaf model for canopy temperature, photosynthesis, and stomatal conductance. J. Clim. 17 (12), 2281–2299. Doherty, J.E., Hunt, R.J., 2010. Approaches to Highly Parameterized Inversion: a Guide to Using PEST for Groundwater-model Calibration. US Department of the Interior, US Geological Survey. Farquhar, G.D., von Caemmerer, S., Berry, J.A., 1980. A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species. Planta 149 (1), 78–90. Jansson, P.-E., 2012. CoupModel: model use, calibration, and validation. Trans. Asabe 55 (4), 1337–1344. Ju, W., Wang, S., Yu, G., Zhou, Y., Wang, H., 2010. Modeling the impact of drought on canopy carbon and water fluxes for a subtropical evergreen coniferous plantation in southern China through parameter optimization using an ensemble Kalman filter. Biogeosciences 7 (3), 845–857. Kaminski, T., Knorr, W., Scholze, M., Gobron, N., Pinty, B., Giering, R., Mathieu, P.-P., 2012. Consistent assimilation of MERIS FAPAR and atmospheric CO2 into a terrestrial vegetation model and interactive mission benefit analysis. Biogeosciences 9 (8), 3173–3184. Kattge, J., Knorr, W., Raddatz, T., Wirth, C., 2009. Quantifying photosynthetic capacity and its relationship to leaf nitrogen content for global-scale terrestrial biosphere models. Glob. Chang. Biol. 15 (4), 976–991. Khu, S.T., Madsen, H., 2005. Multiobjective calibration with Pareto preference ordering: an application to rainfall-runoff model calibration. Water Resour. Res. 41 (3). Krinner, G., Viovy, N., de Noblet‐Ducoudré, N., Ogée, J., Polcher, J., Friedlingstein, P., Ciais, P., Sitch, S., Prentice, I.C., 2005. A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system. Global Biogeochem. Cycles 19 (1). Larsen, M.A., Refsgaard, J.C., Jensen, K.H., Butts, M.B., Stisen, S., Mollerup, M., 2016. Calibration of a distributed hydrology and land surface model using energy flux measurements. Agric. For. Meteorol. 217, 74–88. Li, X., Koike, T., Pathmathevan, M., 2004. A very fast simulated re-annealing (VFSA) approach for land data assimilation. Comput. Geosci. 30 (3), 239–248. Li, X., et al., 2009. Watershed allied telemetry experimental research. J. Geophys. Res. Atmos. (D22), 114. Lindeman, R., Merenda, P., Gold, R., 1980. Introduction to Bivariate and Multivariate Analysis. Foresman and Co., London. Lu, X., Wang, Y.-P., Ziehn, T., Dai, Y.J., 2013. An efficient method for global parameter sensitivity analysis and its applications to the Australian community land surface model (CABLE). Agric. For. Meteorol. 182, 292–303. MacBean, N., Peylin, P., 2014. Biogeochemistry: agriculture and the global carbon cycle. Nature 515 (7527), 351. Metzger, C., Nilsson, M.B., Peichl, M., Jansson, P.-E., 2016. Parameter interactions and sensitivity analysis for modelling carbon heat and water fluxes in a natural peatland, using CoupModel v5. Geosci. Model. Dev. 9 (12), 4313. Milly, P., Shmakin, A., 2002. Global modeling of land water and energy balances. Part II: land-characteristic contributions to spatial variability. J. Hydrometeorol. 3 (3), 301–310. Morris, M.D., 1991. Factorial sampling plans for preliminary computational experiments.

305

Agricultural and Forest Meteorology 271 (2019) 295–306

M. Wu, et al.

European climate space. Glob. Ecol. Biogeogr. 10 (6), 621–637. Su, Y., Zhu, G., Miao, Z., Feng, Q., Chang, Z., 2009. Estimation of parameters of a biochemically based model of photosynthesis using a genetic algorithm. Plant Cell Environ. 32 (12), 1710–1723. Verbeeck, H., Samson, R., Verdonck, F., Lemeur, R., 2006. Parameter sensitivity and uncertainty of the forest carbon flux model FORUG: a Monte Carlo analysis. Tree Physiol. 26 (6), 807–817. Voßbeck, M., Clerici, M., Kaminski, T., Lavergne, T., Pinty, B., Giering, R., 2010. An inverse radiative transfer model of the vegetation canopy based on automatic differentiation. Inverse Probl. 26 (9), 095003. Wang, Y.-P., Jarvis, P.G., 1990. Description and validation of an array model—MAESTRO. Agric. For. Meteorol. 51 (3-4), 257–280. Wang, J., Li, X., Lu, L., Fang, F., 2013. Estimating near future regional corn yields by integrating multi-source observations into a crop growth model. Eur. J. Agron. 49, 126–140. Wania, R., Ross, I., Prentice, I.C., 2010. Implementation and evaluation of a new methane model within a dynamic global vegetation model: LPJ-WHyMe v1. 3.1. Geosci. Model. Dev. 3 (2), 565–584. Wood, S., Sebastian, K.L., Scherr, S.J., 2000. Agroecosystems. International Food Policy Research Institute and World Resources Institute, Washington (DC). Wu, S., Jansson, P.-E., Zhang, X., 2011a. Modelling temperature, moisture and surface heat balance in bare soil under seasonal frost conditions in China. Eur. J. Soil Sci. 62 (6), 780–796. Wu, S.H., Jansson, P.-E., Kolari, P., 2011b. Modeling seasonal course of carbon fluxes and evapotranspiration in response to low temperature and moisture in a boreal Scots pine ecosystem. Ecol. Modell. 222 (17), 3103–3119. Wu, M., Jansson, P.-E., Tan, X., Wu, J., Huang, J., 2016. Constraining parameter uncertainty in simulations of water and heat dynamics in seasonally frozen soil using limited observed data. Water 8 (2), 64. Wu, M., Scholze, M., Voßbeck, M., Kaminski, T., Hoffmann, G., 2019. Simultaneous assimilation of remotely sensed soil moisture and FAPAR for improving terrestrial carbon fluxes at multiple sites using CCDAS. Remote Sens. 11 (1), 27. Zhu, G.-F., Li, X., Su, Y.-H., Lu, L., Huang, C.-L., 2011. Seasonal fluctuations and temperature dependence in photosynthetic parameters and stomatal conductance at the leaf scale of Populus euphratica Oliv. Tree Physiol. 31 (2), 178–195. Zhu, Q., Zhuang, Q., 2014. Parameterization and sensitivity analysis of a process‐based terrestrial ecosystem model using adjoint method. J. Adv. Model. Earth Syst. 6 (2), 315–331.

Technometrics 33 (2), 161–174. Nylinder, J., Stenberg, M., Jansson, P.-E., Klemedtsson, Å.K., Weslien, P., Klemedtsson, L., 2011. Modelling uncertainty for nitrate leaching and nitrous oxide emissions based on a Swedish field experiment with organic crop rotation. Agric. Ecosyst. Environ. 141 (1-2), 167–183. Parton, W., 1996. The CENTURY Model, Evaluation of Soil Organic Matter Models. Springer, pp. 283–291. Piao, S., et al., 2010. The impacts of climate change on water resources and agriculture in China. Nature 467 (7311), 43. Pinty, B., et al., 2001. Radiation transfer model intercomparison (RAMI) exercise. J. Geophys. Res. Atmos. 106 (D11), 11937–11956. Pinty, B., et al., 2007. Retrieving surface parameters for climate models from moderate resolution imaging spectroradiometer (MODIS)-multiangle imaging spectroradiometer (MISR) albedo products. J. Geophys. Res. Atmos. (D10), 112. Pinty, B., Andredakis, I., Clerici, M., Kaminski, T., Taberner, M., Verstraete, M., Gobron, N., Plummer, S., Widlowski, J.L., 2011. Exploiting the MODIS albedos with the two‐stream inversion package (JRC‐TIP): 1. effective leaf area index, vegetation, and soil properties. J. Geophys. Res. Atmos. (D9), 116. Quillet, A., Frolking, S., Garneau, M., Talbot, J., Peng, C., 2013. Assessing the role of parameter interactions in the sensitivity analysis of a model of peatland dynamics. Ecol. Modell. 248, 30–40. Ramankutty, N., Evan, A.T., Monfreda, C., Foley, J.A., 2008. Farming the planet: 1. Geographic distribution of global agricultural lands in the year 2000. Global Biogeochem. Cycles 22 (1). Sándor, R., Barcza, Z., Hidy, D., Lellei-Kovács, E., Ma, S., Bellocchi, G., 2016. Modelling of grassland fluxes in Europe: evaluation of two biogeochemical models. Agric. Ecosyst. Environ. 215, 1–19. Santaren, D., Peylin, P., Viovy, N., Ciais, P., 2007. Optimizing a process‐based ecosystem model with eddy‐covariance flux measurements: a pine forest in southern France. Global Biogeochem. Cycles 21 (2). Scholze, M., Kaminski, T., Knorr, W., Blessing, S., Vossbeck, M., Grant, J., Scipal, K., 2016. Simultaneous assimilation of SMOS soil moisture and atmospheric CO2 in-situ observations to constrain the global terrestrial carbon cycle. Remote Sens. Environ. 180, 334–345. Senapati, N., Jansson, P.-E., Smith, P., Chabbi, A., 2016. Modelling heat, water and carbon fluxes in mown grassland under multi-objective and multi-criteria constraints. Environ. Model. Softw. 80, 201–224. Smith, B., Prentice, I.C., Sykes, M.T., 2001. Representation of vegetation dynamics in the modelling of terrestrial ecosystems: comparing two contrasting approaches within

306