Temporal Dynamics of Canopy Properties and Carbon and Water Fluxes in a Temperate Evergreen Angiosperm Forest

Abstract

The forest–atmosphere exchange of carbon and water is regulated by meteorological conditions as well as canopy properties such as leaf area index (LAI, m² m⁻²), photosynthetic capacity (PC μmol m⁻² s⁻¹), or surface conductance in optimal conditions (G_s,opt, mmol m⁻² s⁻¹), which can vary seasonally and inter-annually. This variability is well understood for deciduous species but is poorly characterized in evergreen forests. Here, we quantify the seasonal dynamics of a temperate evergreen eucalypt forest with estimates of LAI, litterfall, carbon and water fluxes, and meteorological conditions from measurements and model simulations. We merged MODIS Enhanced Vegetation Index (EVI) values with site-based LAI measurements to establish a 17-year sequence of monthly LAI. We ran the Community Atmosphere Biosphere Land Exchange model (CABLE-POP (version r5046)) with constant and varying LAI for our site to quantify the influence of seasonal canopy dynamics on carbon and water fluxes. We observed that the peak of LAI occurred in late summer–early autumn, with a higher and earlier peak occurring in years when summer rainfall was greater. Seasonality in litterfall and allocation of net primary productivity (F_NPP) to leaf growth (a_f, 0–1) drove this pattern, suggesting a complete renewal of the canopy before the timing of peak LAI. Litterfall peaked in spring, followed by a high a_f in summer, at the end of which LAI peaked, and PC and G_s,opt reached their maximum values in autumn, resulting from a combination of high LAI and efficient mature leaves. These canopy dynamics helped explain observations of maximum gross ecosystem production (F_GEP) in spring and autumn and net ecosystem carbon loss in summer at our site. Inter-annual variability in LAI was positively correlated with Net Ecosystem Production (F_NEP). It would be valuable to apply a similar approach to other temperate evergreen forests to identify broad patterns of seasonality in leaf growth and turnover. Because incorporating dynamic LAI was insufficient to fully capture the dynamics of F_GEP, observations of seasonal variation in photosynthetic capacity, such as from solar-induced fluorescence, should be incorporated in land surface models to improve ecosystem flux estimates in evergreen forests.

1. Introduction

Forest–atmosphere exchanges of carbon and water are driven by atmospheric conditions, soil properties, and forest canopy structural variables such as leaf area index (LAI) and height [1,2,3]. Remotely sensed LAI has improved modeling of terrestrial carbon and water fluxes in temperate and subtropical coniferous and broadleaved forests in China [4] but the relationship between LAI and primary production weakened in mesic evergreen broadleaved forests [5]. Inappropriate representation of canopy structure may contribute uncertainties to predictions and simulations of ecosystem fluxes [6,7]. Understanding the interactions between forest canopy dynamics, carbon and water fluxes, and climate will provide new insights that are required to improve land surface models (LSMs) for understanding past and projecting future changes in coupled climate–ecosystem behavior [8,9,10,11].

In many temperate forests, especially deciduous forests, LAI has a clear seasonal pattern, with a period of high LAI in the (summer) growing season and a leafless period in the (winter) dormant season. Deciduous canopy phenology has been studied extensively and shown to have been altered by recent climate change [12,13,14]. For instance, growing season length is increasing with rising air temperature (T_air) in temperate deciduous forests, playing a key role in the interannual variability of carbon fluxes [15,16]. The lengthening of the growing season is also linked to the enhanced amplitude of the seasonal cycle of atmospheric [CO₂] [17,18,19,20]. In evergreen forests, in contrast, less is known about seasonal variation in LAI and its impact on fluxes [21,22,23]. In some evergreen forests, LAI is relatively stable, with environmental conditions and leaf age driving the seasonality of carbon and water exchange [24,25,26,27]. However, in other evergreen forests, particularly those with a relatively short leaf lifespan, seasonal variation in LAI can be significant [28,29]. For example, in Amazon tropical forests, variations in LAI and leaf physiology were both important for explaining the seasonal variability of photosynthetic capacity and gross ecosystem productivity (F_GEP) [25,30,31]. Moreover, a land surface model predicted an inverse relationship between leaf area and carbon fluxes in an Amazon rainforest where greening occurs during the dry season until canopy phenology was constrained with satellite-derived LAI observations [9].

Changes in LAI over time are driven by the timing of leaf growth and litterfall. Leaf growth depends on net primary production (F_NPP) and fractional allocation of F_NPP to leaves (a_f), while litter fall primarily depends on leaf age but can also be influenced by environmental stress (such as heat and drought stress, wind, and insect attacks) [32,33]. The underlying mechanisms controlling carbon allocation in forests are not well understood [34], leading to large discrepancies in simulated LAI among LSMs [35]. Currently, LSMs use empirical relationships, such as fixed coefficients, allometric relations, or optimization theory, to determine carbon allocation to leaf growth [36]. The phenology of leaf allocation has rarely been considered, particularly in evergreen forests (but see [37]). It is widely assumed that leaf phenology is driven by seasonal variation in F_NPP alone. Similarly, few models consider the phenology of litterfall in evergreen forests, although there is evidence of environmental controls (e.g., temperature and photoperiod) on the timing of litterfall. Leaf abscission can be triggered by water stress in response to the production of abscisic acid [38]; increased leaf life span is linked to soil nutrient limitation and lower mean annual temperature [39]. A better understanding of the timing of allocation and litterfall would improve our understanding of ecosystem flux phenology.

LAI can be measured on the ground with a range of methods [28,40]. It can also be tracked with remotely sensed spectral vegetation indices (greenness) used as proxies, such as the MODIS (Moderate Resolution Imaging Spectroradiometer) normalized difference vegetation index (NDVI) or MODIS enhanced vegetation index (EVI) [41,42,43,44]. These proxies are useful to extrapolate from ground-based measurements that are limited in time and space, as these remote sensing products have been available globally since the early 2000s. Long satellite-derived time series can help gain a clearer understanding of LAI dynamics but need to be calibrated against on-ground measurements. Greenness indices have also been shown to be correlated with other ecosystem functions such as F_GEP [45,46], net ecosystem production (F_NEP) [47], and photosynthetic capacity (PC, the maximum productivity of a canopy, i.e., F_GEP in optimal meteorological conditions) [24]. The time-integrated MODIS EVI products, at a spatial resolution of 500 m, can generate proxies calibrated on flux tower data. At the global scale, EVI is a good predictor of F_GEP in cropland and deciduous (grassland) biomes but is typically less so in evergreen forests where EVI varies less seasonally [24,48]. In Australia, where high variability in precipitation leads to strong variations in soil water availability, EVI has been found to be a good predictor of F_GEP and PC in monsoonal savannas and in semi-arid mulga woodland. However, the strength of the correlation varies for eucalypt forests with subtle canopy dynamics. EVI correlation with PC and LAI can be poor in wet-sclerophyll eucalypt forests (AU-Tum; [24]) but good in temperate evergreen forests (AU-Cum; [49]).

Our primary aim in this study was to quantify the dynamics of canopy properties in an evergreen forest at the Cumberland Plain Terrestrial Ecosystem Research Network (TERN) site near Sydney, Australia. This endangered ecosystem hosts unique biodiversity and provides ecosystem services for the growing human population but is under threat from increasingly extreme climatic conditions and urbanization. We used a combination of measurements and the CABLE-POP land surface model to address three questions. (1) How did LAI respond to environmental drivers over a 17-year period (2001–2017)? (2) What was the seasonality of weather, canopy properties, and ecosystem fluxes? (3) What did we learn by comparing measurements with simulations with varying or constant LAI? Our aim is to improve understanding and modeling of the dynamics of canopy properties and thus fluxes of evergreen forests.

2. Materials and Methods

2.1. Site Description

The Cumberland Plain site (‘AU-Cum’ in FLUXNET http://sites.fluxdata.org/AU-Cum/, accessed on 30 June 2018; ‘Cumberland Plain’ http://www.ozflux.org.au/monitoringsites/cumberlandplain/index.html, accessed on 30 June 2018), is located near Sydney, Australia, (latitude 33.61518; longitude 150.72362), with mean annual air temperature (T_air) of 17.7 °C (24.3 °C mean maximum temperature and 11.1 °C mean minimum temperature for the 1993–2018 period) and mean annual precipitation of 800 mm yr⁻¹ for the 1881–2018 period [29]. The dry-sclerophyll forest is dominated by Eucalyptus moluccana Wall ex Roxb. and E. fibrosa F. Muell., with >70% of trees infected by mistletoe toward the end of the study period (Amyema miquelii Lehm. ex Miq.); the forest also has a mid-canopy stratum of Melaleuca decora (Salisb.) ex Britten, a shrub layer dominated by Bursaria spinosa Cav., and a sparse understory composed of forbs, grasses, and ferns [29]. The canopy height is ~24 m, with individual taller trees (airborne lidar survey, November 2015). A 30-m tall flux tower, with instruments at 29 m height, is installed near the geographical center of the remnant Cumberland Plain woodland.

2.2. Meteorological Data

Environmental conditions were measured on site during the 4-year period of 2014–2017. The detail of the methods can be found in [49]. In brief, meteorological data were logged every 30 min to match the time step of the tower flux data, including T_air and relative humidity (HMP155A, Vaisala, Vantaa, Finland sensors at 29 m), photosynthetic photon flux density (F_PPFD, LI190SB, LI-COR Inc., Lincoln, NE, USA), and incoming and outgoing shortwave and longwave radiation (NR01 radiometer, Hukseflux, Delft, The Netherlands). Rainfall was measured at an open area with a tipping bucket (TB4, Campbell Scientific Inc., Logan, UT, USA) ~2 km away from the studied forest site. Soil water content (SWC) was averaged between two CS616 water content reflectometers (Campbell Scientific Inc., Logan, UT, USA) installed horizontally at 5 cm below the soil surface.

Forcing data for the CABLE-POP model simulations for the period 2000–2017 comprised gridded rainfall, T_air; vapor pressure deficit; and solar irradiance surfaces. Daily means of meteorological data gridded at 0.05° spatial resolution were retrieved from the Bureau of Meteorology’s Australian Water Availability Project data set (BoM AWAP), [50,51]. Data for the study site were extracted and then down-scaled to an hourly time resolution [52].

2.3. Leaf Area Index

We estimated LAI_EVI over the period 2001–2017 using the following steps: first, we estimated LAI at a high temporal resolution based on canopy light transmittance during the year 2014–2017 [28] at three locations within the tower footprint, covering approximately 0.1-ha. However, this method underestimates LAI due to the extinction caused by non-leaf elements. Therefore, in a second step, we adjusted the magnitude of LAI using upward canopy photography, taken at six-month intervals at 121 locations across a 1-ha grid [53]. In a final step, in order to extrapolate LAI during the time period 2001–2017, we used a relationship between the previously adjusted LAI based on canopy photography and EVI [54]. These steps can be visualized in Figure A1.

2.4. Litterfall and Leaf Growth

We collected monthly litterfall (F_Lfall, g m⁻² month⁻¹) from July 2015 to December 2017 in the tower footprint approximately once per month, using three litter traps (0.14 m⁻² ground area) located within 5 m of each understory F_PPFD sensor (nine traps in total). For each month, we partitioned the litter into eucalypt leaves, mistletoe leaves, and other (mostly woody) components. We estimated the leaf mass per area (L_MA) of eucalyptus and mistletoe leaves by sampling approximately 50 fresh leaves of each species in June 2017 (177 g m⁻² for eucalypt, 248 g m⁻² for mistletoe). We used this L_MA to estimate the leaf litter production (P_leaf) in m² m⁻² month⁻¹ of eucalypt, mistletoe, and total as the sum of both. Then, we estimated leaf growth (F_Lgrowth, m² m⁻² month⁻¹) as the sum of the net change in LAI (ΔL) and P_leaf.

We estimated the canopy leaf carbon pool (C_leaf, g C m⁻²) as

$${\mathrm{C}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{f}}=\mathrm{L}\mathrm{A}\mathrm{I}\times {\mathrm{L}}_{\mathrm{M}\mathrm{A}}·0.5$$

(1)

where L_MA is leaf mass per area (g m⁻²) and 0.5 converts the dry mass of leaves to carbon mass. L_MA is approximated as the average of eucalypt and mistletoe L_MA, based on the observed equal contributions of each to total leaf litterfall.

2.5. Net Ecosystem Exchange of H₂O and CO₂: Observation

Surface–atmosphere exchanges of carbon (F_NEE) and water (F_ET) were continuously measured during the 4-year period of 2014–2017 using the eddy-covariance method, as described by [49]. In brief, F_ET was measured as the vertical turbulent exchange (change in storage of water were assumed to be negligible) and F_NEE was measured as the sum of vertical turbulent exchange (F_CT) and rate of change in storage (F_CS), as follows:

$${\mathrm{F}}_{\mathrm{N}\mathrm{E}\mathrm{E}}={\mathrm{F}}_{\mathrm{C}\mathrm{T}}+{\mathrm{F}}_{\mathrm{C}\mathrm{S}}$$

(2)

Net ecosystem production (F_NEP) was estimated as—F_NEE; positive values of F_NEP thus indicate uptake of CO₂ by the ecosystem.

Density of CO₂ or water vapor (open-path IRGA (LI-7500A, LI-COR Inc., Lincoln, NE, USA)) and vertical wind speed (CSAT 3D sonic anemometer (Campbell Scientific Inc., Logan UT, USA)) were measured at 29 m height, approx. 5 m above the canopy, at high frequency (10 Hz), and processed using EddyPro (EddyPro^® Software (Version 6.2) LI-COR, Inc., Lincoln, NE, USA) to calculate F_CT and F_ET and a quality check (qc) flag (0: best quality) of F_CT accounting for stationarity tests and turbulence development tests as presented in [55]. Furthermore, data with a friction velocity (u_*) below 0.2 m s⁻¹ were discarded, leading to the retention of 8800 points or 51% of the best quality raw data [49].

The change in storage of CO₂ was calculated using data from a profiler system (CO₂ measured at 8 heights), following the equation [56]

$${{\mathrm{F}}_{\mathrm{C}\mathrm{S}}=\left(\frac{\mathsf{\Delta }\mathrm{C}}{\mathsf{\Delta }\mathrm{t}}\right)}_{\mathrm{k}=1}\times {\mathrm{z}}_{\mathrm{k}=1}+{\sum }_{\mathrm{k}=2}^{\mathrm{n}}\left\{\left[{\left(\frac{\mathsf{\Delta }\mathrm{C}}{\mathsf{\Delta }\mathrm{t}}\right)}_{\mathrm{k}}+{\left(\frac{\mathsf{\Delta }\mathrm{C}}{\mathsf{\Delta }\mathrm{t}}\right)}_{\mathrm{k}-1}\right]\times \frac{{\mathrm{z}}_{\mathrm{k}}–{\mathrm{z}}_{\mathrm{k}-1}}{2}\right\}$$

(3)

where F_CS is the change in storage flux of CO₂, C is the CO₂ concentration (μmol m⁻³) measured at the height z_k (m), Δt (s) is the time between two measurements (1800 s or 30 min), and z_k is the height of a layer (between two inlets sampling CO₂).

Missing and discarded F_NEP were filled using a neural network in PyFluxPro v1.1 [57], allowing for monthly and annual budgets by integration. The neural network predicted F_NEP as a function of short-wave radiation, vapor pressure deficit, T_air, soil temperature, and soil water content. F_NEP was partitioned into F_GEP and F_ER using a neural network with 10 nodes and 500 iterations to predict ecosystem respiration (F_ER) as a function of soil temperature, T_air, and SWC [49,57]. F_GEP was calculated as F_NEP + F_ER.

Soil respiration (F_SR) was measured continuously at 30 min resolution and at three locations and used in CABLE parameterization (see below). CO₂ concentration inside a closed chamber was measured at 1 Hz over 2 min using an IRGA, from which F_SR was calculated (LI-8100, LI-COR Inc., Lincoln NE, USA) [58].

2.6. Net Ecosystem Exchange of H₂O and CO₂: Modeling

We used the CABLE-POP terrestrial biosphere model [59], driven with the following meteorological inputs: T_air, rainfall, short-wave downwelling radiation, long-wave down radiation, surface specific humidity, wind speed, surface air pressure, and CO₂ concentration. We first calibrated the model for the site, optimizing five parameters against flux data (F_NEP anomalies, F_ET, and F_SR quality filtered observations from 2014–2017) using the PEST optimization package (model-independent parameter estimation and uncertainty analysis, http://www.pesthomepage.org/, accessed on 30 June 2018). For the optimization, the model was forced with measured meteorological data measured at the site and LAI_EVI. The optimized parameter values are given in

Table 1. List of symbols and acronyms.

Acronym	Short Definition	Unit
EC	Eddy-covariance	-
Cv	CABLE-POP varying LAIEVI	-
Cc	CABLE-POP constant LAIEVI	-
FNEP	Net ecosystem production	Monthly or annual: (g C m−2) Half-hourly: (μmol m−2 s−1)
FET	Evapotranspiration	Monthly or annual: (g C m−2) Half-hourly: (μmol m−2 s−1)
FGEP	Gross ecosystem productivity	Monthly or annual: (g C m−2) Half-hourly: (μmol m−2 s−1)
FER	Ecosystem respiration	Monthly or annual: (g C m−2) Half-hourly: (μmol m−2 s−1)
LAIEVI	Leaf area index	(m2 m−2)
PC	Photosynthetic capacity	(μmol m−2 s−1)
FLfall	Litter fall	(g C m−2 month−1)
FLgrowth	Leaf growth	(g C m−2 month−1)
af	Allocation of FNEP to leaves	fraction
kf	Canopy leaf turnover	fraction month−1
Gs,opt	Surface conductance	(mmol m−2 s−1)
FPPFD	Photosynthetic Photon Flux Density	(μmol m−2 s−1)
Tair	Air temperature	(°C)
D	Vapour pressure deficit	(kPa)
Precip	Precipitation	(mm month−1)
SWC	Soil water content	(%)

and

Table 2. CABLE-POP parameters fitted to measured data using PEST. Net ecosystem production and evapotranspiration (F_ET) are from eddy-covariance data and soil respiration (F_SR) from auto-chambers.

Parameter	Description	Estimation and 95% CI	Unit
α (mol mol−1)	Quantum yield of electron transport [60]	0.14 ± 0.002	Mol mol−1
kn (dimensionless)	Extinction coefficient for leaf nitrogen with canopy depth	0.61 ± 0.015	-
g1 (kPa0.5)	Stomatal conductance parameter [61]	5.36 ± 0.085	kPa0.5
γ (dimensionless)	Sensitivity of stomatal conductance and root water uptake to SWC [62]	2.80 × 10−3 ± 0.600 × 10−3	-
Vcmax_scalar (dimensionless)	Scaling factor on prior estimate of maximum catalytic activity of Rubisco, as prescribed by [35]	0.69 ± 0.010	-

.

In order to test the influence of dynamic LAI as a regulator of ecosystem fluxes independently from weather, we then conducted two different runs of CABLE-POP: one with a fixed value of LAI_EVI (equal to the average of observed LAI_EVI) at all times and one with time-varying, observed LAI_EVI. We compared these two model fluxes outputs to the measured data (quality-filtered F_NEP and F_ET) over the 2014–2017 period. We used the two simulations of CABLE-POP (with varying LAI or fixed LAI) to estimate the seasonal and interannual variability of modeled gross ecosystem productivity (F_GEP). This allowed us to quantify how much variation in LAI and weather contribute, either alone or combined, to anomalies in modeled F_GEP. In what follows, flux outputs from the two simulations are subscripted “Cv” for variable LAI (e.g., F_GEP,Cv) and “Cc” for constant LAI (F_GEP, Cc).

We further extracted two measures of canopy performance from both model simulations and from flux data to evaluate the canopy phenology. We defined monthly photosynthetic capacity (PC) as the average of F_GEP under “optimal” environmental conditions (F_PPFD > 1000 umol m⁻² s⁻¹, vapor pressure deficit (D) < 1.3 kPa, SWC > median (SWC)) [49]. For the period 2014–2017, we calculated PC from the flux data and from the two model scenarios to compare observed seasonal and interannual variation in PC with simulations from both model runs. Similarly, we calculated the surface conductance under optimal conditions as follows. We calculated surface conductance (G_s) by solving the Penman–Monteith equation, then optimal surface conductance (G_s,opt) was calculated by sub-setting G_s under the following conditions: D between 0.9 and 1.3 kPa (to avoid stomatal closure), F_PPFD between 800 and 1200 umol m⁻² s⁻¹ (light saturated condition for photosynthesis), precipitation in the last 24 h below 0.5 mm, and precipitation in the last 12 h below 0.2 mm (to minimize evaporation contribution to F_ET).

Finally, we used modeled values of F_NPP in conjunction with leaf growth data to estimate the foliage allocation fraction. In CABLE-POP, F_NPP is estimated as the difference in F_GEP and autotrophic respiration [59]. To estimate the allocation of predicted F_NPP to leaf growth, we used a simple difference model of leaf biomass (C_leaf) [63] (Equation (4)), which assumes that the variation in C_leaf over a period of time is equal to leaf growth minus litter fall (F_Lfall) (Equation (5)). We assume that leaf growth (F_Lgrowth) is equal to the allocation to leaves multiplied by F_NPP and that F_Lfall is equal to the fraction of the initial C_leaf senescing over this period of time (Equation (6)).

$${{\mathrm{C}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{f}}}_{\mathrm{t}+1}={{\mathrm{C}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{f}}}_{\mathrm{t}}+\frac{\mathrm{d}{\mathrm{C}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{f}}}{\mathrm{d}\mathrm{t}}$$

(4)

$$\frac{\mathrm{d}{\mathrm{C}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{f}}}{\mathrm{d}\mathrm{t}}={\mathrm{F}}_{\mathrm{L}\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}-{\mathrm{F}}_{\mathrm{L}\mathrm{f}\mathrm{a}\mathrm{l}\mathrm{l}}$$

(5)

$$\frac{\mathrm{d}{\mathrm{C}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{f}}}{\mathrm{d}\mathrm{t}}={\mathrm{a}}_f{\mathrm{F}}_{\mathrm{N}\mathrm{P}\mathrm{P}}-{\mathrm{k}}_f{\mathrm{C}}_{\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{f}}$$

(6)

where F_Lgrowth is leaf growth (gC m⁻² month⁻¹), F_Lfall is litter fall (gC m⁻² month⁻¹), C_leaf is leaf biomass (gC m⁻²), t is time (month), a_f is the allocation of F_NPP to leaf growth (fraction, 0 to 1), and k_f is the turnover rate of leaf biomass due to senescence (fraction month⁻¹).

The only unknowns in these equations are the allocation of F_NPP to leaf growth (a_f = F_Lgrowth/F_NPP) and the turnover rate due to senescence (k_f = C_leaf/F_Lfall), which were calculated directly using the measured estimates of F_Lgrowth and F_Lfall, and the modeled estimate of F_NPP from the CABLE-POP simulation with varying LAI. Notably, this approach to estimating a_f assumes that the carbon used for leaf growth is supplied by that month’s F_NPP, rather than stored non-structural carbohydrates (NSC). We take this approach in order to illustrate seasonal patterns simply but acknowledge that a stored NSC pool could also play a role in the seasonal dynamics of leaf growth [64].

3. Results

3.1. Response of LAI to Environmental Drivers

There was large interannual and seasonal variability in LAI (Figure 1a). The monthly mean of LAI_EVI ranged from 0.8 to 2.0 m² m⁻² over the 4-year period of record with measurements (2014–2017), while the mean annual LAI_EVI varied from 1.03 m² m⁻² in 2003 to 1.64 m² m⁻² in 2016 (

Table 3. Annual precipitation, air temperature, leaf area index, and net ecosystem production over the period 2001–2017. Inter-annual variability of meteorology (precipitation and T_air) and leaf area index estimated from enhanced vegetation index over the 17-year period of 2000–2017. Notably, the large inter-annual range (max–min) of the leaf area index relative to its average magnitude (0.6 compared to 1.3). Bold values indicate maxima and minima.

Year	Precipitation (mm)	Tair (°C)	LAIEVI (m2 m−2)	FNEP,Cv (g C m−2)	FNEP,Cc (g C m−2)	FNEP,data (g C m−2)
2001	752.5	17.4	1.58	114.3	114.4	NA
2002	626.3	17.5	1.09	−66.0	97.0	NA
2003	651.1	17.2	1.03	−42.0	84.5	NA
2004	648.8	17.6	1.21	8.3	−1.6	NA
2005	702.4	17.7	1.36	48.6	42.6	NA
2006	479.2	17.6	1.23	15.8	−18.8	NA
2007	1022.4	17.7	1.32	23.7	38.0	NA
2008	810.6	16.8	1.44	90.1	112.1	NA
2009	691.5	18.0	1.31	−53.7	18.1	NA
2010	911.4	17.5	1.26	−11.8	29.5	NA
2011	782.9	17.2	1.38	14.6	30.2	NA
2012	880.8	17.0	1.19	33.3	120.8	NA
2013	789.2	18.0	1.39	54.6	50.1	NA
2014	714.2	18.0	1.31	19.1	58.9	101.4
2015	996.9	17.7	1.55	51.9	0.8	171.2
2016	779.5	18.3	1.64	125.1	61.9	392.7
2017	731.0	18.2	1.49	62.0	58.9	148.9
Average	763.0	17.6	1.34	28.7	52.8	203.5
Max–min	543.2	1.5	0.62	191.1	139.6	291.3

). Monthly mean LAI-EVI peaked each year between February and April (Figure 1 and Figure 2).

Precipitation is highly variable at the site. Monthly precipitation ranged from 0 mm (June 2001) to 325 mm (March 2017) (Figure 1), with summer months generally wetter than winter months (an average of 100 mm month⁻¹ in December through February and an average of 40 mm month⁻¹ in June through August, Figure 1). There was also significant interannual variability in rainfall, with annual totals ranging from 479 mm in 2006 to 1022 mm in 2007. Tair is also strongly seasonal but less variable from year to year. Monthly mean T_air ranged from 9.1 °C in July 2002 to 26.9 °C in January 2017 (Figure 1). Interannual variability in temperature is relatively small, with mean temperature ranging from 16.8 °C in 2008 to 18.3 °C in 2016.

Precipitation appeared to be the major driver of variability in LAI, although there was also a role for temperature. Monthly LAI_EVI increased most rapidly during wet summer months; the rate of increase correlated significantly with precipitation (r² = 0.20, p < 0.005, n = 204) and T_air (r² = 0.24, p < 0.005, n = 204) but was most strongly predicted by the combination of T_air and precipitation (r² = 0.29, p < 0.005, n = 204) (Figure A2). The magnitude of the yearly LAI_EVI,max was negatively correlated with the day of the year on which the maximum occurred (i.e., years with a low LAI_EVI,max had a late LAI_EVI,max, r² = 0.3, p = 0.04, n = 14, Figure 2). If summer precipitation was low, the day of the year of maximum LAI_EVI occurred later; r² = 0.35, p = 0.025, n = 14, Figure 2. Overall, high summer precipitation led to a high and early LAI_EVI,max whereas low summer precipitation led to a low and late LAI_EVI,max. On the annual time-scale, LAI_EVI had a high coherence with precipitation on an annual period, with no lag (Figure A3), and the yearly max of LAI_EVI was correlated with precipitation over the previous 13 months (r² = 0.47, n = 14, Figure A4).

3.2. Seasonality of Weather, Canopy Properties, and Ecosystem Fluxes

In Figure 3, we examine the seasonal patterns of drivers and fluxes at the site. The timing of maximum and minimum F_PPFD, T_air, and D followed the expected seasonal pattern, i.e., maximum in summer and minimum in winter (Figure 3k–m). The monthly average precipitation was maximum in summer and autumn, while SWC was maximum in autumn and winter (Figure 3n,o).

As noted above, LAI peaked in late summer (Figure 3f). The annual F_Lfall was 157 g C m⁻² or 1.57 m² m⁻² of leaves, which is close to the mean LAI for the site, implying a leaf life span of about a year. F_Lgrowth and F_Lfall both showed peaks in spring and autumn and reduced to near-zero values during winter (Figure 3h) but the peak in F_Lgrowth was larger. The peak in LAI is thus caused by a stronger seasonality in leaf growth than in leaf litterfall.

Photosynthetic capacity, PC, estimated using eddy-covariance data, generally followed by LAI_EVI with a time lag, peaking in autumn, i.e., one or two months after LAI_EVImax, but also showed a second smaller peak in spring (Figure 3f,g). Eddy-covariance optimal surface conductance, G_s,opt, had a seasonal pattern similar to those of LAI_EVI and SWC (Figure 3f,j,o).

High LAI and high PC in autumn led to a maximum of F_GEP,EC in autumn, with a lower peak in spring, and minimum values in summer and winter (Figure 3c). F_ER,EC peaked in summer and had a larger seasonal amplitude than F_GEP,EC. As a result, F_NEP,EC was positive in winter and negative in summer (Figure 3a). Interestingly, F_ET,EC was in anti-phase with F_NEP,EC, as F_ET,EC was maximum in summer and minimum in winter (Figure 1a,b and Figure 3). This pattern results from seasonal F_NEP following the seasonality of F_ER, while F_ET follows the seasonality of F_GEP owing to the coupling between photosynthetic uptake and transpiration.

3.3. Effects of Varying versus Constant LAI

After parameterization using PEST (see Table 2 for parameter descriptions and values and

Table A1. Parameter correlation matrix (Pearson’s r). Refer to Table 2 for the description and units.

	α	kn	g1	γ	Vcmax
α	1	−0.12	−0.54	4.453 × 10−3	0.29
kn	−0.12	1	0.34	0.11	0.58
g1	−0.54	0.34	1	0.26	−0.32
γ	4.45 × 10−3	0.11	0.26	1	0.12
Vcmax_scalar	0.29	0.58	−0.32	0.12	1

for parameter correlation matrix), we simulated the fluxes using CABLE-POP, first with varying and then with constant mean LAI. Modeled half-hourly values of F_NEP and F_ET are compared against non-gap-filled observations for the period 2014–2017 in Figure A5. Both sets of simulations were well correlated with observations and showed little bias. The varying LAI simulations showed a marginally higher r² (for F_NEP, 0.595 vs. 0.581, and for F_ET, 0.62 vs. 0.61).

Seasonal cycles of CABLE-POP F_ET and F_NEP with varying LAI were similar to those in eddy-covariance observations (anti-phase, with maximum F_ET in summer and maximum F_NEP in winter (Figure 3a,b)) but CABLE-POP did not entirely reproduce the observed variability in the data. Extremes, particularly in winter, were not reproduced by the model (Figure 1b,c). The simulated seasonal cycles of the component fluxes F_ER and F_GEP were generally similar to those observed in eddy-covariance data but the eddy-covariance derived values of both fluxes (F_GEP,EC and F_ER,EC) were generally higher in magnitude than the modeled values, particularly in autumn (Figure 3c,d). For example, in March, F_GEP,EC was 147.7 g C m⁻² month⁻¹ on average while F_GEP,Cv was only 103.8 g C m⁻² month⁻¹. The discrepancies between modeled and observed F_GEP and F_ER tend to cancel each other out, except in winter, when the model overestimates F_ER and underestimates F_GEP, leading to an underestimation of the winter peak F_NEP. Modeled F_GEP also produced a dual peak, in spring and autumn, but monthly F_GEP,EC was maximum in March, shortly after LAI_EVI peaked, while F_GEP,Cv was maximum in September, coinciding with one of the lowest values of monthly LAI_EVI (Figure 3c,f).

CABLE-POP correctly simulated the sensitivity of light-saturated F_NEP to D and SWC: the slope of F_NEP versus D was similar between the model and the observations (Figure A6). However, the modeled F_NEP was lower at low D. The magnitude of PC_Cv was lower than the magnitude of PC_EC and the variability of monthly PC_Cv (1.5 μmol m⁻² s⁻¹) was much lower than that of PC_EC (4.9 μmol m⁻² s⁻¹) (Figure 3g). In contrast, the variability and seasonal patterns of measured and modeled optimal surface conductance, G_s,opt, were similar, albeit lagged in time (Figure 3j).

Modeled F_NPP,Cv is given by the difference between F_GEP,Cv and autotrophic respiration. F_NPP,Cv had two peaks during the year, a first peak in spring and a second lower peak in autumn (Figure 3e) as autotrophic respiration was high during the summer owing to the high temperatures. The timing of these peaks of F_NPP,Cv differs strongly from the peak of F_Lgrowth, which was maximum in summer (Figure 3h). As a result, the inferred foliage allocation fraction a_f was not constant throughout the year but peaked strongly in summer, ranging from 0.009 in August to 0.57 in January (Figure 3i). In contrast, the inferred turnover rate, k_f, was less strongly seasonal (ranging from 0.02 in July to 0.15 in October).

Given the strong seasonal pattern of LAI_EVI, we anticipated that driving CABLE-POP with constant LAI would degrade the performance of the model. We found relatively small differences between the simulations with constant and variable LAI; nonetheless, we found improved agreement between simulations and observations when the variable LAI was applied rather than the constant LAI (Figure A5). For F_NEP, the variable LAI simulation had a lower seasonal amplitude compared to the constant LAI simulation (27.8 vs. 33.1 g C m⁻² month⁻¹, Figure 3a) but it also had a markedly larger inter-annual variability (Figure 4, Table 3). F_NEP,Cc ranged from −18.8 to 120.8 g C m⁻² yr⁻¹ while F_NEP,Cv ranged from −66 to 125.1 g C m⁻² yr⁻¹. On an annual timescale, none of F_NEP,Cv, F_NEP,Cc or F_NEP,EC correlated with climate variables (annual precipitation or average T_air) but both F_NEP,Cv and F_NEP,EC were correlated with mean annual LAI_EVI (Figure 4, F_NEP,Cv vs. LAI_EVI: r² = 0.72, n = 17 F_NEP,EC vs. LAI_EVI: r² = 0.68, n = 4). These results suggest that interannual variations in F_NEP were more strongly regulated by interannual LAI than by the climate variables investigated.

Over the four-year period 2014–2017, PC_EC was significantly correlated with LAI_EVI on a monthly time scale (slope = 8.0 ± 1.4, r² = 0.40, p < 0.001, n = 48) (Figure 5). PC_Cv was also significantly correlated with LAI_EVI but with a lower slope (slope = 3.5 ± 0.4, r² = 0.65, p < 0.001, n = 48). In contrast, PC_,Cc was not correlated with LAI_EVI, emphasizing the importance of variable LAI as a major contributor to the observed variability in PC.

Across the entire 17-year period of simulation, the simulated mean annual F_NEP,Cv was lower than the simulated F_NEP,Cc (Table 3). In years where mean LAI was similar to the long-term average LAI used in the Cc simulations, F_NEP,Cv was lower than F_NEP,Cc. This difference results because optimal conditions for photosynthesis occur in September–October, when LAI is lowest in the Cv simulation.

Notably, F_NEP,Cc and F_NEP,Cv fluxes appear similar in Figure 1 but, when integrated over an annual budget (Table 3), they are actually quite different, often sink vs. source of similar magnitude (e.g., 2002, 2003, 2006, and 2010 are −66.0 vs. 97.0, −42.0 vs. 84.5, 15.8 vs. −18.8, and −11.8 vs. 29.5 g C m⁻² for F_NEP,Cv and F_NEP,Cc, respectively). This is probably because F_NEP is large and positive in winter and large and negative in summer but its annual sum is close to 0.

4. Discussion

4.1. LAI Responses to Environmental Drivers

LAI and a_f in evergreen broadleaved forests are often assumed to be constant throughout the year in LSMs because observed changes can be subtle. New leaves typically emerge in spring and summer in Australian eucalypt forests but new growth can be balanced by concurrent leaf losses [65,66,67,68,69,70,71]. In our analysis, we also found that leaf growth occurred in the summer months (December through February) and we further observed that climate drivers could give more information on the precise timing and magnitude of leaf growth in summer. LAI_EVI yearly maxima were driven by summer precipitation, with a high and early peak when summer precipitation was high (Figure 2). Furthermore, legacy effects of precipitation during the past year were also important in driving the yearly maximum of LAI_EVI (Figure A4). Extreme climate events may also affect changes in LAI [72], such as the 2018–2020 drought associated with the Black Summer bushfires, and La Nina-induced flooding events in subsequent years; however, such anomalous conditions were not observed during this study.

Our estimate of F_NPP allocation to leaf growth was very high in summer (0.57) and low in winter (0.009) (Figure 3i). This allocation pattern differs strongly from the default pattern in CABLE-POP, which is a constant value of 0.3. Understanding the dynamics of carbon allocation is a key challenge to improving vegetation modeling [34,36,73,74,75,76]. Identifying the mechanisms regulating this observed pattern will be important to improve the representation of evergreen forests in land surface models. It is important to note, however, that trees may be using carbon reserves for leaf growth in summer, rather than changing the allocation of new photosynthates to leaf growth, as our estimate of a_f relied on the assumption that only recent F_NPP could be used for leaf growth. Moreover, ecosystem-scale carbon uptake from flux tower observations was not clearly associated with carbon allocation to stem and canopy growth in a different dry sclerophyll eucalypt forest [67], consistent with observations from forests more broadly [77]. Detailed estimates of carbon allocation to all pools (leaves, roots, wood, and storage) are required including where carbon is being sourced to grow new tissues (recently acquired carbon or old stored carbon) [64,67].

We also found that litterfall rates varied seasonally. We observed that leaf turnover was highest in spring and continued through summer but was negligible in autumn and winter. This means that nutrient availability from decomposing litter was mostly available in summer.

A further question of interest is whether the dynamics of LAI are adapted to environmental constraints [78,79]. According to our simulations, a higher F_NEP would be achieved if LAI were to be constant rather than varying over time (Table 3). The simulations also demonstrate that the optimal conditions for carbon uptake occur in spring, which coincides with the lowest value of LAI. This mismatch between the timing of F_NEP and LAI suggests that the timing of leaf flush is not optimized to maximize carbon uptake but may rather be constrained by other factors, such as nutrient or moisture availability. For example, the soils at the site have very low phosphorus availability [80] and leaf production may be timed to coincide with litterfall in order to maximize the trees’ ability to capture nutrients released from decomposing litter. Alternatively, [81] demonstrates that the seasonal pattern of LAI and photosynthetic capacity may enable the trees to avoid loss of conductivity during dry periods. Further exploration of the costs and benefits of the phenology of leaf area dynamics in evergreen forests is warranted.

4.2. Seasonality of Weather, Canopy Properties, and Ecosystem Fluxes

Previous studies in other sclerophyll forests in Southeast Australia found the ecosystems to be a net carbon sink all year long, particularly in summer [67,82,83,84]. In contrast, the Cumberland Plain site shows a different seasonal pattern, acting as a net source of CO₂ in summer and as a net sink in winter [49]. This unusual seasonality of F_NEP was simulated well by CABLE-POP, whereas the seasonality of F_GEP was not (Figure 3), suggesting that some canopy properties were not adequately represented in either version of the model. We discuss this in further detail below.

4.3. Effects of Varying versus Constant LAI

Application of the CABLE-POP model allowed us to test whether an LSM can reproduce the unusual pattern of F_NEP at our site and what the key mechanisms causing it are. In particular, we tested whether variable LAI was necessary to capture this pattern.

After parameter optimization, CABLE-POP with varying LAI successfully reproduced the seasonal F_NEP pattern of the Cumberland Plain site (Figure 3a). Half-hourly F_NEP and F_ET flux responses to D and SWC were similar between the observation and model (Figure A6), albeit with a small underestimation by the model at low D (<1.3 kPa). This model result showed the capacity of the optimal stomatal conductance model to capture these responses [61]. This seasonal pattern was also shown by CABLE-POP with constant LAI, indicating that the seasonal variation in LAI was not a major driver of the switch from sink to source in summertime.

There were some important discrepancies between the data and the model during specific seasons, particularly in autumn, with F_GEP being underestimated by the model, irrespective of whether LAI was constant or varying (Figure 3c). This seemed to be caused by a variation in PC not captured in the model (Figure 3g and Figure 5). A possible reasonable explanation would be the absence of age-dependent leaf-level photosynthetic capacity, which was observed to be highest in late summer when leaves were mature, at a site within 2 km of the flux tower [85]. The importance of the age-dependent physiology of leaves to determine F_GEP seasonality also contributed to seasonal variations in carbon fluxes in a tropical evergreen forest [25]. Understanding and incorporating this age dependence of leaf physiology could potentially improve LSMs [86,87,88,89].

Although variable LAI was not important for seasonal patterns, it was key to understanding the inter-annual variability of F_NEP in the model (Figure 4). Annual F_NEP, modeled with both varying and constant LAI, was not correlated with precipitation or T_air but only F_NEP modeled with varying LAI was correlated with LAI_EVI. The four years of observation had similar results, albeit with lower significance due to the low number of annual eddy-covariance observations. The stronger slope of annual F_NEP,EC vs. LAI_EVI was consistent with the higher slope of PC_EC vs. LAI_EVI compared to the model (Figure 5), which may result from age-dependent leaf physiology parameters not incorporated in the model [25,90,91]. These results suggest that model performance at the inter-annual time scale cannot be described adequately by comparing half-hourly fluxes as we did in Figure A5. However, with only four years of data, we could not assess the inter-annual performance of the model.

5. Conclusions

We found that LAI varied seasonally and this dynamic (timing and magnitude) was correlated with summer precipitation. We observed that litterfall was the highest in spring, followed by high leaf growth during summer. The inferred allocation of F_NPP to leaf growth varied seasonally (from 0 in winter to 0.6 in summer). As a result, LAI peaked at the end of summer, with canopy properties (PC and G_s,opt) peaking in early autumn.

Canopy properties in evergreen forests are often oversimplified in LSMs. We used a combination of measurements (meteorological drivers, CO₂ and water fluxes, litterfall, and LAI) and modeling (CABLE-POP, with variable or constant LAI) from the Cumberland Plain ecosystem to study the inter-annual and seasonal dynamics of canopy properties (LAI, PC, and G_s,opt) of an evergreen forest and demonstrated that the temporal variability of canopy properties of this forest type cannot be assumed to be negligible in LSMs.

The CABLE-POP model [59] with varying LAI reproduced interannual variability but was not sufficient to reproduce seasonal variation in fluxes in the Cumberland Plain site. Based on our findings, we suggest several improvements to be made to LSMs: (1) variable C allocation patterns to leaves to capture seasonal and interannual variation in the timing of leaf growth; (2) seasonal variation in leaf turnover rates; and (3) variation in photosynthetic capacity with leaf age [25]. It would be valuable to determine if similar variations are observed in other temperate evergreen forests to develop general parameterizations of these processes for LSMs. Recent developments in satellite remote sensing of solar-induced fluorescence (SIF) could also be promising to track the phenology of evergreen broadleaf forest canopy properties in the future when long-time series become available [92].