Step-wise modifications of the Vegetation Optimality Model

The Vegetation Optimality Model (VOM, Schymanski et al., 2009, 2015) is an optimality-based, coupled water-vegetation model that predicts vegetation properties and behaviour based on optimality theory, rather than calibrating vegetation properties or prescribing them based on observations, as most conventional models do. In order to determine whether optimality theory can alleviate common shortcomings of conventional models, as identified in a previous model inter-comparison study along the 5 North Australian Tropical Transect (NATT) (Whitley et al., 2016), a range of updates to previous applications of the VOM have been made for increased generality and improved comparability with conventional models. To assess in how far the updates to the model and input data would have affected the original results, we implemented them one-by-one while reproducing the analysis of Schymanski et al. (2015). The model updates included extended input data, the use of variable atmospheric CO2-levels, modified soil properties, 10 implementation of free drainage conditions, and the addition of grass rooting depths to the optimized vegetation properties. A systematic assessment of these changes was carried out by adding each individual modification to the original version of the VOM at the flux tower site of Howard Springs, Australia. The analysis revealed that the implemented changes affected the simulation of mean annual evapo-transpiration (ET) and gross primary productivity (GPP) by no more than 20%, with the largest effects caused by the newly imposed free drainage 15 conditions and modified soil texture. Free drainage conditions led to an underestimation of ET and GPP, whereas more finegrained soil textures increased the water storage in the soil and resulted in increased GPP. Although part of the effect of free drainage was compensated for by the updated soil texture, when combining all changes, the resulting effect on the simulated fluxes was still dominated by the effect of implementing free drainage conditions. Eventually, the relative error for the mean annual ET, in comparison with flux tower observations, changed from an 8.4% overestimation to an 10.2% underestimation, 20 whereas the relative errors for the mean annual GPP stayed similar with a change from 17.8% to 14.7%. The sensitivity to free drainage conditions suggests that a realistic representation of groundwater dynamics is very important for predicting ET and GPP at a tropical open-forest savanna site as investigated here. The modest changes in model outputs highlighted the robustness of the optimization approach that is central to the VOM architecture. 1 https://doi.org/10.5194/gmd-2021-151 Preprint. Discussion started: 18 June 2021 c © Author(s) 2021. CC BY 4.0 License.

1 Introduction 25 Novel modelling approaches that are able to explicitly model vegetation dynamics, may lead to an overall improved understanding of flux exchanges with the atmosphere. Recent model inter-comparison studies also reveal that novel model approaches are needed, especially related to vegetation dynamics (e.g. Whitley et al., 2016). Therefore, we use here optimality theory to predict the variation and dynamics of vegetation cover, root systems, water use and carbon uptake without the need for sitespecific input about vegetation properties. The theory is based on the premise that the net carbon profit (NCP), which is the 30 difference between carbon assimilated by photosynthesis and carbon expended on construction and maintenance of all the plant tissues needed for photosynthesis and water uptake and storage, is an appropriate measure of plant fitness, given that assimilated carbon is a fundamental resource of plant growth, development, survival and reproduction. The theory further assumes that construction and maintenance costs of plant organ functionality are general and therefore transferable between species and sites. Hence, the costs and benefits at different sites are determined in a consistent way, leading to vegetation properties that 35 solely depend on physical conditions, such as meteorological forcing, soils and hydrology. As a result, this leads to a systematic and consistent explanation of vegetation behaviour under different external conditions at different sites.
These optimality principles were employed in the Vegetation Optimality Model (VOM, Schymanski et al., 2009Schymanski et al., , 2015. The VOM is a coupled water-vegetation model that optimizes vegetation properties to maximize the Net Carbon Profit (NCP) in the long-term (20-30 years) for given climate and physical properties at the site under consideration. The NCP is defined as 40 the difference between the total carbon amount assimilated by photosynthesis, and the total carbon costs for the maintenance of leaf area, photosynthetic capacity and root surface area, as described in Schymanski et al. (2007Schymanski et al. ( , 2008b. The VOM has been previously applied by Schymanski et al. (2009) and Schymanski et al. (2015) at Howard Springs, a flux tower site in the North Australian Tropical Transect (NATT, Hutley et al., 2011). The NATT consists of multiple flux tower sites along a precipitation gradient from north to south, which allows for a more systematic testing of the VOM under different climatological 45 circumstances. The NATT has been used previously in an intercomparison of terrestrial biosphere models (TBMs) by Whitley et al. (2016), which revealed that lacking or wrong vegetation dynamics and incorrect assumptions about rooting depths have a strong influence on the performance of state-of-the-art TBMs. In contrast to these TBMs, the VOM predicts rooting depths and vegetation dynamics, and provides therefore a novel approach for the simulation of these savanna sites.
To assess in how far the optimality-based simulation of rooting depths, tree cover and vegetation dynamics may alleviate the 50 shortcomings of TBMs identified by Whitley et al. (2016), we propose to run the VOM using the same input data and similar physical boundary conditions at the different sites as in Whitley et al. (2016). In the simulations by Whitley et al. (2016), all TBMs were run under the assumption of a freely draining soil column, even though studies suggest an influence of groundwater on the resulting fluxes (York et al., 2002;Bierkens and van den Hurk, 2007;Maxwell et al., 2007). Free draining conditions can be mimicked in the VOM by setting drainage parameters to very fast drainage and the critical water table for the onset of 55 drainage very low (see below).
In previous applications of the VOM (Schymanski et al., 2009(Schymanski et al., , 2015, atmospheric CO 2 concentrations were assumed constant over the entire modelling period, whereas for inter-comparison with other models, we should use measured CO 2 -levels, which have increased considerably over the past years (Keeling et al., 2005). Previous applications of the VOM also prescribed a grass rooting depth of 1 m, arguing that due to the presence of a hard pan at Howard Springs, only tree roots could 60 penetrate into deeper layers. Since we do not know if such features exist at the other sites along the NATT, rooting depth of seasonal vegetation should be optimised in a similar way as that of perennial vegetation.
To assess to what extent the various changes influence the VOM-results, a new set-up of the VOM was applied to the same flux tower site in Australia, Howard Springs, as in Schymanski et al. (2009Schymanski et al. ( , 2015. This technical note describes the changes to the VOM since its last application by Schymanski et al. (2015), and how they affect the results of the VOM one-by-one and 65 in combination.

Methodology
All steps in the process, from pre-and post-processing to model runs, were done in an open science approach using the RENKU 1 platform. The workflows including code and input data can be found online 2 . In the following, we briefly describe the study site, the VOM, and the various modifications done in this study, compared to Schymanski et al. (2015).

Study site
The study site is Howard Springs (How-AU), which was previously used by Schymanski et al. (2009) andSchymanski et al. (2015) and provides a long record of carbon dioxide and water fluxes starting from 2001 (Beringer et al., 2016). Howard Springs is the wettest site (average precipitation of 1747 mm/year (SILO Data Drill, Jeffrey et al., 2001, calculated for 1980-2017 along the North Australian Tropical Transect (NATT, Hutley et al., 2011), which has a strong precipitation gradient from north 75 to south, with a mean annual precipitation around 500 mm/year at the driest site. The vegetation at Howard Springs consists of a mostly evergreen overstorey (mainly Eucalyptus miniata and Eucalyptus tetrodonta) and an understorey dominated by annual Sorghum and Heteropogon grasses. The soils at Howard Springs are well-drained red and grey kandosols, and have a high gravel content and a sandy loam structure.

80
The Vegetation Optimality Model (VOM, Schymanski et al., 2009Schymanski et al., , 2015 is a coupled water and vegetation model, that optimizes vegetation properties by maximizing the Net Carbon Profit. The model code and documentation can be found online 34 and version v0.5 5 of the model was used here. A general description is given below, whereas a more detailed description can be found at Schymanski et al. (2009Schymanski et al. ( , 2015.

Vegetation model 85
The VOM schematizes the ecosystem as two big leaves, one representing the seasonal vegetation (grasses) and one representing the perennial vegetation (trees). Photosynthesis was modelled according to Schymanski et al. (2007), who simplified the canopy-gas exchange model of von Caemmerer (2000) for C 3 -plants. The model computes CO 2 -uptake as a function of irradiance, atmospheric CO 2 -concentrations, temperature, photosynthetic capacity, projected foliage cover and stomatal conductance, whereby photosynthetic capacity, projected foliage cover and stomatal conductance are optimized dynamically in a 90 way to maximize the overall Net Carbon Profit (NCP) of the vegetation over the entire simulation period. Optimization is possible due to the carbon costs associated with each of these variables: photosynthetic capacity is linked to maintenance respiration, projected cover to the turnover and maintenance of leaf area, while stomatal conductance is linked to transpiration (depending on the atmospheric vapour pressure deficit) and hence root water uptake costs and limitations. Root water uptake is modelled following an electrical circuit analogy, where the water potential difference between the plant and each soil layer drives the 95 flow. Here, the root surface area and soil hydraulic conductivity in each soil layer determine the resistance (Schymanski et al., 2008b). The root surface area, in return, generates carbon costs for maintenance and the vertical distribution in the soil profile is optimized in a way to satisfy the canopy water demand with the minimum possible total root surface area.  Hutley et al. (2011), Hutley (2015 and Whitley et al. (2016), with Eucalyptus (Eu.), Erythrophleum (Er.), Hetropogan (He.). Meteorological data is taken from the SILO Data Drill (Jeffrey et al., 2001) for the model periods of 1-1-1980 until 31-12-2017, with the potential evaporation calculated according to the FAO Penman-Monteith formula (Allen et al., 1998). The ratio of the net radiation Rn with the latent heat of vaporization λ mutiplied with the precipitation P , is defined here as the aridity Rn/λP . Tree cover is determined as the minimum value of the mean monthly projective cover based on fPAR-observations (Donohue et al., 2013). The maximum grass cover was found by subtracting the tree cover from the remotely sensed projective cover.

Long-term optimization
The rooting depths of the perennial trees and the seasonal grasses (y r,p and y r,s ) as well as the foliage projected cover of 100 the perennial vegetation (M A,p ) are derived by optimizing these properties for the long-term, assuming that these do not vary significantly during the simulation period (20-30 years). Similarly, water use strategies of both the perennial and the seasonal vegetation component are assumed to be a result of long-term natural selection for a given site, and are also optimized in order to maximize the Net Carbon Profit. To do so, the water use strategy was expressed as a functional relation between the marginal water cost of assimilation (Cowan and Farquhar, 1977), represented by λ p and λ s for perennial and seasonal where c λf,s , c λe,s , c λf,p and c λe,p are the optimized parameters, while i r,p and i r,s represent the number of soil layers 110 reached by perennial and seasonal roots, respectively. After the establishement of the optimized water use parameters in Eqs. 1 and 2 (i.e. the long-term relation between soil water marginal water costs), the values of λ p and λ s are calculated for each day separately and then used to simulate the diurnal variation in stomatal conductance using Cowan-Farquhar optimality (Cowan and Farquhar, 1977;Schymanski et al., 2008a). The values of c λf,s , c λe,s , c λf,p and c λe,p essentially express how quickly plants reduce water use as soil water suction increases during dry periods. The slower λ s and λ p are reduced in response to 115 drying soil, the larger the root costs are as the root systems are adjusted to satisfy the canopy water demand. The parameters tions of the seasonal and perennial vegetation component, are allowed to vary on a daily basis to reflect their dynamic nature.
Their values are hence optimized from day to day in a way to maximize the daily NCP. This is done by using three different values for each of these vegetation properties, the actual value and a specific increment above and below this value every day, and at the end of the day the combination of values that would have achieved the maximum NCP on the present day is selected for the next day. See Schymanski et al. (2009Schymanski et al. ( , 2015 for details.

Carbon cost functions
As mentioned above, different carbon cost functions are used to quantify the maintenance costs for different plant organs. The carbon cost related to foliage maintenance is based on a linear relation between the total leaf area and a constant leaf turnover cost factor: where L AIc is the clumped leaf area (set to 2.5 (Schymanski et al., 2007)), c tc is the leaf turnover cost factor (set to 0.22 (Schymanski et al., 2007)) and M A,p is the perennial vegetation cover fraction.
The costs for root maintenance were defined as (Schymanski et al., 2008b): where c Rr is the respiration rate per fine root volume (0.0017 mol s −1 m −3 ), r r the root radius (set to 0.3*10 −3 m) and S A,r 135 the root surface area per unit ground area (m 2 m −2 ). The values used in this parameterization stem from observations on citrus plants, as described by Schymanski et al. (2008b).
Water transport costs are assumed to depend on the size of the transport system, from fine roots to the leaves. The canopy height is not modelled in the VOM, and the transport costs are therefore just a function of rooting depth and vegetated cover: where c rv is the cost factor for water transport (mol m −3 s −1 ), M A the fraction of vegetation cover (−), and y r the rooting depth (m). The costfactor c rv was set to 1.0 µmol m −3 s −1 by Schymanski et al. (2015) after a sensitivity analysis for Howard Springs, which is also adopted here. a saturated zone, overlaying an impermeable bedrock with a prescribed drainage level. The model simulates a variable water table based on the vertical fluxes between horizontal soil layers and a drainage flux computed as a function of the water table elevation. Here, the thickness of soil layers was prescribed to 0.2 m, based on a sensitivity analysis for Howard Springs, see also Supplement S2.
The hydrological parameters that determine the drainage outflow and groundwater tables are a hydrological length scale for 150 seepage outflow, channel slope and drainage level z r , as defined by Schymanski et al. (2015) and based on Reggiani et al. (2000). The seepage outflow is determined by the elevation difference between groundwater table and drainage level, divided by a resistance term that uses the hydrological length scale and channel slope (Eq. 10, Schymanski et al., 2008b). Originally, the hydrological length scale and channel slope were adopted from Reggiani et al. (2000) and set to 10 m and 0.033 rad, respectively, in absence of more detailed knowledge about these parameters. At the same time, Schymanski et al. (2015) set 155 the drainage level z r and total soil thickness c z to 10 m and 15 m, respectively, based on the local topography around the flux tower site (Schymanski et al., 2008b). This hydrological schematization resulted in groundwater tables around 5 m below the surface.
Here, the hydrological parameters were set in a way to resemble freely draining conditions, i.e. avoiding a significant influence of groundwater in Figure 2, for consistency with other model applications (e.g. Whitley et al., 2016), with a total soil 160 thickness c z of 30 m, a fast drainage parameterization with a drainage level z r of 5 m, a length scale for seepage outflow set to 2 m and a channel slope set to 0.02 rad. As illustrated in Figure 2, when precipitation falls on this soil block, it either causes immediate surface runoff or infiltrates. Once infiltrated, it can be taken up by roots and transpired, or it can evaporate at the soil surface, or move downwards until it drains away at a depth of 30 m, well below the rooting zone (i.e. parameterized to represent freely draining conditions for comparison with other models). The simulation of soil evaporation and vertical fluxes 165 in the unsaturated zone are described in Schymanski et al. (2008bSchymanski et al. ( , 2015.

Model optimization
The VOM uses the Shuffled Complex Evolution algorithm (SCE, Duan et al., 1994) to optimize the vegetation properties listed in Table 3 for maximum Net Carbon Profit (NCP) over the entire simulation period (37 years for the new VOM set-up, from 1-1-1980 until 31-12-2017). The SCE-algorithm uses first an initial random seed, subdivides the parameter sets into complexes 170 and performs a combination of local optimization within each complex and mixing between complexes to converge to a global optimum. Here, we set the initial number of complexes to 10.

Model input and data
A relatively long timeseries of meteorological inputs is required to run and optimize the VOM. the new VOM set-up. The meteorological data includes time series of daily maximum and minimum temperatures, shortwave radiation, precipitation, vapour pressure and atmospheric pressure. Atmospheric CO 2 -levels were originally assumed constant by Schymanski et al. (2015), but in the new set-up, these were taken from the Mauna Loa CO 2 -records (Keeling et al., 2005).
Observed atmospheric CO 2 -levels at the flux tower were not used due to the required length of the timeseries for the VOM (20-30 years). The measured meteorological variables at the flux tower sites were only used to verify the SILO meteorological 180 data, which revealed only minor differences in the resulting fluxes of the VOM when the SILO-data was replaced for the days that flux tower observations were available (max. 6%, see Figure S4.3 in Supplement S4). See also Supplement S3, Figure S3.1 for the time series of meteorological data.
The soils, originally assumed vertically homogeneous by Schymanski et al. (2015), from Carsel and Parrish (1988). See also Table 2 for the soil parameterization.
At Howard Springs, a flux tower that is part of the regional FLUXNET network OzFlux (Beringer et al., 2016), provides  (Asrar et al., 1984;Lu, 2003) allowed for the calculation of FPC by dividing the fPAR-values by the maximum value of 0.95.

Modifications to the VOM set-up
In comparison with previous applications of the VOM (Schymanski et al., 2009(Schymanski et al., , 2015, several changes were made regarding the input data and process representation. Each individual change was added to the reference set-up of Schymanski et al. (2015) 205 to assess the sensitivity of the model results for that change (see also Supplement S1). Briefly, the changes were assessed by the following cases:

210
-Re-run SCE: the VOM was re-optimized with the same settings and input data as (Schymanski et al., 2015), in order to assess wheter the optimization algorithm converges to the same results as (Schymanski et al., 2015). The model was run -Optimized grass rooting depth: the rooting depth of grasses was prescribed at 1.0 m by Schymanski et al. (2009Schymanski et al. ( , 2015, which is roughly the position of a hard pan in the soil profile at Howard Springs. In this study, grass rooting depth is -Modified soil properties and hydrology: the modified soils and hydrology, as described above, will strongly interact. For that reason, their combined effect was assessed by implementing both changes in the set-up of Schymanski et al. (2015), keeping everything else constant.

Effects of modifications to the VOM
To compare previous simulations using the VOM (Schymanski et al., 2015) with the new VOM set-up that includes the modifications as outlined in Sect. 2.2.8, each modification was applied to the previous setup in a stepwise manner to quantify the influence of each change in isolation. The resulting simulations were compared with those presented in Schymanski et al. 255 (2015) for the site Howard Springs. In general, sensitivities varied between +20% and -25% in total GPP and ET, and are summarized in Figure 3. See also Supplement S1 for detailed time series.  The updated meteorological input data, for the runs until 31-12-2005 as well as the extended runs until 31-12-2017, hardly influenced the outcomes, with less than 10% relative change in the resulting fluxes (Figure 3a and b). However, a higher contribution of the perennial vegetation in the fluxes can be observed, related to an increase in perennial vegetation cover 260 ( Figure 3i). This happened as well for re-running the SCE-algorithm, pointing at a relatively large uncertainty in the predicted perennial cover, with a range of values that results in similar fluxes.
Changing the fixed atmospheric CO 2 -levels (350 ppm) in the set-up of Schymanski et al. (2015) to variable atmospheric CO 2 -levels had a relatively large influence on perennial vegetation, yielding values of GPP for perennial vegetation that were up to 21.0% higher (Figure 3d). Note that the CO 2 -levels of the Mauna Loa records have a mean of 369 ppm and a maximum 265 of 410 ppm during the modelling period, i.e. mostly higher than the 350 ppm prescribed in Schymanski et al. (2015), who also simulated an increase in GPP by perennial plants in response to elevated CO 2 . See Figure S3.1f in Supplement S3 for more details about the CO 2 -levels used here.
Changing the vertical soil discretization of 0.5 m in Schymanski et al. (2015) to a finer resolution of 0.2 m had a minor influence, with a change of 2.6 % in the resulting GPP and 0.3 % in ET (Figure 3b, a). Similarly, when the grass rooting depths 270 were optimized instead of the prescribed grass rooting depth of 1 m (everything else being the same as in Schymanski et al. (2015)), simulated GPP and ET were changed by 1.0% and -2.3% respectively (Figure 3b and a). The optimization led to shallower grass roots of 0.5 m (incurring lower carbon costs) and therefore to reductions in GPP and ET.
transpiration by seasonal vegetation. At the same time, simulated soil evaporation was increased, relating to an increased soil water storage and pointing at a reduced ability of the roots to take up water due to reduced hydraulic conductivity in the soil. The increase in simulated GPP was largely due to increased GPP by perennial vegetation, which at the same time slightly increased its transpiration. These changes were connected to a largely increased perennial vegetation cover and reduced rooting depth compared to the original simulations (vegetation cover went up to 0.51 from 0.31, while rooting depth went down to 3.5 280 m from 4 m). Overall, the perennial vegetation benefited from the finer soil texture due to larger soil moisture storage capacity and carry-over of soil moisture into the dry season, whereas the seasonal vegetation suffered from reduced root water uptake due to lower hydraulic conductivity and increased soil evaporation during the wet season.
The implementation of free draining conditions had strong effects on the simulated fluxes as well, with lower values of both ET and GPP (-20.4% and -6.9% respectively, Figure 3a and b). However, here especially the simulated ET of perennial 285 vegetation was reduced, whereas the transpiration by seasonal vegetation stayed relatively similar (Figure 3c and e). This is because in the original simulations, capillary rise from the water table was most important during the dry season, when seasonal vegetation is inactive, and a change in the water table due to free draining conditions affects therefore mostly the perennial and not so much the seasonal vegetation.
Combining the new soils with the new hydrological settings still resulted in a reduction in ET by -11.5%, whereas their 290 combined effect on GPP led to only a small reduction by -1.1% (Figure 3a and b). Here, the reduction of ET occured mainly during the wet season, and related to reductions in the perennial transpiration (Figure 4a-c), whereas the GPP stayed relatively similar ( Figure 4e).
These findings are in accordance with the isolated effects of the new soils and hydrology, where free drainage conditions resulted in a large reduction in ET and GPP, while finer soil texture resulted in a small reduction in ET but large increase in 295 GPP. Hence, the finer soil texture largely compensated the effect on GPP, but not on ET.
The other changes had smaller effects on the fluxes, so we do not discuss them here individually. Instead, we will perform a more in-depth analysis of the differences in model results when all changes were combined.

Resulting differences
After incorporating all the changes, the relative error for mean annual evapo-transpiration (ET) changed from an overestimation 300 by 8.4% to an underestimation of -10.2%, whereas the relative error for the mean annual GPP changed from 17.8% to 14.7%.
The ensemble years in Figure 5 revealed that the evapo-transpiration (ET) was most strongly underestimated by the VOM during the dry season at Howard Springs. The observed groundwater tables (Figure 6a) ranged from 5-15 m depth seasonally, whereas the VOM was parameterized now to keep groundwater tables close to 25 m depth, for consistence with free drainage conditions in other models. Schymanski et al. (2015) originally assumed a much shallower drainage level at Howard Springs, 305 which led to groundwater tables around 5 meters depth, and better correspondence with the observed fluxes (Fig. 5). The simulated soil moisture in the top soil layer, as illustrated in Fig. 6b, remained similar to the soil moisture values of Schymanski et al. (2015). The higher vertical resolution in the new model runs (20 cm cf. 50 cm soil layers) resulted in stronger surface soil moisture spikes around rainfall events, which makes the red line appear generally more noisy than the green line in Fig. 6b. Observed soil moisture in the upper 5 cm was generally lower than the simulated soil moisture in the top soil layer,

Conclusions
The Vegetation Optimality Model has undergone several changes regarding model set-up and input data since its last application by Schymanski et al. (2015). The modifications consisted of updated and extended input data, the use of variable atmospheric CO 2 -levels, modified soil properties, modified drainage levels as well as the addition of grass rooting depths to the optimized vegetation properties. The changes were applied to the VOM in a step-wise manner, by applying each modification to the 320 previous set-up of Schymanski et al. (2015) in isolation to evaluate its effect on the results, before combining all modifications and analyzing their effect in combination. This analysis revealed that updated soil textures and a changed hydrological schematization had a strong influence on the results. An underestimation of dry season ET at Howard Springs was much more apparent when compared to the results of Schymanski et al. (2015), where the drainage parameterization maintained a water simulations was partly buffered by a more fine-grained soil texture below 0.4 m (sandy clay loam instead of sandy loam), which resulted in an increase in water storage at otherwise similar water potentials in the top 5 m of soil compared to the simulations by Schymanski et al. (2015). The use of variable atmospheric CO 2 -levels also had a strong influence on the results, which is especially important as the model time period has been extended in this study. This was mainly due to generally higher levels of atmospheric CO 2 in recent observations, compared to the constant values used by Schymanski et al. (2015).

330
The stepwise analysis of modifications to a model setup and comparison against a benchmark dataset proved very helpful for identifying sensitivities of simulations to the different changes that might otherwise remain undiscovered due to compensating effects of the various modifications during model development. In this way, we found that the neglect of a varying water table may have a strong effect on simulated surface fluxes, especially when soils are highly permeable. The common assumption of free draining conditions in modelling studies should be revised, and if such an assumption is necessary, due to a lack of better 335 hydrological understanding of a given site, or for comparison with other model simulations using this assumption (e.g. Whitley et al., 2016), potential bias in simulation results has to be acknowledged.