Optimizing shrub parameters to estimate gross primary production of the sagebrush ecosystem using the Ecosystem Demography ( EDv 2 . 2 ) model

Gross primary production (GPP) is one of the most critical processes in the global carbon cycle, but is difficult to 10 quantify in part because of its high spatiotemporal variability. Direct techniques to quantify GPP are lacking, thus, researchers rely on data inferred from eddy covariance (EC) towers and/or ecosystem dynamic models. The latter are useful to quantify GPP over time and space because of their efficiency over direct field measurements and applicability to broad spatial extents. However, such models have also been associated with internal uncertainties and complexities arising from distinct qualities of the ecosystem being analyzed. Widely distributed sagebrush-steppe ecosystems in western North America are threatened by 15 anthropogenic disturbance, invasive species, climate change, and altered fire regimes. Although land managers have focused on different restoration techniques, the effects of these activities and their interactions with fire, climate change, and invasive species on ecosystem dynamics are poorly understood. In this study, we applied an ecosystem dynamic model, Ecosystem Demography (EDv2.2), to parameterize and predict GPP for sagebrush-steppe ecosystems in the Reynolds Creek Experimental Watershed (RCEW), located in the northern Great Basin. Our primary objective was to develop and parameterize a sagebrush 20 (Artemisia spp.) shrubland Plant Functional Type (PFT) for use in the EDv2.2 model, which will support future studies to model estimates of GPP under different climate and management scenarios. To accomplish this, we employed a three-tiered approach. First, to parameterize the sagebrush PFT, we fitted allometric relationships for sagebrush using field-collected data, gathered information from existing sagebrush literature, and borrowed values from other PFTs in EDv2.2. Second, we identified the five most sensitive parameters out of eleven that were found to be influential in GPP prediction based on previous 25 studies. Third, we optimized the five parameters using an exhaustive search method to predict GPP, and performed validation using observations from two EC sites in the study area. Our modeled results were encouraging, with reasonable fidelity to observed values, although some negative biases (i.e., seasonal underestimates of GPP) were apparent. We expect that, with further refinement, the resulting sagebrush PFT will permit explicit scenario testing of potential anthropogenic modifications of GPP in sagebrush ecosystems, and will contribute to a better understanding of the role of sagebrush ecosystems in shaping 30 global carbon cycles. Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2018-264 Manuscript under review for journal Geosci. Model Dev. Discussion started: 13 December 2018 c © Author(s) 2018. CC BY 4.0 License.


Introduction
Terrestrial gross primary production (GPP) is a major driver of the global carbon cycle as it plays an important role in regulation of atmospheric CO2 by offsetting anthropogenic CO2 emissions.GPP quantifies the rate of carbon uptake from the atmosphere through photosynthesis and is often presented as one of the most critical elements in global carbon research because of its high spatiotemporal variability (Chen, 2012;Zhao and Running, 2010).Since direct techniques to quantify GPP are lacking, researchers derive estimates using observations from eddy covariance (EC) towers or using ecosystem dynamic models (Dong et al., 2017;Turner et al., 2006;Zhao et al., 2005).Ecosystem models have been widely used to estimate terrestrial carbon flux and to project ecosystem dynamics over time and space (Dietze et al., 2014;Fisher et al., 2017), largely due to their efficiency over direct field measurements and their applicability to broader spatial scales.However, these models have also been associated with high levels of uncertainty and complexity associated with their applicability to distinct ecosystems.
Mainly, three different types of errors have been associated with these models: (1) process error arising from formulations in the model and associated parameters, (2) forcing error related to the quality of meteorological data, and (3) the initial ecosystem state at the start of the simulation (Antonarakis et al., 2014;Huntzinger et al., 2012).Initialization error is generally not an issue for long-term simulations, and researchers can minimize both forcing and initialization errors by using observational data rather than reanalysis data (Antonorakis et al., 2014;Medvigy et al., 2009).Indeed, process error is the most problematic, as it can mask uncertainties caused by forcing errors and can create potential bias in predictions.Fortunately, process error can be quantified and minimized by systematically comparing model predictions with observable ecosystem metrics (Braswell et al., 2005;Dietze et al., 2014;Medvigy et al., 2009).
Another critical limitation to widely applying ecosystem dynamic models is their suitability for a unique ecosystem for which they have not been parametrized.Semi-arid, non-forest ecosystems provide an excellent example of this limitation, including sagebrush (Artemisia spp.) ecosystems, one of the most widespread community types in North America.Sagebrush ecosystems hold high ecological and social value, but have been reduced to nearly half of their historical range and are declining at an alarming rate (Knick et al., 2003;Schroeder et al., 2004).Various factors have contributed to this decline, including land clearing, invasion of nonnative species such as cheatgrass (Bromus tectorum), and climate change, that have collectively altered vegetation composition, hydrological function, and wildfire frequency (Bradley, 2010;Connelly et al., 2004;McArthur and Plummer, 1978;Schlaepfer et al., 2014).In order to restore the sagebrush ecosystem, land managers have focused on suppressing fire, reducing flammable vegetation, controlling invasive species, and seeding native species (Chambers et al., 2014;McIver and Brunson, 2014).There are relatively few studies that have evaluated carbon flux in sagebrush ecosystems in response to prescribed fire or restoration activities, and most of them used observational data from EC stations.Previous studies identified temporal variation in net carbon exchange rates after restoration treatments, including documenting increases in carbon uptake in the early years after prescribed fire (caused by re-sprouting shrubs and fast growing grasses), followed by an eventual levelling off to pre-fire conditions (Cleary et al., 2010;Fellows et al., 2018).However, because of the paucity in EC station sites, coupled with the large spatial extent of the sagebrush ecosystems in the Great Basin, In this study, we applied an ecosystem dynamic model, Ecosystem Demography (EDv2.2) (Medvigy et al., 2009;Moorcroft et al., 2001), to parameterize and predict GPP for sagebrush ecosystems in an experimental watershed located in the northern Great Basin.The Great Basin is a ~500,000 km2 cold-desert region dominated by expansive, shrub-steppe ecosystems.Our primary objective was to develop preliminary sagebrush Plant Functional Type (PFT) parameters in the EDv2.2 model, based on sensitivity analysis and optimization, with respect to GPP prediction.EDv2.2 was originally developed for tropical forest ecosystems (Moorcroft et al., 2001), and it has been tested in boreal forests (Trugman et al., 2016), temperate forests (Antonarakis et al., 2014;Medvigy et al., 2009;Medvigy et al., 2013), and tundra (Davidson et al., 2011), however, its application to semi-arid shrubland ecosystems has not been explored.Preliminary parameterization of sagebrush PFT, from this study, will be a first step towards further studies on shrubland ecosystem using EDv2.2 or any other terrestrial models.

Ecosystem Demography (EDv2.2) model
The Ecosystem Demography (EDv2.2) is a process-based terrestrial biosphere model, which occupies a mid-point on the continuum of models from individual-based (or gap) models to area-based (or big-leaf) models (Fisher, 2010;Smith et al., 2001).Area based models like LPJ-DGVM (Lund-Potsman-Jena Dynamic Vegetation Model) (Sitch et al., 2003), and BIOME BGC (Running and Hunt, 1993as cited in Bond-Lamberty et al., 2014) represent plant communities with area-averaged representation of a plant function type (PFT) for each grid cell.The simplification and computational efficiency of these models make them widely applicable for regional ecosystem analysis, however, this advantage comes in trade-off with their failure to capture light competition, competitive exclusion, and disturbances (Fisher, 2010;Bond-Lamberty et al., 2014;Smith et al, 2001).On the other hand, individual-based models such as JABOWA (Botkin et al., 1972), and SORTIE (Pacala et al., 1993) represent vegetation at individual plant level thus making it possible to incorporate competition, coexistence, and disturbances.
The disadvantage of IBMs is that they are often confined to limited spatial and temporal scales due to added computational burden.EDv2.2 is a cohort based model where individual plants with similar properties, in terms of size, age, and function, are grouped together to reduce the computational cost while retaining most of the dynamics of IBMs.Each cohort is defined by a PFT, number of plants per unit area, and dimensions of a single representative plant like diameter, height, structural biomass, and live biomass (Fisher et al., 2010).In EDv2.2, land surface is composed of a series of gridded cells, which experience meteorological forcing from a corresponding gridded data or from a coupled atmospheric model (Medvigy, 2006).The mechanistic scaling from individual to the region is achieved through size and age structured partial differential equations that closely approximate mean behaviour of stochastic gap model (Medvigy et al., 2009;Moorcroft et al., 2001).Each grid cell is subdivided into a series of dynamic horizontal tiles, which represent locations that had experienced similar disturbance history and has an explicit vertical canopy structure.This mechanism helps capture both vertical and horizontal distribution of vegetation structure and compositional heterogeneity compared to area-based models (Kim et al., 2012;Moorcroft et al., 2001;Moorcroft et al., 2003;Sellers et al., 1992).EDv2.2 consists of multiple sub-models for plant growth and mortality, phenology, disturbance, biodiversity, hydrology, land surface biophysics, and soil biogeochemistry, to predict short-term and long-term ecosystem flux and to represent natural and anthropogenic disturbances (Kim et al., 2012;Medvigy et al., 2009;Zhang et al., 2015).Sub-models in EDv2.2 rely mostly on many PFT-specific parameters, representing unique attributes of that particular group of species, to resolve the stated biological processes (Knox et al., 2015).EDv2.2 has parameters defined for 17 different PFTs including grasses (C3 & C4), conifers, and deciduous trees (temperate & tropical), and agricultural crops.In this study, we identified parameters for sagebrush (shrub) ecosystem to simulate it in the model as a new PFT.Furthermore, since we focussed on GPP prediction, we selected eleven different parameters related to plant ecophysiology and biomass allocation from overall to conduct sensitivity and optimization study.We used similar studies (Dietze et al., 2014;Fisher et al., 2010;LeBauer et al., 2013;Medvigy et al., 2009;Mo et al., 2008;Pereira et al., 2017), preliminary sensitivity analyses, and consultation with other developers and users of the EDv2.2 model to select these parameters.These included maximum photosynthetic capacity at 15⁰ C (Vm0), specific leaf area (SLA), fine root turnover rate, leaf turnover rate, storage turnover rate, slope of stomatal conductance-photosynthesis relationship (stomatal slope), ratio of fine roots to leaves (Q), water conductance, cuticular conductance, growth respiration factor (GRF), and leaf width.
As we can find detailed descriptions of sub-models of EDv2.2 in the existing literature (Medvigy et al., 2009;Moorcroft et al., 2001), here we are trying to describe the ones related to the parameters we have used in this study.The ecophysiological sub-model has a coupled photosynthesis and stomatal conductance scheme developed by Farquhar and Sharkey (1982) and Leuning (1995) respectively takes care of leaf-level carbon and water fluxes.Leaf-level carbon demand of C3 plants is determined by the minimum of light-limited rate (Je) and Rubisco-limited rate (Jc), and Vm0 controls the later as given by Eq.(1) after being scaled to a given temperature.
where,   (  ) is the maximum capacity of Rubisco to perform carboxylase function at a given temperature Tv scaled from Vm0,   is the intercellular CO2 concentration, Г is the compensation point for gross photosynthesis, K1 is the Michaelis-Menten coefficient for CO2, and K2 is proportional to the Michaelis-Menten coefficient for O2.Stomatal conductance which is modeled using Leuning (1995), a variant of Ball Berry model (Eq.2), is influenced by stomatal slope and cuticular conductance.where,   is photosynthetic rate,   is stomatal conductance for water,  is stomatal slope,  is cuticular conductance,  0 is empirical constant,   is water vapour deficit and CO2 concentration within leaf boundary, and Г is as described above.
Stomatal control is also affected by soil moisture supply term, which is a function of soil moisture, fine root biomass, and water conductance.Water and CO2 concentrations within the leaf boundary layer is influenced by leaf width along with other factors like wind speed, leaf area index, and molecular diffusivity of heat.Specific leaf area (SLA) has units of leaf are per unit leaf carbon and is used scale up leaf-level fluxes to canopy-level fluxes.Relationship between growth respiration and net photosynthesis is controlled by the growth respiration factor.In EDv2.2, while leaf biomass is determined by allometric equations based on diameter, fine root biomass is defined by a ratio of leaves to fine roots.Leaf turnover and fine root turnover rates together influence overall litter turnover rate, even though in deciduous trees dropping of leaves also affects this rate.
Turnover rate of stored leaf pool and storage respiration depends on storage turnover rate, size of stored leaf pool, and storage biomass.

Study area
We initialized and performed parameter optimization for sagebrush ecosystems in the EDv2.Seyfried et al., 2000) and is characterized as having light cattle grazing (AmeriFlux, 2018).Another AmeriFlux tower, US-Rws, is located at 43.1675 N and 116.7132W in the Nancy Gulch drainage, within about 2.5 km distance to the northeast from the study polygon.This area is dominated by Wyoming big sagebrush (A.tridentata ssp.wyomingensis) and bluebunch wheatgrass (Pseudoroegneria spicata) (Stephenson, 1970).Hereafter, these two sites are designated as LS (for low sagebrush) and WBS (for Wyoming big sagebrush), respectively.(Homer, et al., 2015) in the background.The inset map shows the location of the study area within the Great Basin (LCC, 2018).

Inventory and EC tower data
A field inventory dataset of sagebrush shrubs from RCEW recorded in 2014 (Glenn et al., 2017) was used to fit the allometric equations in EDv2.2, which were developed for temperate PFTs, and to initialize the ecosystem structure for the model simulation.Variables used to fit allometric equations for the sagebrush shrub included volume, crown diameter, height, and stem diameter.EDv2.2 was originally developed for tropical forests, and thus typically specifies allometric relationships in terms of diameter at breast height (dbh).However, this length-scale variable has limited application to shrubs of the sagebrushsteppe ecosystem, which rarely exceed 1.5 m in height.Thus, we developed a substitute length-scale variable for dbh that effectively corresponds to shrub volume.To accomplish this, shrub volume was first calculated using crown area (characterized as an ellipse, and approximated with semi-major and semi-minor axis lengths) and height, and the cube root of this volume was then used as the characteristic length-scale variable required to parameterize allometric relationships in EDv2.2.We identified the coefficients in allometric equations (Table 1) for shrub height, leaf biomass, structural biomass, canopy area, and wood area index as a function of this cube-root of volume (used as DBH in the equation).GPP data for 2015 and 2016 water years were derived from the LS and WBS EC stations (Fellows et al., 2017) using the REddyProc software in R (Reichstein et al., 2005) to fill and partition net ecosystem exchange (NEE) into ecosystem respiration and GPP.

Meteorological forcing data
Outputs from a long-term, high resolution climate reanalysis obtained from the Weather Research and Forecast (WRF) model (Skamarock et al., 2008) were used to provide meteorological forcing data for the EDv2.2 model (Table 2).The WRF outputs correspond to atmospheric temperature and specific humidity at 2 m height, wind speed at 10 m height, downward shortwave radiation and long-wave radiation at ground surface, surface pressure and accumulated precipitation (Flores, et al., 2016).The spatial and temporal resolutions of the data are 3 km and 3 hours, respectively.The EDv2.2 model then partitions shortwave radiation into direct and diffuse, visible and near-infrared components as summarized by Weiss and Norman (1985).We obtained these forcing data for nine years from 2006 to 2014 at the WRF pixel corresponding to the polygon study area and bounding the LS EC site location.For each year of EDv2.2 simulation, a random year of meteorological forcing data was chosen from the available range of data.

Initial parameterization and sensitivity analysis
We identified initial shrub PFT parameters based on field allometric equations, previous research studies on the sagebrush ecosystem (Cleary et al., 2010;Comstock and Ehleringer, 1992;Gill and Jackson, 2000;Li et al., 2009;Olsoy et al., 2016;Qi et al., 2014), and existing PFT parameters in the EDv2.2 model for C3 grass, northern pines, and late conifers (Table S1 in the Supplement).The initial ecosystem state for the model run was designated to be a single sagebrush cohort with an average cube root volume (diameter) of 0.6 m, average height of 0.52 m, and density of 1 plant/m 2 representing 2014 field inventory data.The soil column was configured to be 2.3 m deep with 9 vertical layers and a free-drainage lower boundary.
Corresponding to a gravelly loam soil in the study site (USDA, 2018), we used a soil texture with 55% sand, 25% silt, and 20% clay.Initial soil moisture was set to near saturation with no temperature offset, and the initial atmospheric carbon dioxide level was set at 370 ppm.The EDv2.2 model was then run with these initial settings and initial shrub PFT parameters values for an 8-year simulation period.
We used sensitivity index (SI) suggested by Hoffman and Gardner (1983) (Eq. 3) to perform preliminary one at a time sensitivity analysis and rank the parameters.Since, this index is highly affected by the extreme values of parameters being studied, it is recommended that the parameter range cover the entire range of possible values.SI has been used in different areas of studies including ecology (Waring et al., 2016) and hydrology (Wambura et al., 2015), mostly to assess the effect of parameters on target variables, and sometimes to reduce the number of variables for further analysis.
where,  is sensitivity index,   is the value of GPP corresponding to the simulation with the maximum value of a parameter, and   is the value of GPP corresponding to the simulation with the minimum value of a parameter.We identified minimum and maximum possible values for each of the selected parameters based on previous sensitivity and optimization studies, the range of parameters for other PFTs in EDv2.2, and our preliminary sensitivity analyses (Table 3).
EDv2.2 was then run for an eight year period with both minimum and maximum values of each parameter while keeping all other parameters constant.The average daily GPP outputs throughout the simulation years for maximum and minimum values of parameters were used to get   and   respectively.Based on SI, we limited optimization and validation to the five most sensitive parameters from the list of eleven, to keep time and computing performance manageable.

Optimization and validation
In the third step, an optimization of the five selected parameters was performed using an exhaustive search method within the specified range of values.This was performed to identify the best values for five selected parameters for each EC station in predicting GPP.We ran 900 simulations with a unique combination of parameter values for 15 years, at which point it was assumed to reach an equilibrium with climate.EDv2.2 simulations were configured to allow for growth of the C3 grass, northern pines, and late conifers together with the shrub PFT.This was done because although the vegetation assemblage in the footprints of flux sites is primarily composed of sagebrush and grasses, conifers are present in some parts of the experimental watershed (Seyfried et al., 2000).For each simulation, we calculated a skill score, Nash-Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970), to compare the final year simulated GPP with those derived from both the LS and WBS EC stations for 2016.Although, NSE is closely related to root mean square error (RMSE) (or mean square error, MSE), the skill score from it can be interpreted as comparative ability of the model over a baseline model, which is the mean of site observations in this case.While the RMSE value depends on the unit of predicted variables, which can vary from 0 to infinity, the NSE is dimensionless and varies from negative infinity to 1 (Krause et al, 2005;Gupta et al, 2009).NSE is calculated using Eq. ( 4): where,   is observation,   is predicted value,  ̅ is mean of observation, and n is number of observations.For both EC stations, we selected the 10 best simulations based on NSE scores and computed ensemble means of all five parameter values.
We then ran the EDv2.2 model using the ensemble mean parameter values for both EC sites.The simulated GPP from ensemble mean parameter values and the best case (highest NSE) were then compared against EC site data from 2015, which was withheld from the optimization as a means of providing an independent validation.

Initial parameterization and sensitivity analysis
For the model run based on the initial values of parameters (Table S1 of Supplement), the 8-year simulations produced an annual cycle in GPP that decreases in amplitude during the initial 1-3 years, and remains at a level of approximately 0.07 KgC/m 2 /yr in the remaining years (Fig. 2a).Observed GPP in 2016 were 0.51 KgC/m 2 /yr and 0.38 KgC/m 2 /yr for the LS and WBS sites, respectively.This result was significantly lower than the observed GPP from either of the EC sites, and thus warranted changes in initial values for one or more parameters.
Based on the SI ranking,  0 , , stomatal slope, fine root turnover rate, and Q-ratio were identified as the top five sensitive parameters compared to the other parameters explored (Fig. 2; Table 3).Related studies (Dietze et al., 2014;Medvigy et al., 2009;Pereira et al., 2017;Zaehle et al., 2005) have also identified similar model parameters being important in estimating GPP.In our study, higher parameter values of SLA,  0 , and stomatal slope, resulted in higher GPP estimates (Fig. 2b and d), whereas for fine root turnover rate and Q-ratio, higher parameter values produced lower GPP (Fig. 2e and f).The impact of shifts in SLA,  0 , and stomatal slope values are observed from the very beginning of the simulations, while changes in fine root turnover rate and Q-ratio parameters start to show differences from roughly 3-4 years after the initial model run.Although not ranked in the top five, leaf turnover rate, cuticular conductance, and growth respiration factor also had considerable influences over GPP (Table 3).Cuticular conductance (µmolm -2 s -1 ) 10 3 10 2 10 7 0.501 Barnard and Bauerle (2013): Duursma et al., (2018) Water Conductance (ms -1 kgCroot -1 ) 1.9 × 10 -5 1.9 × 10 -6 1.9 × 10 -4 0.

Optimization and validation
For our exhaustive search of parameter values, we limited search domains for parameters based on previous studies and the result of our sensitivity analysis.SLA search limits were largely based on Olsoy et al. (2016), who suggested a range of 3 to 6 m 2 /Kg for sagebrush SLA, but who also hinted at variation due to regional and seasonal variation.Similarly, limits for  0 were extended slightly beyond Comstock and Ehleringer's (1992) recommendations for Great Basin shrubs, and the upper limit for stomatal slope was extended slightly beyond that used by Oleson et al. (2013) for a shrub PFT in the Community Land Model (CLMv4.5).We set search domains for Q-ratio based on a leaf and root biomass study of sagebrush by Cleary et al. (2010), and fine root turnover ratio was based on results from a study on Artemisia ordosica in a semi-arid region of China (Li et al, 2009).Interval distances (or 'steps') were calculated to equally space out the range between the maximum and minimum of each parameter for a given number of intervals (Table 4).Parameters identified as exerting more control on GPP prediction were assigned higher number of steps, resulting in the following: five steps of SLA and  0 , four steps for stomatal slope, and three steps for Q-ratio and fine root turnover rate.Among 900 simulations for unique parameter value combinations, 180 cases which did not provide model optimization results because of numerical instabilities (with GPP approaching zero), were excluded from subsequent analysis.We selected ten simulations with the best NSE scores for both LS and WBS sites (Table S2 and Fig. S1 in the Supplement) and determined ensemble means of parameter values for these sites (Table 5).We then ran EDv2.2 to predict GPP using ensemble mean parameter values for each of the EC stations (hereafter, the 'ensemble case').Among the ten best simulations selected for each EC sites, two of them were commonly selected in both sites.One of them was ranked top with highest NSE score (hereafter, the 'best case') and the other one was ranked as top fifth for LS and top second for WBS site.Both of these common simulations, showed gradual growth of C3 grass through the simulation years even though we initialized model with only shrub PFT.In the final year of simulation, the 'best case' had about 51% of total GPP coming from C3 grass and the other one had about 43% from it (Table S2 and Fig. S2 in the Supplement).Optimized parameter values were considerably different between the best case and ensemble cases for both stations, possibly suggesting interaction effects among the parameters (Table 5).When parameters from ensemble means between two stations were compared, mean  0 for LS was higher by 3 µmolm -2 s -1 than that for WBS.Mean parameter for stomatal slope was lower by 1 for LS than for WBS site.Another clear difference was with ensemble mean for Q-ratio which was higher by 0.70 for LS than for WBS.But, we did not find such differences for SLA and fine root turnover rate.Figure 3 compares simulated GPP from the final modeled year (from the ten best simulations and the ensemble case) with the observed GPP in 2016 from each EC station.Optimization results for the LS site in Fig. 3a show that simulated GPP matches well with observed data for most days, except during the spring, during which a clear peak in observed GPP was not captured by the simulation results.A lower spring peak in GPP was observed at the WBS site (Fig. 3b) and was far more comparable to simulation results.The impact of the spring mismatch in the LS site was such that, despite some over-prediction during the fall season for WBS simulations, Bias and NSE were better for WBS than for LS data (Table 5).However, for both EC site comparisons, most simulations resulted in a negative bias, and optimization NSEs for the ensemble cases were not as good as for the best cases.For validation of the parameter estimates, we compared the simulated GPPs produced from; (i) the best case, and (ii) the ensemble case, for both LS and WBS sites, with observed GPP from the respective locations in 2015.When comparing Bias and NSE derived from validation with that from optimization, we found that simulations for the WBS site data performed well and with relative parity, while those for the LS site performed poorly, especially given the negative NSE values for both best case (-0.193) and ensemble mean case (-0.183) for validation (Table 6).In contrast, validation results for the WBS site showed the model performed well for both cases, although the NSE from the best case (0.408) was higher than that from ensemble case (0.260).Poor validation results for the LS site could be due to missing data (69 days) from the 2015 observations (i.e., the validation dataset), or possibly due to the higher difference in observed mean annual GPP values between 2015 and 2016 for the LS site than for the WBS site.Results shown in Fig. 4 also indicate that the simulated daily GPP matched well with the validation data for the WBS site, but was a relatively poor match for the LS site.In particular, for the LS site, predicted GPP was significantly lower than observed for most of the spring and summer months (Fig. 4c).The validation data (from 2015) for the LS site had clearly higher GPP values for late summer days (Fig. 4a and c), thus resulting in higher negative Bias and negative NSE compared to the optimization results.Predicted daily GPP for this site during the remaining months, however, is comparable with observed values.In contrast, predicted daily GPP for the WBS site matched well with the validation GPP data through most of the year (Fig. 4b and d), with slight inconsistencies during September (under estimation) and October (over estimation).Compared to 5 the best case, the ensemble case simulation for the LS site performed slightly better for most months, much better for July and August, but worse for October and November (Fig. 4c).For the WBC site (Fig. 4d), the clearest differences between the two cases were for April, May, July, and August, during which the best case simulation strongly outperformed the ensemble case.

Discussion
Our sensitivity analysis results were similar to previous studies (Dietze et al., 2014;LeBauer et al., 2013;Medvigy et al., 2009), wherein parameters  0 , SLA, fine root turnover rate, and stomatal slope were found to be the most influential in determining carbon flux or primary productivity.A high variation in parameter values was observed among the simulations that resulted in the best NSE values (Table S2 in the Supplement).The effects of some parameters on GPP prediction were not the same when altered individually or simultaneously with other parameters.For instance, the sensitivity analysis showed an increase in GPP with an increase in  0 (Fig. 2) and, since our initial GPP prediction was very low, we would expect higher  0 for better prediction.However, when all five parameters were optimized simultaneously, a generally lower value of  0 produced the best NSE (Table 5).Likewise, despite the sensitivity analysis suggesting higher GPP with lower fine root turnover ratio and Q-ratio, the best results were obtained with the same initial values (maximum values in the search domain) for these parameters, Nonetheless, top ten best parameter combinations exhibited greater variation of these parameters for either EC site, and resulted in slightly lower mean values than the initial ones.This suggests interaction effects are dominant over first order effects of the studied parameters, and there is likely nonlinear dependence among them.
GPP simulations from this study demonstrated a closer similarity with observations from the WBS site compared to the LS site, even though the WBS EC station is outside of the WRF pixel used in this analysis.This could be due to slight differences in soils and hydraulic conditions in the field compared to those conditions used for initialization.Moreover, variation between morphological characteristics of the vegetation at the LS and WBS EC towers (characterized by low sagebrush and Wyoming big sagebrush, respectively), such as differences in common plant heights and flowering seasons, may have resulted in the observed differences in GPP (Howard, 1999;USDA, 2018).Since Wyoming big sagebrush is the dominant species in the Reynold Creek Watershed area (Seyfried et al., 2000), meteorological forcing data used in this simulation, as well as the allometric equations fitted for sagebrush (representing most areas of RCEW), could be favoring the more realistic growth pattern of this species in the model (e.g., Fig. 3 and 4).
Additionally, differences in the phenology of the associated grass species between the two sites could result in differences in seasonal and annual productivity (Cleary et al., 2015).For instance, the perennial grass at the LS site is Sandberg bluegrass, which is photosynthetically active in early spring and senesces by early summer (USDA, 2016), and thus may have contributed to observed spring GPP peak at the LS site.In contrast, the associated grass at the WBS site, bluebunch wheatgrass, does not typically senesce by early summer.Indeed, the best optimization result from this study showed a gradual increase in C3 grass through the years along with the shrubs, with about 51% of GPP coming from the former by the final year.We could not validate how close this result was in terms of the actual species composition and ecosystem dynamics of the EC sites, as we did not have GPP observations for unique PFTs.However, the emphasis of this study was to optimize shrub PFT parameters, rather than C3 grass PFT parameters for the study area.Although we would expect that simultaneous optimization of both grass and shrub PFTs would result in improved representation of the vegetation composition in the study area, it would also increases the number of parameters required, potentially complicating the process of optimization and validation unique to each PFT.Moreover, several studies suggest that the parameters  0 and SLA vary considerably across seasons (Groenendijk et al., 2011;Kwon et al., 2016;Olsoy et al., 2016;Zhang et al., 2014).The mismatch in daily GPP pattern between simulated and flux tower data for specific seasons could be partly attributed to the lack of the model's ability to address these seasonal deviations correctly.Like most other terrestrial biosphere models, EDv2.2 does not incorporate seasonal variation in  0 , SLA, or other model parameters (Medvigy et al., 2009).Finally, we can achieve better results in parameter optimization and GPP prediction of sagebrush ecosystem, by making some advances in our methods.We can adopt some robust sensitivity (including variance decomposition, first order and second order analysis) (Zhang et al., 2017) and optimization (including cost function, gradient descent, and uncertainty analysis) (Richardson et al., 2010) methods to fine tune the sagebrush PFT parameters.Similarly, instead of relying on a single year of observation data for optimization and/or validation, we can use multiple years of data that would take into account the inter annual variability normally observed in ecosystem fluxes.

Conclusions
In this study, the Ecosystem Demography (EDv2.2) model was used to parameterize shrub PFT parameters and predict GPP for a sagebrush ecosystem in the Great Basin.Initial shrub PFT parameters were identified based on allometric equations fitted with field data, previous studies on sagebrush and shrubs in the sagebrush-steppe, and other PFTs (C3 grass, northern pines, and late conifers) in EDv2.2.The WRF model was used to acquire and force simulation with meteorological inputs to predict GPP.The simulation with initial shrub PFT parameters showed annual decline in GPP for 1-3 years and remained at a low level (compared to observed GPP data) for the remaining simulation period.Sensitivity analysis suggested Vm0, SLA, stomatal slope, fine root turnover rate, and Q-ratio ranked top five in influencing GPP prediction, which agrees with previous studies.
An exhaustive search was performed over constrained domains to explore the optimum combination of parameters to predict GPP.This led to identification of parameter values for best case and ensemble mean (of the ten best cases) cases optimized for the LS and WBS EC sites, using the NSE.Even though the model predicted daily GPP quite well, mostly negative bias was observed in predictions, and there was mismatch during the spring months.Validation results showed better performance by parameters optimized for WBS site than those done for LS site in GPP prediction.The difference in the local site vegetation community and the overall dominance of Wyoming big sagebrush in the study area and in the Great Basin may explain why the GPP predictions were closest to the WBS site.Similarly, the limitation of EDv2.2 in incorporating seasonal variation of parameters like Vm0 and SLA, could also be attributed to its poor predictions for spring seasons.
Our identification of coefficients for allometric equations coupled with the other parametrization of a semiarid shrub PFT for EDv2.2 will permit exploration of additional research questions.For instance, we can run EDv2.2 at regional scales with optimized parameters to model the spatiotemporal dynamics of sagebrush community composition and ecosystem flux, under different climate and ecological restoration scenarios.Another direction is to optimize C3 grass PFT parameters in EDv2.2, simultaneously with shrub PFT parameters, by using multiple years of observation data to characterize inter annual variation.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-264Manuscript under review for journal Geosci.Model Dev. Discussion started: 13 December 2018 c Author(s) 2018.CC BY 4.0 License. the collective effects of restoration activities, fire, climate change, and invasive species on the spatiotemporal dynamics of structure, composition, and function of ecosystem are poorly understood.
2 model using field data and two EC station sites located in the Reynolds Creek Experimental Watershed (RCEW) and Critical Zone Observatory (CZO), operated by the USDA Agricultural Research Service (Fig. 1).We used a 200 m x 200 m polygon with center location of 43.15 N and 116.72 W and a mean elevation of 1583 m.The AmeriFlux US-Rls EC station, located at 43.1439 N and 116.7356W within the Lower Sheep Creek drainage is approximately 0.7 km from the center of our study site.The area within the footprint of this tower is dominated by low sagebrush (Artemisia arbuscula) and Sandberg bluegrass (Poa secunda)(Stephenson, 1970;

Figure 2 .
Figure 2. Simulated daily GPP outputs from 1-8 years for the study location with (a) initial values of all five parameters, and 5 (b-f) maximum (green), minimum (blue), and initial (red) parameter values for SLA,  0 , stomatal slope, Q-ratio, and fine root turnover rate.

Figure 3 .
Figure 3. Optimization of daily GPP (KgC/m 2 /yr) for water year 2016 based on EC station observation data from a) LS and b) WBS EC towers.Spring is shown as roughly 170-255 DOY.

Figure 4 .
Figure 4. Validation of GPP (KgC/m 2 /yr) using best case and ensemble case against respective EC station observation data 10 from water year 2015.a) daily GPP for LS, b) daily GPP for WBS, c) monthly GPP for LS, and d) monthly GPP for WBS.
Information about the range comes from range of values for other PFTs in EDv2.2, and our preliminary sensitivity analysis

Table 4 .
Minimum value, maximum value, interval size, and number of steps for each parameter used in optimization.

Table 5 .
Optimized parameter values from best cases (highest NSEs) and ensemble means (mean of top 10 simulations) for LS and WBS EC stations.

Table 6 .
Bias and NSE for optimization and validation of GPP using parameter values from the best case and the ensemble case for both EC stations.