Influence of modifications (from AoB2015 to v0.5) in 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. Several updates to previous applications of the VOM have been made for the study in the accompanying paper of Nijzink et al. (2021), where we assess whether optimality 5 theory can alleviate common shortcomings of conventional models, as identified in a previous model inter-comparison study along the North Australian Tropical Transect (NATT) (Whitley et al., 2016). Therefore, we assess in this technical paper how the updates to the model and input data would have affected the original results :: of ::::::::::::::::::::: Schymanski et al. (2015), and we implemented them ::::: these :::::: changes : one at a timewhile 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 :: in :::::::::: comparison :::: with :: the :::::: results ::: of :::::::::::::::::::: Schymanski et al. (2015), whereas more fine-grained 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 20 from an 8.4% overestimation to an 10.2% underestimation, whereas the relative errors for the mean annual GPP stayed similar with a change :::::::: remained :::::: similar :::: with ::: an :::::::::::: overestimation :::: that ::::::: slightly ::::::: reduced : 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. 25

common TBMs, in the accompanying paper by Nijzink et al. (2021) the VOM was applied to the same sites along the NATT and systematically compared with the previous simulations presented by Whitley et al. (2016). 60 In order to understand if predicted optimality-based rooting depths and vegetation cover result in better simulations, Nijzink et al. (2021) ran the VOM both with predicted and prescribed rooting depths and vegetation cover while systematically comparing simulated fluxes with observations and the output of the other TBMs. For that reason, Nijzink et al. (2021) made several changes to the VOM set-up of Schymanski et al. (2015) in order to use the same input data and similar physical boundary conditions at the different sites as the TBMs in Whitley et al. (2016). In the remainder of this paper, the new set-up in Nijzink et al. (2021) will be referred to as , in contrast to VOM-AoB2015 for the set-up of Schymanski et al. (2015).
First, in the simulations by Whitley et al. (2016), all TBMs were run under the assumption of a freely draining soil column.
In contrast, the VOM-AoB2015 used a hydrological schematization based on the local topography around the flux tower site (Schymanski et al., 2008b), which resulted in groundwater tables varying around 5 m below the surface. For better comparability with Whitley et al. (2016), the boundary conditions of the VOM were adjusted to resemble freely draining conditions. 70 Assessment of the influence of this change could lead to additional insights about the influence of groundwater on the resulting carbon and water fluxes, which can be significant (York et al., 2002;Bierkens and van den Hurk, 2007;Maxwell et al., 2007).
Another modification concerns the prescribed atmospheric CO 2 -concentrations. The VOM-AoB2015 assumed constant atmospheric CO 2 -concentrations over the entire modelling period, whereas the VOM-v0.5 used measured CO 2 -concentrations, which have increased considerably over the past years (Keeling et al., 2005). The previously documented influence of at-75 mospheric CO 2 -concentrations on the water and carbon fluxes simulated by the VOM (Schymanski et al., 2015) calls for a systematic assessment of this change.
In the VOM-AoB2015, the grass rooting depth was prescribed to a value of 1 m, arguing that only tree roots could penetrate into deeper layers due to the presence of a hard pan at Howard Springs, which not necessarily exists at the other sites along the NATT. In order to make the VOM applicable to all sites along the NATT, the grass rooting depth was not prescribed but also 80 optimized for NCP in the VOM-v0.5, but an evaluation of the effect of this change on the original simulations was not included in Nijzink et al. (2021).
In addition to the changes in the boundary conditions, model code and parametrization mentioned above, higher computational power, allowing for finer discretizations, and updated forcing data may also affect simulation results. In general, reproducing benchmark datasets is often seen as necessary in order to be confident about the model and the numerical imple-85 mentation (Blyth et al., 2011;Abramowitz, 2012;Clark et al., 2021). Here we argue that it is also important to assess effects of individual changes one at a time, to avoid obtaining "the right results for the wrong reasons" as two errors may compensate for each other when comparing the final results with a given benchmark. Therefore, we assess to what extent the various changes influence the VOM-results, by adding these changes one at a time to the VOM-AoB2015, for the same flux tower site in Australia, Howard Springs, as in Schymanski et al. (2009Schymanski et al. ( , 2015. This 90 technical note describes the nine changes to the VOM since its last application by Schymanski et al. (2015), and how they affect the results of the VOM individually and in combination. In this way, this work should point at important model decisions and sensitivities of the VOM, as well as TBMs in general. At the same time, this works :::: work showcases a systematic evaluation of model updates and changes, which we deem necessary in model applications.

95
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 : ,3 . 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 used by Schymanski et al. (2009) and Schymanski et al. (2015) is Howard Springs (AU-How), which was 100 therefore used in this analysis as well. At the same time, the flux tower site at Howard Springs provides a long record of carbon dioxide and water fluxes starting from 2001 (Beringer et al., 2016). Howard Springs is the wettest site with an average precipitation of 1747 mm/year (SILO Data Drill, Jeffrey et al., 2001Jeffrey et al., , calculated for 1980Jeffrey et al., -2017 along the North Australian Tropical Transect (NATT, Hutley et al., 2011), which has a strong precipitation gradient from north 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 105 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.

155
, :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: (values taken for Eucalyptus pauciflora, see Schymanski et al., 2007;Medlyn et al., 2002) : . ::::::::::: Furthermore, ::::::: J max,25 :: is the electron transport capacity at 25 o C J max,25 (mol m −2 s −1 ) and :::: T opt the optimal temperature T opt (K), set to the mean monthly daytime temperature at the site (305 K, Schymanski et al., 2007). In the equations above, the electron transport capacity at 25 o C J max,25 , the projected foliage cover M a and stomatal conductance G s are optimized dynamically in a way to maximize the overall NCP of the vegetation over the entire simulation period. Optimization is possible due to the carbon costs associated 160 with each of these variables: photosynthetic capacity is linked to maintenance respiration, projected cover is linked to foliage turnover and maintenance costs, 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 :::: (Q r,i , ::::: m/s), is modelled following an electrical circuit analogy, where the water potential difference between the plant and each soil layer drives the flow: 165 with S A,r the root surface area (m 2 m −2 ), h r,i the hydraulic head in the roots (m), h i the hydraulic head in the soil (m), Ω r the radial root resistivity (s) and Ω s,i the soil resistivity (s), with subscript i denoting the specific soil layer. The root surface area, S A,r , is optimized in a way to satisfy the canopy water demand with the minimum possible total root surface area.

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 the perennial vegetation (M A,p ) are derived by optimizing these properties for the long-term, assuming that these do not vary 195 significantly during the simulation period (20-30 years). Similarly, water use strategies of both the perennial and the seasonal vegetation components are assumed to be a result of long-term natural selection for a given site, and are also optimized in order to maximize the NCP. 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 (mol mol −1 m −1 ) for perennial and seasonal vegetation respectively, and the sum of water suction heads (h i ) in all soil layers within the root zone: where c λf,s (mol mol −1 m −1 ), c λe,s (-), c λf,p (mol mol −1 m −1 ) and c λe,p (-) are the optimized parameters, while i r,p and i r,s represent the number of soil layers reached by perennial and seasonal roots, respectively. Note that Cowan and Farquhar (1977) proposed that λ should decline with declining soil water content, whereas Schymanski et al. (2009) argued that plants 205 more likely sense the soil suction head than the total available water. Equations 11 and 12 formulated λ as an explicit but flexible function of the average suction head in the root zone, where the shape of the function (determined by the two optimized parameters) represents a specific water use strategy.
After the establishement of the optimized water use parameters in Eqs. 11 and 12, 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 210 (Cowan and Farquhar, 1977;Schymanski et al., 2008a). The values of c λf,s , c λe,s , c λf,p and c λe,p express how quickly plants reduce water use as soil water suction increases during dry periods. The parameters (c λf,s , c λe,s , c λf,p and c λe,p ) are optimized and constant in the long-term, along with y r and M A,p , to maximize the total NCP over the entire simulation period.

Short-term optimization
Some vegetation properties, such as seasonal vegetation cover (M A,s ), the electron transport capacities at 25 o C for the seasonal 215 and perennial vegetation (J max25,p and J max25,p ) and root surface area distributions of the seasonal and perennial vegetation (S Adr,i,s and S Adr,i,p ), are allowed to vary on a daily basis to reflect their dynamic nature. Here, the root surface area distributions are optimized day by day in a way to satisfy the canopy water demand. :: In : a :::: first :::: step, :::: this :: is :::: done ::: by :::::::::: determining :: a :::::::: coefficient ::: of :::::: change ::: k r , :::::: defined ::: as:

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 240 cost factor: where L AIc is the clumped leaf area (set to 2.5), c tc is the leaf turnover cost factor (set to 0.22 µmol −1 s −1 m −2 ) and M A,p is the perennial vegetation cover fraction.
The costs for root maintenance were defined as: 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), which 245 were both derived by Schymanski et al. (2008b) from experimental data on citrus plants. Here, we present a sensitivity analysis of these parameters in Supplement S5. S A,r represents the root surface area per unit ground area (m 2 m −2 ).
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 Springs, which is also adopted here.

Water balance model
The soil is schematized as a permeable block containing an unsaturated zone and a saturated zone (see also  The hydrological parameters that determine the drainage outflow and groundwater tables are a hydrological length scale for 260 seepage outflow, channel slope and drainage level z r , 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): with γ 0 the average slope of the seepage face (rad), Λ s a hydrological length scale (m), y s the groundwater table (m), K sat the saturated hydraulic conductivity (m s −1 ), ω 0 the saturated surface area fraction and z r the drainage level (m).

265
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 defined by the drainage level z r and total soil thickness c z . In the top soil layer, soil evaporation is determined as a function of the soil moisture, global radiation and vegetation cover.

Model optimization 270
The VOM uses the Shuffled Complex Evolution algorithm (SCE, Duan et al., 1994) to optimize the vegetation properties listed in Table 3 for maximum NCP over the entire simulation period (37 years for the VOM-v0.5, from 1-1-1980 until 31-12-2017).
The SCE-algorithm uses first an initial random seed, subdivides the parameter sets into complexes 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. The VOM uses a variable time step, where the target time step length of 1 hour is reduced 275 if any state variable in the model changes by more than 10% per time step.

Meteorological data
A relatively long timeseries of meteorological inputs is required to run and optimize the VOM. The necessary meteorological data includes daily time series of maximum and minimum temperatures, shortwave radiation, precipitation, vapour pressure and atmospheric pressure, which were taken from the Australian SILO Data Drill . In addition, the VOM 280 requires information about atmospheric CO 2 -levels, which were provided as a constant in the VOM-AoB2015 version, whereas for the VOM-v0.5 used here, we used the Mauna Loa CO 2 -records (Keeling et al., 2005), see also Sect. 2.2.9. Eventually, these daily time series are converted in the VOM to hourly time series, as described by ::: with :: a :::::: diurnal :::::::: variation ::: that :::: was ::::::: imposed ::: for :: the :::::: global :::::::: radiation, :::::::::: temperature ::: and :::::: vapour ::::::: pressure :::::: deficit, :: as ::::::::: described :: in ::::: detail :: in :::::::: Appendix :  Observed atmospheric CO 2 -levels at the flux tower were not used due to the required length of the timeseries for the VOM 285 (20-30 years). The measured meteorological variables at the flux tower sites were only used to verify the SILO meteorological 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.

290
At Howard Springs, a flux tower that is part of the regional FLUXNET network OzFlux (Beringer et al., 2016), provides time series of net ecosystem exchange (NEE) of carbon dioxide and latent heat flux (LE) for model evaluation. The Dingo-algorithm (Beringer et al., 2017) was applied to the data for a gap-filled estimation of gross primary productivity (GPP) and latent heat flux (LE). LE was converted to evapo-transpiration (ET), defined here as the sum of all evaporation and transpiration processes, even though these processes are different in nature (Savenije, 2004). Eventually, the gap-filled observations of GPP and ET  (Asrar et al., 1984;Lu, 2003) allowed for the calculation of FPC by dividing the fPARvalues by the maximum value of 0.95.
To assess the effect of the modifications to the VOM set-up from version VOM-AoB2015 to VOM-v05, each individual change was added to the reference set-up of the VOM-AoB2015 (see also Supplement S1). We define here 12 model cases, including the 9 changes to the VOM-AoB2015, the reproduction of the VOM-AoB2015, the re-optimization of the VOM-AoB2015 and the final VOM-v0.5 (see also Table 4): -1. Reproduction of the results of Schymanski et al. (2015) 310 The model code of the VOM-v0.5 was run with the same vegetation parameters and input data as the VOM-AoB2015.
The model was run from 1-1-1976 until 31-12-2005. This was done in order to check the reproducibility of the results of Schymanski et al. (2015).
The VOM was re-optimized with the same settings and input data as the VOM-AoB2015, and run with the same model 315 period from 1-1-1976 until 31-12-2005. Also here, the specific goal was to reproduce the results of Schymanski et al. (2015), as the optimization algorithm should converge to the same solutions.
-3. Variable CO 2 -levels Atmospheric CO 2 -levels were originally assumed constant in the VOM-AoB2015 with CO 2 -concentrations of 350 ppm, but in the VOM-v0.5, these were taken from the Mauna Loa CO 2 -records (Keeling et al., 2005). Therefore, the VOM-

320
AoB2015 was run with variable CO 2 -levels, and optimized for the period 1-1-1976 until 31-12-2005. This was done in order to assess the sensitivity of the model to variable CO 2 -levels.  -1-1980. This was performed to assess the importance of precise atmospheric pressure data for the VOM simulations.
The rooting depth of grasses was prescribed at 1.0 m in the VOM-AoB2015, which is roughly the position of a hard pan in the soil profile at Howard Springs. In the VOM-v0.5, grass rooting depth is optimized along with the tree rooting depth in order to let also the grass rooting depths adapt to local conditions. To assess the effect of an optimized grass rooting depth separately, we also optimized it in the VOM-AoB2015 simulations for the period 1-1-1976 until 31-12-2005. In the VOM-AoB2015, the average slope of the seepage face γ 0 and hydrological length scale Λ s were adopted from Reggiani et al. (2000) and set to 0.033 rad and 10 m, respectively, in absence of more detailed knowledge about these parameters. At the same time, the drainage level z r and total soil thickness c z were set to 10 m and 15 m, respectively, based on the local topography around the flux tower site (Schymanski et al., 2008b). This hydrological schematization 355 resulted in groundwater tables around 5 m below the surface.
The hydrological parameters for the VOM-v0.5 were set in a way to resemble freely draining conditions, i.e. avoiding a significant influence of groundwater in Figure  into one of the soil textural groups of Carsel and Parrish (1988) based on the fractions of sand, silt and clay. Eventually, the parameters for the soil water retention model of van Genuchten (1980) and the hydraulic conductivity were taken from the accompanying tables 7 from Carsel and Parrish (1988). See also Table 2 for the soil parameterization in the 370 VOM-v0.5. As a result, the soil profile at Howard Springs is now assumed to consist of sandy loam in the top 0.4 m and sandy clay loam below, whereas VOM-AoB2015 used a soil of sandy loam in the entire soil profile.
The VOM-AoB2015 was optimized here with the a modified soil profile, using the soil discretization of the VOM-AoB2015 of 0.5 m, with a sandy loam structure in the top layer and sandy clay loam below. This was done for the period 1-1-1976 until 31-12-2005, in order to assess the changes due to the different soils.

375
-11. Modified soil properties and hydrology The modified soils and hydrology, as described above, will strongly interact. Free draining conditions are expected to reduce soil water storage, while finer soil texture is expected to increase water storage. In order to better understand in how far these changes compensated for each other, they were both implemented together in the VOM-AoB2015, while keeping everything else unmodified.

385
To compare previous simulations using the VOM-AoB2015 with the VOM-v0.5 set-up that includes the modifications as outlined in Sect. 2.2.9, each modification was applied to the previous setup in a one-step-at-a-time approach to quantify the influence of each change in isolation. The resulting simulations were compared with those presented in Schymanski et al. (2015) for the site Howard Springs. In general, sensitivities varied between +20% and -25% in total GPP and ET after optimizing the VOM with the new changes, and are summarized in Figure 3. Without optimizing the VOM, but still including the modifica-390 tions, the sensitivities in GPP and ET were even smaller (see Supplement S1, Figure S1.51), except for a strong increase in soil evaporation after changing the soils. See also Supplement S1 for detailed time series.
Besides the changes in the fluxes, several changes in the optimized vegetation properties were observed as well ( Figure 4 and Figure S1.49 in Supplement S1) . Similar to the fluxes, the changes in vegetation parameters for reproducing VOM-AoB2015 and re-running the optimization algorithm remained small ( Figure 4).

465
Interestingly, the implementation of variable CO 2 -concentrations led to a large increase in c λe,p , i.e. one of the water use strategy parameters of the perennials (Fig. 4a). Another effect was a larger perennial vegetation cover (Figure 4e), while the seasonal cover was on average reduced, which relates to the generally elevated CO 2 -levels and the long-term optimization of the perennial vegetation cover. Hence, perennial vegetation cover benefited, and the grass cover, optimized on a daily basis, reduced on average (Figure 4f).

470
The reduced soil layer thickness mainly affected the water use parameters c λe,p and c λe,s (Figure 4a and b, respectively).
At the same time, the root depths for the perennials increased (Figure 4g) , which compensated for the change in water use, resulting in only minor changes in the fluxes ( Figure 3).
The variable atmospheric pressure led to changes in the vegetation parameters as well, but this stayed limited to a maximum of 25 % for c λe,p (Figure 4a). However, this is stronger than the observed changes in the resulting fluxes of ET and GPP ( Figure   475 3), which remained much smaller, related also to the non-linear relationship between c λ , ET and GPP.
Optimized grass rooting depths were strongly reduced compared to the prescribed 1 m, and this was accompanied by increased c λe,s and reduced c λf,s (Figure 4b and d), pointing at a more efficient water use strategy with less water transpired per assimilated CO 2 (Figure 3h).
The updated meteorology, as well as the updated and extended meteorology, led to changes in the vegetation properties as 480 well. The projected cover for the perennials increased (24.8 %, Figure 4) from 31.7 % to 39.6 %, whereas the projected cover for the grasses decreased on average (-6.7 % , Figure 4). Also the water use strategy parameters (4a-d) changed, with especially a strong change for c λf,p with -27.8 % (Figure 4c), but the resulting total water use efficiency remained similar (Figure 3h).
The modified hydrology led to larger changes in the vegetation properties, with especially a strong decrease in c λf,p ( Figure   4c), but as well for c λf,s (Figure 4d). Hence, the modified hydrology led to a more efficient water use (Figure 3h). At the same 485 time, the root depths for the perennial trees strongly reduced as well (Figure 4g). The modified soils, in isolation, led to similar effects for the water use strategy parameters (Figure 4a-d), which resulted in a more efficient water use (Figure 3h), see also Supplement S1, Figure S1.38d). The projective cover of the perennial trees changed strongly increased by 60 % (Figure 4g), i.
The ensemble years in Figure 6 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 7a) ranged from 5-15 m depth seasonally, 505 whereas the VOM was parameterized now to keep groundwater tables close to 25 m depth, required for the accompanying paper of Nijzink et al. (2021). Schymanski et al. (2015) originally assumed a much shallower drainage level at Howard Springs, which led to groundwater tables around 5 meters depth, and better correspondence with the observed fluxes (Fig. 6).

Conclusions
As models, input data and parameterizations evolve, the effects of individual improvements are usually not systematically 530 evaluated. Here, we analyse the effects of changes to the Vegetation Optimality Model (VOM) between version VOM-AoB2015 and VOM-v0.5, the latter of which is the basis of a companion paper (Nijzink et al., 2021)). Some of the modifications were done for improved realism, while others for better comparability with other models and general applicability across the different sites investigated.
The modifications consisted of updated and extended input data, the use of variable atmospheric CO 2 -levels, modified soil 535 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 one-step manner for the flux tower site Howard Springs, by applying each modification to the 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.   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).
Hence, our approach led to the insights that the neglect of a varying water table may have a strong effect on simulated surface fluxes, especially when soils are highly permeable. This is in line with other studies (e.g. York et al., 2002;Bierkens 550 and van den Hurk, 2007;Maxwell et al., 2007), and shows that the common assumption of free draining conditions in modelling studies should be revised. Interestingly, the deficiencies in TBMs related to water access and tree rooting depth as identified by Whitley et al. (2016) strongly relate to this, as root depths also depend on groundwater tables. At least, if a free draining assumption is necessary, due to a lack of better hydrological understanding of a given site, or for comparison with other model simulations, i.e. for the study in the accompanying paper of Nijzink et al. (2021), potential bias in simulation results has to be 555 acknowledged.
In addition to this, we can more generally conclude that optimality-based modelling is able to provide robust modelling results. The sensitivities for the different changes remained rather limited and varied only between +20% and -25% in total GPP and ET after re-optimizing, but implementing the changes without re-optimizing the vegetation properties, resulted in even smaller changes. Therefore, we gained confidence that the VOM will also provide reliable and robust results for other 560 study sites along the North Australian Tropical Transect, with similar boundary conditions. At the same time, the new boundary conditions made the VOM comparable to the other TBMs of Whitley et al. (2016).
To conclude, with our analysis we identified more conceptual issues, e.g. the influence of the hydrological schematization, but also found confirmation that the optimality-based modelling approach provides a robust result. Our method is in line with other developments in terrestrial biosphere modelling, where benchmark testing is seen as a necessary step in the modelling 565 process (Blyth et al., 2011;Abramowitz, 2012;Clark et al., 2021). Here, we provided an example on how to perform a systematic benchmark test in a one-step-at-a-time approach, i.e. applying one change in isolation to the benchmark model. This 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. Our work also highlights the importance of open source code and data.

570
The availability of the original VOM-AoB2015 code was a pre-requisite of conducting this analysis. The long-term storage of the new VOM-v0.5, all input data, parameterization and data analysis workflows used in the present study will ensure that the effects of future modifications can again be compared step-by-step to the results presented here. This is a pre-requisite for maintaining generality of a model, as opposed to models that are highly customized for a particular site and hence unable to provide general insights.
Author contributions. SJS and RN designed the set-up of the study. Model code was originally developed by SJS, but updated and modified by RN. RN set-up the repositories for the pre-processing, modelling and post-processing on renkulab. LH and JB provided site-specific 580 knowledge and data. The main manuscript and supplementary information were prepared by RN, together with input from SJS. LH and JB provided corrections, suggestions and textual inputs for the main manuscript.