Articles | Volume 15, issue 2
Development and technical paper
01 Feb 2022
Development and technical paper |  | 01 Feb 2022

Influence of modifications (from AoB2015 to v0.5) in the Vegetation Optimality Model

Remko C. Nijzink, Jason Beringer, Lindsay B. Hutley, and Stanislaus J. Schymanski

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. (2022), where we assess whether optimality 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 these changes one at a time.

The model updates included extended input data, the use of variable atmospheric CO2 levels, modified soil properties, 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 evapotranspiration (ET) and gross primary productivity (GPP) by no more than 20 %, with the largest effects caused by the newly imposed free drainage 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 from an 8.4 % overestimation to an 10.2 % underestimation, whereas the relative errors for the mean annual GPP 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.

1 Introduction

Novel modelling approaches that are able to explicitly model vegetation dynamics, such as vegetation cover or root surfaces, may lead to an overall improved understanding of carbon and water flux exchanges with the atmosphere. At the same time, terrestrial biosphere models (TBMs) often rely on (remotely sensed) observations of vegetation properties, such as the leaf area index (LAI) or vegetation cover, complicating the ability of these models to make predictions for future scenarios. In addition, this makes the models rely on the quality of the data and does not enhance our understanding of the vegetation dynamics, which is highly important regarding the feedbacks between the land and the atmosphere in a changing climate. Recent model inter-comparison studies also confirm that models with explicit vegetation dynamics are needed, as Whitley et al. (2016) showed that prescribing rooting depths and the lack of dynamic representations of LAI led to highly variable performances for a selection of TBMs. Optimality theory predicts the variation and dynamics of vegetation cover, root systems, water use and carbon uptake without the need for site-specific input about vegetation properties by optimizing these properties for a certain objective, such as maximizing the carbon gain by photosynthesis (Hikosaka2003; Raupach2005; Buckley and Roberts2006) or minimizing water stress (Rodríguez-Iturbe et al.1999a, b). The theory used here is based on the premise that the net carbon profit (NCP, Schymanski et al.2007, 2008a, 2009), which is the difference between carbon assimilated by photosynthesis and carbon expended on construction and maintenance of all the plant tissues needed for photosynthesis, 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 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.2009, 2015). The VOM is a coupled water–vegetation model that optimizes vegetation properties to maximize the NCP in the long term (20–30 years) for given climate and physical properties at the site under consideration. 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 circumstances. The NATT has been used previously in an inter-comparison 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. Previous studies show that rooting depths vary considerably with precipitation (Schenk and Jackson2002) and, thus, also along the North Australian Tropical Transect (Williams et al.1996; Ma et al.2013). 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. In order to find out whether this novel approach can help to overcome the shortcomings of common TBMs, in the accompanying paper by Nijzink et al. (2022), the VOM was applied to the same sites along the NATT and systematically compared with the previous simulations presented by Whitley et al. (2016).

In order to understand whether predicted optimality-based rooting depths and vegetation cover result in better simulations, Nijzink et al. (2022) 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. (2022) 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. (2022) will be referred to as VOM-v0.5, 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. Assessment of the influence of this change could lead to additional insights into the influence of groundwater on the resulting carbon and water fluxes, which can be significant (York et al.2002; Bierkens and van den Hurk2007; Maxwell et al.2007).

Another modification concerns the prescribed atmospheric CO2 concentrations. The VOM-AoB2015 assumed constant atmospheric CO2 concentrations over the entire modelling period, whereas the VOM-v0.5 used measured CO2 concentrations, which have increased considerably over the past years (Keeling et al.2005). The previously documented influence of atmospheric CO2 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 does not necessarily exist 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 was also optimized for NCP in the VOM-v0.5.

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 implementation (Blyth et al.2011; Abramowitz2012; 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. (2009, 2015). This 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 to important model decisions and sensitivities of the VOM as well as TBMs in general. At the same time, this work showcases a systematic evaluation of model updates and changes, which we deem necessary in model applications.

2 Methodology

All steps in the process, from pre- and post-processing to model runs, were done in an open science approach using the RENKU1 platform. The workflows, including code and input data, can be found online2. 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).

2.1 Study site

The study site used by Schymanski et al. (2009) and Schymanski et al. (2015) is Howard Springs (AU-How), which was 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 yr−1 (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 to south, with a mean annual precipitation around 500 mm yr−1 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.

Figure 1Location of the Howard Springs site, together with the other flux tower sites that are part of the North Australian Tropical Transect in the Northern Territory of Australia, with the mean annual precipitation shown on the blue colour scale (SILO Data Drill, Jeffrey et al.2001, calculated for 1980–2017).

Table 1Characteristics of the Howard Springs site. Vegetation data from Hutley et al. (2011) and Whitley et al. (2016), with Eucalyptus (Eu.), Erythrophleum (Er.), and Hetropogan (He.). Meteorological data are taken from the SILO Data Drill (Jeffrey et al.2001) for the model periods of 1 January 1980 until 31 December 2017, with the reference crop evapotranspiration calculated according to the FAO Penman–Monteith formula (Allen et al.1998). The ratio of the net radiation Rn to the latent heat of vaporization λ multiplied by 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.

Download Print Version | Download XLSX

2.2 Vegetation Optimality Model

The Vegetation Optimality Model (VOM, Schymanski et al.2009, 2015) is a coupled water and vegetation model that optimizes vegetation properties by maximizing the NCP. The model code and documentation can be found online3, and version v0.54 of the model was used here. The processes and parameterizations are not modified in VOM-v0.5 unless explicitly stated in Sect. 2.2.9. Therefore, more details of the VOM can be found in Schymanski et al. (2009, 2015), whereas detailed descriptions about the root processes can be found in Schymanski et al. (2008b) and the canopy processes in Schymanski et al. (2007). Nevertheless, a general description of the VOM is given below for completeness.

Figure 2Schematization of the Vegetation Optimality Model as two big leaves, with MA,p and MA,s the fractional cover of perennial trees and seasonal grasses respectively, yr,p and yr,s the rooting depths of the perennial trees and seasonal grasses respectively, ΔZ the soil layer thickness, CZ the total soil depth, and Zr the drainage depth.


2.2.1 Water balance model

The soil is schematized as a permeable block containing an unsaturated zone and a saturated zone (see Fig. 2) 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. The vertical fluxes between soil layers are determined using a discretization of the Buckingham–Darcy equation (Radcliffe and Rasmussen2002), resulting in the 1-D Richards' equation of steady flow. The matric suction heads and unsaturated conductivities were determined with the model of van Genuchten (1980).

The hydrological parameters that determine the drainage outflow and groundwater tables are a hydrological length scale for seepage outflow, channel slope and drainage level zr, 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):

(1) Q s , f = K sat ω 0 ( y s - z r ) 2 cos ( γ 0 ) Λ s ,

with γ0 the average slope of the seepage face (rad), Λs a hydrological length scale (m), ys the groundwater table (m), Ksat the saturated hydraulic conductivity (m s−1), ω0 the saturated surface area fraction and zr the drainage level (m).

As illustrated in Fig. 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 zr and total soil thickness cz. In the top soil layer, soil evaporation is determined as a function of the soil moisture, global radiation and vegetation cover (Eqs. 18, 19, Sect. 3.3.3, Schymanski et al.2009).

2.2.2 Vegetation model

The VOM schematizes the ecosystem as two big leaves (see Fig. 2), one representing the seasonal vegetation (grasses) and one representing the perennial vegetation (trees). Currently, the VOM does not explicitly consider leaf area dynamics, but the photosynthesis of both the seasonal and perennial vegetation was modelled with a simplified canopy–gas exchange model of von Caemmerer (2000) for C3 plants. This was done by Schymanski et al. (2009) for a greater generality of the VOM, even though it may not correctly represent photosynthesis of the C4 grasses at the site. The model computes CO2 uptake as a function of irradiance, atmospheric CO2 concentrations, photosynthetic capacity and stomatal conductance:

(2) A g = 1 8 4 C a G s + 8 Γ * G s + J e - 4 R l - 4 G s C a - 2 Γ * 2 + 16 G s 8 C a G s + J e + 8 R l Γ * 1 2 ,

with Je the electron transport rate (mol m−2 s−1), Gs stomatal conductance (mol m−2 s−1), Rl leaf respiration (mol m−2 s−1), Ca the mole fraction of CO2 in the air and Γ* the CO2 compensation point (mol CO2 mol−1 air). The electron transport rate Je was calculated as

(3) J e = 1.0 - e 0.3 I a J max J max M A ,

with Ia the irradiance (mol m−2 s−1), Jmax the electron transport capacity (mol m−2 s−1) and MA the projected cover of vegetation (dimensionless fraction). The leaf respiration Rl, as used in Eq. (2), is defined as

(4) R l = M A c R l J max C a - Γ * 8 C a + 2 Γ * ,

with cRl a constant set to 0.07 (dimensionless), as defined by Schymanski et al. (2007) and based on the results of Givnish (1988), who found that leaf respiration equals 7 % of the maximum photosynthetic capacity for a range of different species. The electron transport capacity Jmax in Eqs. (3) and (4) is determined in the following way:

(5) J max = J max , 25 h a e - h d T opt - 298.0 273 T opt R + 25.0 - 1.0 + h d e h a T a - 25.0 273.0 T opt R + T opt - 273.0 T a + 273.0 T opt R 273.0 T opt R + 25.0 h a e h d T a - T opt + 273.0 T a + 273.0 T opt R - 1.0 + h d ,

with ha the rate of exponential increase in the function below the Jmax,25 and hd the rate of exponential decrease in the function above Jmax,25, set to 43.79 and 200 kJ mol−1, respectively (values taken for Eucalyptus pauciflora; see Schymanski et al.2007; Medlyn et al.2002). Furthermore, Jmax,25 is the electron transport capacity at 25 C (mol m−2 s−1) and Topt the optimal temperature (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 C Jmax,25, the projected foliage cover MA and stomatal conductance Gs 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 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 (Qr,i, m s−1) is modelled following an electrical circuit analogy, where the water potential difference between the plant and each soil layer drives the flow:

(6) Q r , i = S A , r h r , i - h i Ω r + Ω s , i ,

with SA,r the root surface area (m2 m−2), hr,i the hydraulic head in the roots (m), hi 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, SA,r, is optimized in a way to satisfy the canopy water demand with the minimum possible total root surface area.

2.2.3 Carbon cost functions and net carbon profit

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 (Rf) is based on a linear relation between the total leaf area and a constant leaf turnover cost factor:

(7) R f = L AIc c tc M A , p ,

where LAIc is the clumped leaf area index (LAI of vegetated area, set to a constant 2.5 based on Schymanski et al.2007), ctc is the leaf turnover cost factor (set to 0.22 µmol-1s-1m-2, after an analysis of the Glopnet dataset; Wright et al.2004, by Schymanski et al.2007), and MA,p is the perennial vegetation cover fraction.

The costs for root maintenance (Rr) were defined as

(8) R r = c Rr r r 2 S A , r ,

where cRr is the respiration rate per fine root volume (0.0017 mol s−1 m−3) and rr the root radius (set to 0.3×10-3 m), which 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. SA,r represents the root surface area per unit ground area (m2 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 (Rv) are therefore just a function of rooting depth and vegetated cover:

(9) R v = c rv M A y r ,

where crv is the cost factor for water transport (mol m−3 s−1), MA the fraction of vegetation cover (–), and yr the rooting depth (m). The cost factor crv was set to 1.0 µmolm-3s-1 by Schymanski et al. (2015) after a sensitivity analysis for Howard Springs, which is also adopted here.

Based on the carbon cost functions and the assimilated carbon by photosynthesis (Ag), we can calculate the net carbon profit:

(10) NCP = A g ( t ) - R f ( t ) - R r ( t ) - R v ( t ) d t ,

with t representing the time step.

2.2.4 Long-term optimization

The rooting depths of the perennial trees and the seasonal grasses (yr,p and yr,s) as well as the foliage-projected cover of the perennial vegetation (MA,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 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 Farquhar1977), represented by λp and λs (mol mol−1 m−1) for perennial and seasonal vegetation, respectively, and the sum of water suction heads (hi) 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 ir,p and ir,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 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 establishment 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 (Cowan and Farquhar1977; 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 yr and MA,p, to maximize the total NCP over the entire simulation period.

2.2.5 Short-term optimization

Some vegetation properties, such as seasonal vegetation cover (MA,s), the electron transport capacities at 25 C for the seasonal and perennial vegetation (Jmax25,s and Jmax25,p) and root surface area distributions of the seasonal and perennial vegetation (SAdr,i,s and SAdr,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 kr, defined as

(13) k r = 0.95 M q x - M q , min 0.05 M q x ,

with Mqx the maximum tissue water content (kg m−2) and Mq,min the minimum daily tissue water content (kg m−2). The minimum daily tissue water content is not allowed to be less than 0.9⋅Mqx (i.e. stomata close when it approaches this value) and cannot exceed the maximum tissue water content Mqx. The coefficient of change (kr) ranges between 1 if maximum tissue water depletion was reached during the day (if Mq,min=0.9Mqx) and −1 if tissue water was not depleted at all (Mq,min=Mqx). After a value of kr is calculated at the end of a day, the relative effectiveness of the roots in the different layers kr,eff,i is evaluated:

(14) k r , eff , i = 0.5 Q r , daily , i S A , r , i max Q r , daily , 1 S A , r , 1 , , Q r , daily , n S A , r , n ,

with Qr,daily,i the daily root water uptake in layer i (m s−1), SA,r,i the root surface area in layer i (m2 m−2) and n the number of layers. Eventually, the new root surface areas per layer are determined based on the factors kr and kr,eff and the maximum growth per day (Gr,max, set to a value of 0.1 m2 m−3; see Schymanski et al.2008b):

(15) Δ S A , r , i = G r , max k r k r , eff , i .

The other vegetation properties are 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 simulated day the combination of values that would have achieved the maximum NCP on the present day is selected for the next day. These vegetation parameters always vary on a daily basis, even though the time step of the VOM is usually hourly or sub-hourly. Only the stomatal conductances, as determined by Cowan–Farquhar optimality, are varied over an hourly time step.

2.2.6 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 NCP over the entire simulation period (37 years for VOM-v0.5, from 1 January 1980 until 31 December 2017). The SCE algorithm uses first an initial random seed and 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 h is reduced if any state variable in the model changes by more than 10 % per time step.

2.2.7 Meteorological data

A relatively long time series of meteorological inputs is required to run and optimize the VOM. The necessary meteorological data include 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 (Jeffrey et al.2001). In addition, the VOM requires information about atmospheric CO2 levels, which were provided as a constant in the VOM-AoB2015 version, whereas for VOM-v0.5 used here, we used the Mauna Loa CO2 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, with a diurnal variation that was imposed for the global radiation, temperature and vapour pressure deficit, as described in detail in Appendix A of Schymanski et al. (2009).

Observed atmospheric CO2 levels at the flux tower were not used due to the required length of the time series for the VOM (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 were replaced for the days that flux tower observations were available (max. 6 %; see Fig. S4.3 in Supplement S4). See also Supplement S3, Fig. S3.1 for the time series of meteorological data.

2.2.8 Model evaluation data

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 LE. LE was converted to evapotranspiration (ET), defined here as the sum of all evaporation and transpiration processes, even though these processes are different in nature (Savenije2004). Eventually, the gap-filled observations of GPP and ET were compared with the modelled fluxes. The VOM-AoB2015 was originally run until the end of 2005, and for this reason, the modelled fluxes were evaluated for the overlapping period between model and flux tower observations from 7 August 2001 until 31 December 2005. For consistency, the VOM-v0.5 runs were evaluated for the same time period.

To evaluate the foliage-projected cover (FPC) dynamics of seasonal and perennial vegetation predicted by the VOM, defined as the sum of MA,p and MA,s, we used satellite-derived monthly fractions of photosynthetically active radiation absorbed by vegetation (fPAR) from Donohue et al. (2008, 2013), which were converted into estimates of FPC. The maximum possible value of fPAR was defined as 0.95 by Donohue et al. (2008) and relates to maximum projective cover (i.e. FPC = 1.0). The linear relation of FPC to fPAR data (Asrar et al.1984; Lu2003) allowed for the calculation of FPC by dividing the fPAR values by the maximum value of 0.95.

Table 2Vertical profile of soil characteristics at Howard Springs, based on data from the Soil and Landscape Grid of Australia (Viscarra Rossel et al.2014a, b, c), in addition to field measurements of Jason Beringer and Lindsay B. Hutley. Here, θr refers to the residual moisture content, θs to the saturated water content, α and n to the Van Genuchten soil parameters (van Genuchten1980) and Ksat to the saturated hydraulic conductivity.

Download Print Version | Download XLSX

2.2.9 Systematic analysis of modifications to the VOM set-up

To assess the effect of the modifications to the VOM set-up from versions VOM-AoB2015 to VOM-v05, each individual change was added to the reference set-up of 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)

    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 January 1976 until 31 December 2005. This was done in order to check the reproducibility of the results of Schymanski et al. (2015).

  2. Re-run SCE

    The VOM was re-optimized with the same settings and input data as the VOM-AoB2015 and run with the same model period from 1 January 1976 until 31 December 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 CO2 levels

    Atmospheric CO2 levels were originally assumed constant in the VOM-AoB2015 with CO2 concentrations of 350 ppm, but in the VOM-v0.5, these were taken from the Mauna Loa CO2 records (Keeling et al.2005). Therefore, the VOM-AoB2015 was run with variable CO2 levels and optimized for the period 1 January 1976 until 31 December 2005. This was done in order to assess the sensitivity of the model to variable CO2 levels.

  4. Reduced soil layer thickness

    The soil layer thickness was set to 0.2 m, instead of the 0.5 m used in the VOM-AoB2015, after running a sensitivity analysis with the VOM-v0.5 (see Supplement S2). The VOM-AoB2015 was also run with 0.2 m now and optimized for the period 1 January 1976 until 31 December 2005, in order to assess the influence of different soil layer thicknesses on the VOM-AoB2015 results.

  5. Variable atmospheric pressure

    A new version of the meteorological data from the Australian SILO Data Drill (Jeffrey et al.2001) provided time series of atmospheric pressure starting from 1 January 1980, whereas originally this had been fixed at a level of 1013 hPa for the VOM-AoB2015. The variable atmospheric pressure was included in the VOM-AoB2015 and the model was optimized for the period 1 January 1980 until 31 December 2005 due to the available time series from 1 January 1980. This was performed to assess the importance of precise atmospheric pressure data for the VOM simulations.

  6. Optimized grass rooting depth

    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 also let 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 January 1976 until 31 December 2005.

  7. Updated meteorological data

    A new version of the meteorological data from the Australian SILO Data Drill (Jeffrey et al.2001) was used in the VOM-v0.5, starting from 1 January 1980. Therefore, the VOM-AoB2015 was also optimized with the new meteorological data for the period 1 January 1980 until 31 December 2005. The time series of daily maximum and minimum temperatures, shortwave radiation, precipitation, vapour pressure and atmospheric pressure were all updated but the CO2 concentrations kept fixed at 350 ppm. In this way, it can be assessed to what extent the different data versions lead to different model results.

  8. Updated and extended meteorological data

    The new version of the meteorological data from the Australian SILO Data Drill (Jeffrey et al.2001) includes more recent years, which were included in the VOM-v0.5 simulations. To find out how far the inclusion of more recent meteorological forcing alone affected the results, we also re-ran the VOM-AoB2015 with the extended meteorological forcing but a constant atmospheric CO2 concentration of 350 ppm. Therefore, the model period of the VOM-AoB2015 was extended and the model was optimized from 1 January 1980 until 31 December 2017.

  9. Modified hydrology

    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 the absence of more detailed knowledge about these parameters. At the same time, the drainage level zr and total soil thickness cz were set to 10 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.

    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 Fig. 2, required for a systematic comparison with other model applications in the accompanying paper of Nijzink et al. (2022), with a total soil thickness cz of 30 m, a fast drainage parameterization with a drainage level zr of 5 m (i.e. 25 m below the surface), a length scale for seepage outflow Λs set to 2 m and a channel slope γ0 set to 0.02 rad.

    Therefore, the VOM-AoB2015 was optimized with the new hydrological schematization of the VOM-v0.5 for the period 1 January 1976 until 31 December 2005. In this way, the effect of the hydrological settings on the model results can be assessed.

  10. Modified soil properties

    The soils were assumed vertically homogeneous in the VOM-AoB2015 but were parameterized in the VOM-v0.5 based on field measurements of sand, clay and silt content provided by Jason Beringer and Lindsay B. Hutley in the top 10 cm and the Soil and Landscape Grid of Australia (Viscarra Rossel et al.2014a, b, c) for the deeper layers. The soils were classified 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 tables5 from Carsel and Parrish (1988). See also Table 2 for the soil parameterization in the 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 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 January 1976 until 31 December 2005 in order to assess the changes due to the different soils.

  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 how far these changes compensated for each other, they were both implemented together in the VOM-AoB2015 while keeping everything else unmodified.

  12. VOM-v0.5

    Eventually, all changes were applied to the VOM-AoB2015, resulting in the new VOM-v0.5 simulations, as presented in the accompanying paper by Nijzink et al. (2022).

Table 3Vegetation properties in the Vegetation Optimality Model optimized for maximizing the NCP.

Download Print Version | Download XLSX

(Jeffrey et al.2001)(Jeffrey et al.2001)

Table 4Modifications to VOM-AoB2015.

Download Print Version | Download XLSX

Figure 3Relative changes in the mean annual values of the fluxes for the different (incremental) changes, as described in Sect. 2.2.9, in comparison with the VOM-AoB2015, for (a) ET, (b) GPP, (c) transpiration perennials (trees), (d) GPP perennials (trees), (e) transpiration seasonals (grasses), (f) GPP seasonals (grasses), (g) soil evaporation and (h) combined water use efficiency (WUE) of seasonal and perennial vegetation.


Figure 4Relative changes in the vegetation properties for the different (incremental) changes, as described in Sect. 2.2.9, in comparison with the VOM-AoB2015, for (a) the water use parameter cλe,p for perennial vegetation, (b) the water use parameter cλe,s for seasonal vegetation, (c) the water use parameter cλf,p for perennial vegetation, (d) the water use parameter cλf,s for seasonal vegetation, (e) projected cover perennials (trees), (f) mean projected cover seasonals (grasses), (g) root depth perennials (trees) and (h) root depth seasonals (grasses).


3 Results and discussion

3.1 Effects of modifications to the VOM

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 set-up 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 Fig. 3. Without optimizing the VOM, but still including the modifications, the sensitivities in GPP and ET were even smaller (see Supplement S1, Fig. S1.51), except for a strong increase in soil evaporation after changing the soils. See also Supplement S1 for detailed time series.

As expected, re-running the VOM-AoB2015 (Case 1 in Table 4) with the originally optimized parameters resulted in negligible differences. Re-running the optimization (Case 2) did result in slightly different results (12.6 % higher projective cover for the perennials), but none of the simulated fluxes changed by more than 10 % (Fig. 3). Similarly to the fluxes, the changes in vegetation parameters for reproducing VOM-AoB2015 and re-running the optimization algorithm remained small (Fig. 4).

In contrast, changing the fixed atmospheric CO2 levels (350 ppm) in the VOM-AoB2015 to variable atmospheric CO2 levels (Case 3) had a relatively large influence on perennial vegetation, yielding values of GPP for perennial vegetation that were up to 21.0 % higher (Fig. 3d). Note that the CO2 levels of the Mauna Loa records have a mean of 369 ppm and a maximum of 410 ppm during the modelling period, i.e. mostly higher than the 350 ppm prescribed in the VOM-AoB2015. See Fig. S3.1f in Supplement S3 for more details about the CO2 levels used here. Interestingly, the implementation of variable CO2 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 (Fig. 4e), while the seasonal cover was on average reduced, which relates to the generally elevated CO2 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 (Fig. 4f).

Changing the vertical soil discretization of 0.5 m in the VOM-AoB2015 to a finer resolution of 0.2 m (Case 4) had a minor influence, with a change less than the variability due to re-running the optimization algorithm (i.e. 2.6 % in the resulting GPP and 0.3 % in ET, Fig. 3a and b). The reduced soil layer thickness mainly affected the water use parameters cλe,p and cλe,s (Fig. 4a and b, respectively). At the same time, the root depths for the perennials increased (Fig. 4g), which compensated for the change in water use, resulting in only minor changes in the fluxes (Fig. 3).

The variable atmospheric pressures (Case 5) only had a minor influence as well (5.0 % change in GPP and 0.4 % change in ET), which could also relate to re-running the optimization algorithm. It led to changes in the vegetation parameters as well, but this stayed limited to a maximum of 25 % for cλe,p (Fig. 4a). However, this is stronger than the observed changes in the resulting fluxes of ET and GPP (Fig. 3), which remained much smaller, related also to the non-linear relationship between cλ, ET and GPP.

Similarly, when the grass rooting depths were optimized (Case 6) instead of the prescribed rooting depth of 1 m, simulated GPP and ET were changed by 1.0 % and −2.3 %, respectively (Fig. 3a and b). The optimization led to shallower grass roots of 0.5 m (incurring lower carbon costs) and, therefore, to reductions in GPP and ET. This was accompanied by increased cλe,s and reduced cλf,s (Fig. 4b and d), pointing to a more efficient water use strategy with less water transpired per assimilated CO2 (Fig. 3h).

The updated meteorological input data, for the runs until 31 December 2005 (Case 7) as well as the extended runs until 31 December 2017 (Case 8), hardly influenced the outcomes, with less than 10 % relative change in the resulting fluxes (Fig. 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 (24.8 %, Fig. 4e) from 31.7 % to 39.6 %. This happened as well for re-running the SCE algorithm, and the changes related to the updated meteorological input could be attributed to re-running the optimization algorithm as well. However, the water use strategy parameters (Fig. 4a–d) also changed, with an especially strong change for cλf,p with −27.8 % (Fig. 4c), but the resulting total water use efficiency again remained similar (Fig. 3h).

The implementation of free draining conditions (Case 9) had strong effects on the simulated fluxes, with lower values of both ET and GPP (−20.4 % and −6.9 %, respectively, Fig. 3a and b). However, here especially the simulated transpiration of perennial vegetation was reduced, whereas the transpiration by seasonal vegetation stayed relatively similar (Fig. 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 therefore mostly affects the perennial and not so much the seasonal vegetation. The modified hydrology led to larger changes in the vegetation properties, with a particularly strong decrease in cλf,p (Fig. 4c) and to a lesser degree cλf,s (Fig. 4d). Hence, the modified hydrology led to a more efficient water use (Fig. 3h). At the same time, the root depths for the perennial vegetation strongly reduced as well (Fig. 4g).

Even stronger effects were found for the updated soil texture, which resulted in slightly reduced ET (−1.1 %) but clearly increased GPP by 9.6 % (Fig. 3a and b). The increase in simulated GPP was largely due to increased GPP by perennial vegetation, which at the same time slightly decreased its transpiration. This coincided with a largely increased perennial vegetation cover and reduced rooting depth compared to the original simulations (vegetation cover increased from 0.31 to 0.43, while rooting depth declined from 4 to 3.5 m). Overall, the perennial vegetation benefited from the greater carry-over of soil moisture from the wet season into the dry season (Supplement S1, Fig. 1.40c) due to finer soil texture and associated increased water retention and reduced drainability. More water availability during the dry season allowed an increased perennial vegetation cover to be maintained. On the other hand, the finer soil texture resulted in higher suction heads during the wet season (Supplement S1, Fig. 1.40d), making water uptake more “expensive” as more fine roots would be needed to achieve the same water uptake rates as in coarser-textured soil. This resulted in a new optimum with enhanced perennial vegetation cover and, therefore, increased GPP (Supplement S1, Fig. 1.37f), while perennial water use was only moderately increased in the dry season and even slightly reduced during the wet season (Supplement S1, Fig. 1.37b), leading to largely increased water use efficiency per ground area (Fig. 4h).

Combining the new soils with the new hydrological settings (Case 11) still resulted in a reduction in ET by −11.5 %, whereas their combined effect on GPP led to only a small reduction by −1.1 % (Fig. 3a and b). Here, the reduction of ET occurred mainly during the dry season and was related to reductions in the perennial transpiration (Fig. 5a–c), whereas the GPP stayed relatively similar (Fig. 5e). At the same time, updating the soils and hydrology combined resulted in a more moderate increase in the perennial vegetation cover (Fig. 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 a large increase in GPP. Hence, the finer soil texture largely compensated for the effect on GPP but not on ET.

Figure 5Comparison between the results of the VOM-AoB2015 in green, simulations using new soil and hydrological parameterization (red) and simulations using all changes in combination (black). (a) ET, (b) transpiration by perennials (trees), (c) transpiration seasonals (grasses), (d) soil evaporation, (e) GPP, (f) GPP perennials (trees), (g) GPP seasonals (grasses) and (h) projective cover. Time series in (a)(g) were smoothed using a moving average of 7 d. The daily average quality flags of the flux tower observations are shown as a dashed line in panel (e), with a value of 100 for a completely gap-filled day and 1 for gap-free observations.


Figure 6Ensemble years of evapotranspiration (ET) and gross primary productivity (GPP) for the VOM-v0.5 (dark green), flux tower observations (black) and VOM-AoB2015 (light green), all smoothed by a 7 d moving average. The ensemble years are calculated for the overlapping time periods with the flux tower observations (7 August 2001 until 21 December 2016).


Figure 7Simulated and observed hydrological state variables at Howard Springs. (a) Groundwater depths, with dashed lines representing the prescribed bedrock depths, dotted lines the rooting depths (trees in red and grasses in orange), the VOM-v0.5 in red, the VOM-AoB2015 in green, and observations of three different boreholes in the vicinity of the study site in blue (Northern Territory Government, Australia2018); (b) the volumetric water content in the upper soil layer with the VOM results in red and the results of Schymanski et al. (2015) in green and measurement-based values at 5 cm depth at the flux tower sites in blue; (c) the total water storage in the root zone for the VOM-v0.5 (red, storage until 2.6 m) and the VOM-AoB2015 in green (storage until 4.0 m); (d) the suction heads for the VOM-v0.5; and (e) the suction heads of the VOM-AoB2015. Observed water tables in panel (a) were obtained from the water data portal (Northern Territory Government, Australia2018, accessed 8 March 2019).


3.2 Comparing VOM-v0.5 and VOM-AoB2015: resulting differences and underlying mechanisms

Eventually, all the changes were incorporated into the VOM-AoB2015 (Case 1), resulting in the VOM-v0.5 (Case 12), as used in the accompanying paper of Nijzink et al. (2022). Previously, we identified the isolated effects of these changes, but here we explore the combined effects of the most important modifications.

The relative error for mean annual ET changed for the VOM-AoB2015 from an overestimation by 8.4 % to an underestimation of −10.2 % for the VOM-v0.5, whereas the relative error for the mean annual GPP changed from 17.8 % to 14.7 % overestimation. The GPP especially showed differences during the transition from the wet to dry season (Fig. 6b). Both the VOM-v0.5 and VOM-AoB2015 overestimated the GPP in comparison with the flux tower observations, but the modifications reduced this overestimation for the VOM-v0.5, mainly caused by a reduction of GPP by the perennial trees (Fig. 3d). However, the changes in ET were stronger and the water use efficiency strongly increased (Fig. 3h).

The ensemble years in Fig. 6 revealed that the ET was most strongly underestimated by the VOM during the dry season at Howard Springs. The observed groundwater tables (Fig. 7a) ranged from 5 to 15 m depth seasonally, whereas the VOM was parameterized now to keep groundwater tables close to 25 m depth, required for the accompanying paper of Nijzink et al. (2022). Schymanski et al. (2015) originally assumed a much shallower drainage level at Howard Springs, which led to groundwater tables around 5 m depth and better correspondence to the observed fluxes (Fig. 6).

The simulated soil moisture in the top soil layer, as illustrated in Fig. 7b, 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. 7b. Observed soil moisture in the upper 5 cm showed similar patterns to the simulated soil moisture in the top soil layer. The modelled soil moisture was generally higher, as this represents the integrated soil moisture over 0.5 m (VOM-AoB2015) and 0.2 m (VOM-v0.5), whereas the observed soil moisture in the top 5 cm is expected to be lower due to soil evaporation and percolation. The total water storage in the root zone was higher during the dry season in the new model simulations compared to Schymanski et al. (2015) (Fig. 7c) but lower during the wet season. The water retention curves (Fig. 8) also show a clear shift, especially for the layers below 0.4 m, indicating extra storage.

However, the water suction heads again showed strong similarities between the current model runs and the results of Schymanski et al. (2015) (Fig. 7d and e, respectively) but with differences reflecting the vertical resolutions of the soil domains and simulated rooting depths and generally higher values for the VOM-v0.5 for the deeper layers. Hence, the water storage in the root zone also increased during the transition from the wet to dry periods (Fig. 7c), but resistivity increased as well due to lower hydraulic conductivities. Changing the hydraulic conductivity in isolation clearly leads to reductions in ET for the perennials (Supplement S1, Fig. S1.41b) due to the increased resistivity. Nevertheless, this is in the final VOM-v0.5 compensated for by the newly optimized water use strategy parameters and a higher water use efficiency (Fig. 3h). However, the new soil structure has the largest effect, as can be observed when running the VOM-v0.5 (i.e. with the new soils) with the water use strategy parameters of the VOM-AoB2015 in Fig. S1.56b of Supplement S1.

Figure 8Water retention curves at Howard Springs for the VOM-v0.5 results (red: top two layers, blue: deeper layers) and the results of the VOM-AoB2015 in green. Note that multiple lines are shown for the VOM-v0.5 due to the different soil parameterizations per soil layer in the current model runs, whereas the VOM-AoB2015 used one soil parameterization for all the soil layers. The upper soil layers have however the same soil parameters, leading to overlapping curves (red and dashed green).


4 Conclusions

As models, input data and parameterizations evolve, and the effects of individual improvements are usually not systematically evaluated. Here, we analyse the effects of changes to the Vegetation Optimality Model (VOM) between versions VOM-AoB2015 and VOM-v0.5, the latter of which is the basis of a companion paper (Nijzink et al.2022). Some of the modifications were done for improved realism and 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 CO2 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 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 analysing their effect in combination.

This analysis revealed that updated soil textures and a changed hydrological schematization had a strong influence on the results (see Figs. 3 and 4). An underestimation of dry season ET at Howard Springs (Fig. 6) was much more apparent when compared to the results of Schymanski et al. (2015), where the drainage parameterization maintained a water table depth much closer to the observed water table at this site. The effect of a much deeper groundwater table in the present 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 field capacities (Fig. 8) in the deeper soil layers compared to the simulations by Schymanski et al. (2015). The use of variable atmospheric CO2 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 CO2 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. This is in line with other studies (e.g. York et al.2002; Bierkens and van den Hurk2007; 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. (2022), potential bias in simulation results has to be 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 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 process (Blyth et al.2011; Abramowitz2012; Clark et al.2021). Here, we provided an example of 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 set-up 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. The availability of the original VOM-AoB2015 code was a pre-requisite for 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.

Code and data availability

Model code is available on GitHub (; Schymanski2022). Release v0.5 is used in this study (; Nijzink and Schymanski2020). The full analysis, including all scripts and data is available on renku (; Nijzink2022) and a static version of this repository can be found on Zenodo (; Nijzink and Schymanski2021).


The supplement related to this article is available online at:

Author contributions

SJS and RCN designed the set-up of the study. Model code was originally developed by SJS but updated and modified by RCN. RCN set up the repositories for the pre-processing, modelling and post-processing on renkulab. LBH and JB provided site-specific knowledge and data. The main manuscript and supplementary information were prepared by RCN, together with input from SJS. LBH and JB provided corrections, suggestions and textual inputs for the main manuscript.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This work used eddy covariance data collected by the TERN-OzFlux facility (, last access: 18 January 2022). OzFlux would like to thank the financial support of the Australian Federal Government via the National Collaborative Research Infrastructure Scheme and the Education Investment Fund.

We thank the SILO Data Drill hosted by the Queensland Department of Environment and Science for providing the meteorological data (, last access: 8 March 2019).

We thank the Scripps CO2 programme (, last access: 2 July 2019) for the Mauna Loa Observatory records.

We also thank CSIRO for the Soil and Landscape Grid of Australia (, last access: 23 October 2019) and the Australian monthly fPAR derived from Advanced Very High Resolution Radiometer reflectances – version 5 (, last access: 5 April 2019).

We thank the Northern Territory Water Data WebPortal for the groundwater data (, last access: 3 March 2021).

We would like to thank the three anonymous reviewers and the editor Tomomichi Kato for their constructive feedback and comments, which were very helpful in improving the manuscript.

Financial support

This research is part of the WAVE project funded by the Fonds National de la Recherche (FNR) Luxembourg (grant no. A16/SR/11254288).

Review statement

This paper was edited by Tomomichi Kato and reviewed by three anonymous referees.


Abramowitz, G.: Towards a public, standardized, diagnostic benchmarking system for land surface models, Geosci. Model Dev., 5, 819–827,, 2012. a, b

Allen, R. G., Pereira, L. S., Raes, D., and Smith, M.: Crop evapotranspiration – Guidelines for computing crop water requirements, Irrigation and drainage paper 56, FAO – Food and Agriculture Organization of the United Nations, Rome, 1998. a

Asrar, G., Fuchs, M., Kanemasu, E. T., and Hatfield, J. L.: Estimating Absorbed Photosynthetic Radiation and Leaf Area Index from Spectral Reflectance in Wheat1, Agronom. J., 76, 300,, 1984. a

Beringer, J., Hutley, L. B., McHugh, I., Arndt, S. K., Campbell, D., Cleugh, H. A., Cleverly, J., Resco de Dios, V., Eamus, D., Evans, B., Ewenz, C., Grace, P., Griebel, A., Haverd, V., Hinko-Najera, N., Huete, A., Isaac, P., Kanniah, K., Leuning, R., Liddell, M. J., Macfarlane, C., Meyer, W., Moore, C., Pendall, E., Phillips, A., Phillips, R. L., Prober, S. M., Restrepo-Coupe, N., Rutledge, S., Schroder, I., Silberstein, R., Southall, P., Yee, M. S., Tapper, N. J., van Gorsel, E., Vote, C., Walker, J., and Wardlaw, T.: An introduction to the Australian and New Zealand flux tower network – OzFlux, Biogeosciences, 13, 5895–5916,, 2016. a, b

Beringer, J., McHugh, I., Hutley, L. B., Isaac, P., and Kljun, N.: Technical note: Dynamic INtegrated Gap-filling and partitioning for OzFlux (DINGO), Biogeosciences, 14, 1457–1460,, 2017. a

Bierkens, M. F. P. and van den Hurk, B. J. J. M.: Groundwater convergence as a possible mechanism for multi-year persistence in rainfall, Geophys. Res. Lett., 34, L02402,, 2007. a, b

Blyth, E., Clark, D. B., Ellis, R., Huntingford, C., Los, S., Pryor, M., Best, M., and Sitch, S.: A comprehensive set of benchmark tests for a land surface model of simultaneous fluxes of water and carbon at both the global and seasonal scale, Geosci. Model Dev., 4, 255–269,, 2011. a, b

Buckley, T. N. and Roberts, D. W.: DESPOT, a process-based tree growth model that allocates carbon to maximize carbon gain, Tree Physiol., 26, 129–144,, 2006. a

Carsel, R. F. and Parrish, R. S.: Developing joint probability distributions of soil water retention characteristics, Water Resour. Res., 24, 755–769,, 1988. a, b

Clark, M. P., Zolfaghari, R., Green, K. R., Trim, S., Knoben, W. J. M., Bennett, A., Nijssen, B., Ireson, A., and Spiteri, R. J.: The Numerical Implementation of Land Models: Problem Formulation and Laugh Tests, J. Hydrometeorol., 22, 1627–1648,, 2021. a, b

Cowan, I. R. and Farquhar, G. D.: Stomatal Function in Relation to Leaf Metabolism and Environment, in: Integration of activity in the higher plant, edited by: Jennings, D. H., Cambridge University Press, Cambridge, 471–505, 1977. a, b, c

Donohue, R., McVicar, T., and Roderick, M.: Australian monthly fPAR derived from Advanced Very High Resolution Radiometer reflectances – version 5. v1. CSIRO, Data Collection,, 2013. a, b

Donohue, R. J., Roderick, M. L., and McVicar, T. R.: Deriving consistent long-term vegetation information from AVHRR reflectance data using a cover-triangle-based framework, Remote Sens. Environ., 112, 2938–2949,, 2008. a, b

Duan, Q., Sorooshian, S., and Gupta, V. K.: Optimal use of the SCE-UA global optimization method for calibrating watershed models, J. Hydrol., 158, 265–284,, 1994. a

Givnish, T.: Adaptation to Sun and Shade: a Whole-Plant Perspective, Funct. Plant Biol., 15, 63,, 1988. a

Hikosaka, K.: A Model of Dynamics of Leaves and Nitrogen in a Plant Canopy: An Integration of Canopy Photosynthesis, Leaf Life Span, and Nitrogen Use Efficiency, The American Naturalist, 162, 149–164,, 2003. a

Hutley, L. B., Beringer, J., Isaac, P. R., Hacker, J. M., and Cernusak, L. A.: A sub-continental scale living laboratory: Spatial patterns of savanna vegetation over a rainfall gradient in northern Australia, Agricult. Forest Meteorol., 151, 1417–1428,, 2011. a, b, c

Jeffrey, S. J., Carter, J. O., Moodie, K. B., and Beswick, A. R.: Using spatial interpolation to construct a comprehensive archive of Australian climate data, Environ. Model. Softw., 16, 309–330,, 2001. a, b, c, d, e, f, g, h, i

Keeling, C. D., Piper, S. C., Bacastow, R. B., Wahlen, M., Whorf, T. P., Heimann, M., and Meijer, H. A.: Atmospheric CO2 and 13CO2 Exchange with the Terrestrial Biosphere and Oceans from 1978 to 2000: Observations and Carbon Cycle Implications, in: A History of Atmospheric CO2 and its effects on Plants, Animals, and Ecosystems, edited by: Ehleringer, J. R., Cerling, T. E., and Dearing, M. D., Springer Verlag, New York, 83–113,, 2005. a, b, c

Lu, H.: Decomposition of vegetation cover into woody and herbaceous components using AVHRR NDVI time series, Remote Sens. Environ., 86, 1–18,, 2003. a

Ma, X., Huete, A., Yu, Q., Coupe, N. R., Davies, K., Broich, M., Ratana, P., Beringer, J., Hutley, L. B., Cleverly, J., Boulain, N., and Eamus, D.: Spatial patterns and temporal dynamics in savanna vegetation phenology across the North Australian Tropical Transect, Remote Sens. Environ., 139, 97–115,, 2013. a

Maxwell, R. M., Chow, F. K., and Kollet, S. J.: The groundwater–land-surface–atmosphere connection: Soil moisture effects on the atmospheric boundary layer in fully-coupled simulations, Adv. Water Resour., 30, 2447–2466,, 2007. a, b

Medlyn, B. E., Dreyer, E., Ellsworth, D., Forstreuter, M., Harley, P. C., Kirschbaum, M. U. F., Roux, X. L., Montpied, P., Strassemeyer, J., Walcroft, A., Wang, K., and Loustau, D.: Temperature response of parameters of a biochemically based model of photosynthesis. II. A review of experimental data, Plant Cell Environ., 25, 1167–1179,, 2002. a

Nijzink, R. C.: VOMcases, RenkuLab [code/data], available at:, last access: 27 January 2022. a

Nijzink, R. C. and Schymanski, S. J.: schymans/VOM: Code used for 2020 paper on the NATT (v0.5), Zenodo [code],, 2020. a, b

Nijzink, R. and Schymanski, S.: VOMcases (v0.3), Zenodo [code/data],, 2021. a, b

Nijzink, R. C., Beringer, J., Hutley, L. B., and Schymanski, S. J.: Does maximization of net carbon profit enable the prediction of vegetation behaviour in savanna sites along a precipitation gradient?, Hydrol. Earth Syst. Sci., 26, 525–550,, 2022. a, b, c, d, e, f, g, h, i, j, k

Northern Territory Government, Australia: Water data portal, available at: (last access: 3 March 2021), 2018. a, b

Radcliffe, D. E. and Rasmussen, T. C.: Soil water movement, in: Physics Companion, edited by: Warrick, A. W., CRC Press, Boca Raton, Fla, 85–126, ISBN 9781420041651, 2002. a

Raupach, M. R.: 19 Dynamics and Optimality in Coupled Terrestrial Energy, Water, Carbon and Nutrient Cycles, in: Predictions in Ungauged Basins: International Perspectives on State of the Art and Pathways Forward, edited by: Franks, S., Sivapalan, M., Takeuchi, K., and Tachikawa, Y., IAHS Press, Wallingford, 16, ISBN 978-1901502381, 2005. a

Reggiani, P., Sivapalan, M., and Hassanizadeh, S. M.: Conservation equations governing hillslope responses: Exploring the physical basis of water balance, Water Resour. Res., 36, 1845–1863,, 2000. a, b

Rodríguez-Iturbe, I., D'Odorico, P., Porporato, A., and Ridolfi, L.: On the spatial and temporal links between vegetation, climate, and soil moisture, Water Resour. Res., 35, 3709–3722,, 1999a. a

Rodríguez-Iturbe, I., D'Odorico, P., Porporato, A., and Ridolfi, L.: Tree-grass coexistence in Savannas: The role of spatial dynamics and climate fluctuations, Geophys. Res. Lett., 26, 247–250,, 1999b. a

Savenije, H. H. G.: The importance of interception and why we should delete the term evapotranspiration from our vocabulary, Hydrol. Process., 18, 1507–1511,, 2004. a

Schenk, H. J. and Jackson, R. B.: Rooting depths, lateral root spreads and below-ground/above-ground allometries of plants in water-limited ecosystems, J. Ecol., 90, 480–494, 2002. a

Schymanski, S.: VOM, GitHub [code], available at:, last access: 27 January 2022. a

Schymanski, S. J., Roderick, M. L., Sivapalan, M., Hutley, L. B., and Beringer, J.: A test of the optimality approach to modelling canopy properties and CO2 uptake by natural vegetation, Plant Cell Environ., 30, 1586–1598,, 2007. a, b, c, d, e, f, g

Schymanski, S. J., Roderick, M. L., Sivapalan, M., Hutley, L. B., and Beringer, J.: A canopy-scale test of the optimal water use hypothesis, Plant Cell Environ., 31, 97–111,, 2008a. a, b

Schymanski, S. J., Sivapalan, M., Roderick, M. L., Beringer, J., and Hutley, L. B.: An optimality-based model of the coupled soil moisture and root dynamics, Hydrol. Earth Syst. Sci., 12, 913–932,, 2008b.  a, b, c, d, e, f

Schymanski, S. J., Sivapalan, M., Roderick, M. L., Hutley, L. B., and Beringer, J.: An optimality-based model of the dynamic feedbacks between natural vegetation and the water balance, Water Resour. Res., 45, W01412,, 2009. a, b, c, d, e, f, g, h, i, j, k, l

Schymanski, S. J., Roderick, M. L., and Sivapalan, M.: Using an optimality model to understand medium and long-term responses of vegetation water use to elevated atmospheric CO2 concentrations, AoB Plants, 7, plv060,, 2015. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab

van Genuchten, M. T.: A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Sci. Soc. Am. J., 44, 892–898,, 1980. a, b, c

Viscarra Rossel, R., Chen, C., Grundy, M., Searle, R., Clifford, D., Odgers, N., Holmes, K., Griffin, T., Liddicoat, C., and Kidd, D.: Soil and Landscape Grid National Soil Attribute Maps – Clay (3′′ resolution) – Release 1, CSIRO Data Access Portal,, 2014a. a, b

Viscarra Rossel, R., Chen, C., Grundy, M., Searle, R., Clifford, D., Odgers, N., Holmes, K., Griffin, T., Liddicoat, C., and Kidd, D.: Soil and Landscape Grid National Soil Attribute Maps – Silt (3′′ resolution) – Release 1, CSIRO Data Access Portal,, 2014b. a, b

Viscarra Rossel, R., Chen, C., Grundy, M., Searle, R., Clifford, D., Odgers, N., Holmes, K., Griffin, T., Liddicoat, C., and Kidd, D.: Soil and Landscape Grid National Soil Attribute Maps – Sand (3′′ resolution) – Release 1, , CSIRO Data Access Portal,, 2014c. a, b

von Caemmerer, S.: Biochemical Models of Leaf Photosynthesis, in: Techniques in Plant Sciences, CSIRO Publishing, Collingwood, 2,, 2000. a

Whitley, R., Beringer, J., Hutley, L. B., Abramowitz, G., De Kauwe, M. G., Duursma, R., Evans, B., Haverd, V., Li, L., Ryu, Y., Smith, B., Wang, Y.-P., Williams, M., and Yu, Q.: A model inter-comparison study to examine limiting factors in modelling Australian tropical savannas, Biogeosciences, 13, 3245–3265,, 2016. a, b, c, d, e, f, g, h, i, j

Williams, R. J., Duff, G. A., Bowman, D. M. J. S., and Cook, G. D.: Variation in the composition and structure of tropical savannas as a function of rainfall and soil texture along a large-scale climatic gradient in the Northern Territory, Australia, J. Biogeogr., 23, 747–756,, 1996. a

Wright, I. J., Reich, P. B., Westoby, M., Ackerly, D. D., Baruch, Z., Bongers, F., Cavender-Bares, J., Chapin, T., Cornelissen, J. H. C., Diemer, M., Flexas, J., Garnier, E., Groom, P. K., Gulias, J., Hikosaka, K., Lamont, B. B., Lee, T., Lee, W., Lusk, C., and Midgley, J. J.: The worldwide leaf economics spectrum, Nature, 428, 821–827,, 2004. a

York, J. P., Person, M., Gutowski, W. J., and Winter, T. C.: Putting aquifers into atmospheric simulation models: an example from the Mill Creek Watershed, northeastern Kansas, Adv. Water Resour., 25, 221–238,, 2002. a, b

Short summary
The Vegetation Optimality Model (VOM) is a coupled water–vegetation model that predicts vegetation properties rather than determines them based on observations. A range of updates to previous applications of the VOM has been made for increased generality and improved comparability with conventional models. This showed that there is a large effect on the simulated water and carbon fluxes caused by the assumption of deep groundwater tables and updated soil profiles in the model.