Observing the Whole-Canopy Short-Term Dynamic Response to Natural Step Changes in Incident Light: Characteristics of Tropical and Temperate Forests

Abstract

We examine the physical and biological responses of forest canopies to step changes in light caused by passing low cumulus clouds that intermittently block the direct solar beam. Using data obtained at a tropical rainforest and at a midlatitude deciduous forest, we estimate the course of sensible heat flux, net ecosystem exchange, evapotranspiration, and water-use efficiency in response to the rapid changes in the incident radiative flux. To describe these fluxes during the interval over which the effects of stomatal time delays can be most influential, eddy fluxes are estimated over minute or shorter intervals by invoking a conditional-sampling procedure based on forming a Reynolds-average ensemble. The most important differences between the two forests’ physical responses are in the thermal balances and heat-flux time response constants. During the initial period after the light transition the only mean variables that show appreciable changes are the blackbody and air temperatures, the other scalars being little affected. We find that a distinct transient thermal internal boundary layer appears ≈ 20 m thick above the temperate deciduous forest and ≈ 45 m thick above the tropical rainforest. At each forest, the effective thickness of the inferred thermal outer-canopy ‘big leaf’ is about 1 mm. Twenty minutes after the abrupt change in incident light, ensemble eddy-flux estimates approach those found using conventional time averaging, confirming the validity of the ensemble approach. Previously unrecognized transient maxima in net ecosystem exchange and evapotranspiration are evident 5–10 min following the shadow-to-light transition, longer than the average light interval between shadows observed on partly-cloudy days in each case. Short-term variations in sensible heat flux, net ecosystem exchange, and evapotranspiration approximate an exponential adjustment, implying that first-order time-dependent single-leaf models are adequate to describe whole-canopy processes in these cases, providing an experimental method for determining whole-canopy bulk stomatal time constants. During the sunlit interval (direct and diffuse radiative fluxes combined) net ecosystem exchange is enhanced, while under cloud shadow (only diffuse radiative flux) water-use efficiency increases. This light and shadow alternation provides a mechanism describing the observed enhanced net ecosystem exchange and water-use efficiency under certain types of partly-cloudy sky. We apply empirical flux-response curves to an idealized case of radiative flux varying in a regular on–off light and shadow pattern. For this case, an analytical solution for mean net ecosystem exchange flux as a function of diffuse-radiation fraction yields results that strongly resemble previous findings based on conventional time-averaged fluxes. Our analysis and modelling indicates that a well-known correlation between diffuse radiative flux and enhanced net ecosystem exchange and water-use efficiency on cloudy days is in many cases not causal, but rather to be a consequence of time averaging over light-and-shadow intervals. By linking processes associated with photosynthesis in fluctuating light at the leaf scale to the canopy scale, our efforts facilitate the scaling-up of leaf responses to the ecosystem scale.

1 Introduction

A vegetation canopy is an active, living surface—energy and scalar fluxes respond to a constantly varying lower atmospheric environment. Photosynthesis, and the resultant carbon uptake, has a complex response to changes in the photosynthetically-active radiative flux (PAR), temperature (T), vapour pressure deficit (VPD), and water potential (Ψ) (Wiegand and Swanson 1972; Jarvis 1976; Powell and Thorpe 1977; Farquhar et al. 1980; Kirschbaum and Farquhar 1984; Jarvis and McNaughton 1986; Collatz et al. 1991; Jones 1992; Pearcy et al. 1997; Gates et al. 1998; von Caemmerer 2000; Tan et al. 2017). Photosynthesis operates at a variety of time scales (e.g., Way and Pearcy 2012). Recent interest in both the ecological and the agricultural communities in what is now referred to as dynamic photosynthesis focuses on how plants and canopies overall respond to the naturally varying environment (e.g., Murchie et al. 2018; Matthews et al. 2018; Matsubara 2018). Because plant response depends on the characteristic time scale of environmental forcing relative to relevant plant time scales, Pearcy (1987) and Givnish (1988) argued that the adaptation to the different lighting environment (in the whole canopy or understorey) often plays a major role in the photosynthetic response and capacity of plants, because other factors such as temperature and humidity often change relatively slowly.

One common variable-light environment occurs as passing clouds obscure the sun, offering light and shadow intervals from 1 to 20 min, introducing sharp light-to-shadow and shadow-to-light changes in incident radiative flux. In a matter of 10–30 s, the incident radiative flux changes from being largely diffuse to being largely direct of a much larger intensity (e.g., Tomson 2014, Fig. 2). Light and shadow durations on partly-cloudy days are on the order of a few min (Kivalov and Fitzjarrald 2018), intervals comparable to the stomatal opening and closing time constants for many common tree species (e.g., Woods and Turner 1971; Vico et al. 2011). At the leaf scale, these stomatal time lags prevent plants from responding instantaneously to an abrupt change in lighting (e.g., Barradas et al. 1994; Barradas and Jones 1996; Vico et al. 2011). Forest carbon uptake is observed to increase along with water-use efficiency during partly-cloudy conditions, when compared to clear-day conditions, a fact often attributed to the more complete illumination of the canopy occasioned by the enhanced diffuse radiative flux (Freedman et al. 2001; Gu et al. 2002; Alton 2008; Urban et al. 2012). Because the incident radiative flux is reduced relative to a clear sky, we refer to this observation as the ‘cloudy-day anomaly’.

The light and shadow alternation associated with certain types of partly-cloudy skies provides a mechanism of describing the observed enhanced net ecosystem exchange and water-use efficiency. We aim to show that the cloudy-day anomaly is plausibly a consequence of averaging this type of time series. On such days, there are particular illumination and shadow intervals comparable to the stomatal response times that lead to optimal conditions for photosynthesis to increase. An intermittent-light dependency mechanism outlined below can explain the cloudy-day anomaly independent of the diffuse-radiative-flux enhancement, whose correlation between diffuse radiation and carbon uptake survives as an averaging artefact.

To address these challenges, we assess two methodological hypotheses:

• (a) Valid eddy-covariance fluxes can be obtained using the classical ensemble-averaging approach, for the case of strong ‘impulsive’ radiative-flux forcing of a forest canopy. Over the period when the canopy responds to abrupt changes in light, other environmental changes have a secondary effect. For example, changes in soil moisture do not respond rapidly enough to play a role in the 1–15-min recovery time;

• (b) The course of these short-term fluxes can be approximated using the same exponential adjustment that describes stomatal behaviour at the leaf level. This implies that valid e-folding time constants represent eddy fluxes during the critical first few minutes following transition.

We argue that hypothesis (a) can be verified if the fluxes approach the conventional time-averaged results at sufficient time following the light transition and if statistics of the turbulent fluctuations are comparable to what is known about the equilibrium roughness sublayer just above the forest canopy. Note that verifying methodological hypothesis (b) is important for demonstrating that, for the transient responses common in nature, the big-leaf model is adequate for describing whole-canopy exchange processes.

To summarize, we hypothesize that: a big-leaf model is adequate for describing whole-canopy fluxes in a variable-light environment; during the transient period, the flux time series exhibits exponential responses with characteristic time constants; as a consequence of time averaging of a series of shadow and light intervals leads to a correlation between diffuse radiative flux and carbon uptake.

1.1 Is a Big-Leaf Model Adequate in Fluctuating-Light Conditions?

Leaf to canopy scaling. Scaling-up leaf properties to canopy properties is a modelling approach, done by simply asserting the validity of some type of big-leaf model, treating the entire canopy as an equivalent single ‘big leaf’. The big-leaf extrapolations are based largely on the analogy with leaf-scale observations that the flux responses for different kinds of plants are similar (see Reich et al. 1997; Meinzer 2003). Models carrying elements of leaf-scale plant physiology have been brought nearly intact to the larger scale and used to describe whole-canopy–atmosphere exchange (Jarvis and McNaughton 1986; Ball et al. 1987; Raupach and Finnigan 1988; Jones 1992; Amthor et al. 1994; De Pury and Farquhar 1997; Dai et al. 2004; Medlyn et al. 2017; Lee 2018). For its part, evaporation is most often described using the Penman–Monteith approach (e.g. Raupach 1995). De Pury and Farquhar (1997) noted “…these big-leaf models are not truly scaled models…, in the strict sense that the variables have different definitions and relations with one another at each scale…” Reichstein et al. (2014) repeated this caution seventeen years later: “…one should not always expect a simple relation between ecosystem functions and traits at the organismic level…”

The idea of using a single big-leaf approach seems more plausible if we consider the forest-canopy geometry. Many observations indicate that certain trees (e.g., oak and maple) in mature midlatitude deciduous-forest stands form an effective monolayer of nearly closed canopy (see Kikuzawa 1995). Such species may also have similar stomatal-conductance responses (Woods and Turner 1971) and the forests’ similar leaf area indices (LAI; Horn 1971), the ratio of the total leaf area in a volume to the ground area below (Watson 1947). The leaf area index is often estimated by reference to a light extinction coefficient (e.g., Sakai et al. 1997; Bréda 2003). Closed, thick forest canopies limit ventilation, effectively cutting-off the inner-canopy layers from gas exchange with the atmosphere above (Cionco 1983; Pinker 1983). The top 5-m thick layer of the forest canopy (the outer-canopy surface) absorbs up to 90% of the incoming photosynthetically-active radiative flux (Yoda 1974; Ellsworth and Reich 1993; Parker 1995; Lalic et al. 2013), greatly reducing the light intensity in the inner canopy. Hence, the main contribution of the measured forest fluxes can be directly associated with fluxes originating in the outer-canopy volume and not so much with inner-canopy layers. High rugosity of the outer-canopy surface, a measure of the variability of the canopy-top height (Parker 1995), increases the forest surface exposed to turbulent mixing and to the direct radiative flux up to a factor of four (Parker and Russ 2004). This promotes outer-canopy coupling with the atmospheric flow, reduces occurrence of light saturation and heat stress, and limits the potential for photoinhibition. In summary, since the outer canopy assimilates the bulk of the carbon dioxide, it likely is a dominant flux source or sink region during the day.

1.2 The Eddy-Covariance Approach in a Variable Environment: Consequences of Natural Averaging

The empirical coefficients needed to calibrate big-leaf models are determined using data from eddy-covariance flux measurements (e.g., Miner et al. 2017). These eddy fluxes, taken to represent the whole canopy (Baldocchi et al. 2017), are used in these models as if the system were always nearly at equilibrium, even in a variable environment. Katul et al. (2009) and Vico et al. (2011) noted that to reach such an equilibrium may require up to a 10-min stomatal accommodation time—a situation unlikely to occur during common partly-cloudy conditions.

The eddy-flux method is widely used to obtain typical ecosystem responses for different conditions, but the calculation is most appropriate for a locally steady-state environment. Data consisting of 30-min time-averaged fluxes characteristic of clear-day, cloudy, or partly-cloudy conditions are collected and fit with empirical curves (e.g., Michaelis–Menten process, Thornley 1976) to describe the PAR—carbon-uptake dependency (Amthor et al. 1994; Gu et al. 1999; Goulden et al. 2004; Rocha et al. 2004; Hutyra 2007, Hutyra et al. 2007; Alton 2008; Urban et al. 2012). This is how the possible effect of the increased diffuse fraction in the incident radiative flux with an increase in the forest carbon uptake under partly-cloudy conditions was identified (e.g., Gu et al. 2002; Alton 2008; Urban et al. 2012). On partly-cloudy days, this standard time-averaged eddy-flux approach mixes flux contributions from both shadow and sunny periods that occur during the imposed 30-min averaging period. Are such methods adequate to investigate the consequences of a rapidly-changing incident radiative flux?

To understand whole-canopy sensitivity to a rapidly-changing natural-light environment, accurate flux observations during the first 20 min following the light transition are needed. However, statistical limitations on flux accuracy—notably the convergence-of-averages criterion (e.g., Wyngaard 2010, p. 36)—usually require the eddy flux to be evaluated at 30-min to 1-h intervals. There are scattered reports of eddy fluxes obtained at 5-min intervals. For example, Doughty and Goulden (2008) presented eddy-covariance CO₂ and water vapour fluxes at 5-min intervals at a tropical rainforest following cloud-induced abrupt light changes. They acknowledged that a considerable fraction of the flux (≈ 37%) was missed in this processing due to the loss of the low-frequency contribution (Ogura 1957; Lenschow et al. 1993; Sakai et al. 2001) but did not consider high-frequency losses. Van Kesteren et al. (2013) reported eddy fluxes at 1-min intervals using laser scintillometry, but their method relies on having a uniform surface and requires that traditional flux–gradient relations are applicable at this time scale.

McAusland et al. (2015, 2016) noted that the leaf-scale coupling between assimilation (A) and stomatal conductance (g_s) can be applicable in steady-state conditions, but assimilation and stomatal conductance are uncoupled in intermittent-light regimes, making it difficult to estimate their relation in the conventional manner (e.g., Farquhar and Sharkey 1982),

$$ A = g_{s} \left( {c_{a} - c_{i} } \right), $$

(1)

where c_a and c_i are the ambient and intercellular CO₂ molar fractions. This means that the assimilation dynamics in variable-light environment cannot be simply modelled through the stomatal-conductance variation, the Vico et al. (2011) approach,

$$ g_{s} \left( t \right) = g_{eq} - \left( {g_{eq} - g_{0} } \right){ \exp }\left\{ { - {\raise0.7ex\hbox{$t$} \!\mathord{\left/ {\vphantom {t {\tau_{g} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${\tau_{g} }$}}} \right\}, $$

(2a)

$$ \tau_{g} = \left\{ {\begin{array}{*{20}c} {\tau_{open} g_{s} \left( t \right) < g_{eq} } \\ {\tau_{close} g_{s} \left( t \right) > g_{eq} } \\ \end{array} } \right. $$

(2b)

Barradas and Jones (1996, Figs. 2 and 3), examined abrupt changes in light in laboratory observations of single leaves. We can infer that the time constant for assimilation is much smaller than that for the stomatal conductance for water vapour, which could be an illustration of this decoupling. We seek empirical evidence to describe how the net ecosystem exchange and the bulk big-leaf conductance at the canopy level are related.

Because there might be decoupling between the stomatal conductance and assimilation in fluctuating-light regimes, we must describe the whole-canopy exchange with the atmosphere more directly. The need for a whole-canopy description is even clearer once one considers the enormous extrapolation required to link properties from the single-leaf to whole-canopy scales: compare the 10³ mm² of typical leaf area with the 0.5–1-mm thickness (see Jones 1992; Vile et al. 2005) with an average seasonal eddy-flux footprint 0.23 km² for Harvard Forest and 1 km² for Tapajós National Forest (see Kim 2015; Hutyra 2007) with the top–down cumulative LAI ≈ 2–3 (see the outer-canopy discussion above).

Here is the conundrum: to understand whole-canopy responses, the eddy-flux approach is ideal, because it naturally averages over a wide area. However, to connect whole-canopy results to the leaf-scale response triggered at the time of the light transition, these eddy fluxes must be found at intervals of minutes or less.

1.3 Objectives and Research Questions

Two issues have limited progress on understanding how the natural variable-light environment affects surface water and carbon exchanges at the canopy scale:

• a. When measuring eddy fluxes, in addition to the footprint sizes, the real-world averaging process depends on the mix of the individual tree species and how deep into the canopy the turbulent mixing extends. Owing to the process uncertainty when scaling-up from the leaf to ecosystem scale, we must show observationally which fluxes (net ecosystem exchange and evapotranspiration) represent the whole canopy and which time constants are appropriate to describing their adjustment to abrupt light changes. To do this, we must identify the time history of the fluxes immediately following the light transition.

• b. Because of statistical constraints, the conventional Reynolds time-averaging approach is unsuitable to obtaining accurate whole-canopy fluxes on the zero to 10-min time scale, precisely the interval needed to assess the dynamic-photosynthetic effects caused by fluctuating light.

To address these issues and to describe the dynamics of eddy fluxes during critical periods when abrupt light changes dominate in forest ecosystems, we developed a practical observational Reynolds ensemble-average method to study non-stationary regimes, similar to the technique commonly used in laboratory experiments (e.g., He and Jackson 2000). Such ensemble averaging avoids the limitations on estimating eddy fluxes at short time intervals (Lumley and Panofsky 1964; Monin and Yaglom 1971). Even though a large amount of data is required to generate a sufficient number of realizations at each time interval to form the ensemble average, eddy-flux towers have been operating for some years. Many sites have archived extensive datasets. With time, increasingly larger datasets will be available, making our proposed method ever more viable. We present results using two such datasets obtained by our group.

To apply the ensemble average, ensemble members must be carefully chosen identifying adequately the similarity conditions needed to isolate a particular factor. The ensemble method is particularly suited to describing well-defined events, such as the shifts between shadow and sunlight studied here. We seek the response of the forest–atmosphere fluxes to a step change in incident shortwave radiative flux, starting with the transition from shadow to full sunlight or vice versa. We assert that, during the rapid response to the abrupt impulse forcing, the change from shadow to light, other relative changes in the physical or plant physiological environment are small relative to the response initiated by the light impulse. It need not necessarily be the case that every detail of the environment be exactly the same for each ensemble realization. If the response to the impulse forcing is sufficiently strong, its signal will overwhelm the effects of small fluctuations characteristic of variations in other variables.

We address three research questions to assess our hypotheses:

• 1. Can reliable dynamic outer-canopy fluxes for typical fluctuating-light conditions be obtained with the proposed eddy-covariance ensemble-flux technique?

• 2. Are the time responses of whole-canopy fluxes compatible with a big-leaf formulation?

• 3. How much do the time-varying canopy conductance and outer-canopy fluxes differ from conventional steady-state values?

First, we show that our ensemble-average fluctuations after environmental offsets are removed obey the same statistical properties as do the conventional time-averaged results (Sect. 3.1) and that our ensemble-flux method yields reliable dynamic outer-canopy fluxes for typical fluctuating-light conditions, addressing research question 1. An additional validation (Sect. 3.3) is provided by the direct comparisons between the ensemble-averaged fluxes for equilibrium conditions and the time-averaged fluxes for steady-state conditions long after the transitions.

Applying the procedure presented in Sect. 2.5, we obtain ensemble fluxes on 10- and 20-s scales (Sect. 3.3). Fitting exponential solutions, we present parametrizations of dynamic fluxes using characteristic time constants and final equilibria, the steady-state values achieved long after the transition. This addresses research question 2. Knowing the dynamic ensemble averages (Sect. 3.2) and fluxes (Sect. 3.3) allows us to estimate temporal changes of whole-canopy conductances, the thickness of any transient thermal internal boundary layer (IBL) that develops above the forests, and estimate the thermally-active big-leaf layer thickness (Sect. 3.4.2). These findings address the big-leaf compatibility issue raised in research question 2 by connecting the thermally-active big-leaf layer thickness with the outer-canopy surface.

Research question 3 is addressed in Sect. 3.3, where we show the differences of the ensemble-averaged and time-averaged fluxes obtained during the transient phase, and in Sect. 4.1, in which the ensemble-flux equilibrium values are compared with the steady-state values of the conventional time-averaged fluxes above similar forests. In Sect. 4.2 we apply our empirical results to examine an ideal case, how canopy fluxes respond to periodic incident radiative forcing. Analytic results illustrate how the empirical-flux responses so obtained successfully reproduce the enhancement in carbon uptake on cloudy days, a finding that supports our fluctuating-light-enhancement hypothesis.

2 Methods

2.1 Sites and Data Description

We analyzed two multiyear datasets. The first is taken at the Harvard Forest (HF) Environmental Measurement Site tower located in an old-growth temperate deciduous forest near Petersham, Massachusetts, USA (42.53°N, 72.17°W). The second is from the Large-Scale Biosphere–Atmosphere (LBA) Experiment km-67 site tower in the Tapajós National Forest, a tropical rainforest located near Santarém, Pará, Brazil (2.86°S, 54.96°W). Kivalov and Fitzjarrald (2018) discussed the radiation instruments deployed at each site, while Sakai et al. (1997) and Hutyra (2007), Hutyra et al. (2007) presented more details of the sensors used for the scalar observations at the HF and LBA sites respectively (see Table 1).

Table 1 Forest-site instrumentation and characteristics: longwave (L↑, L↓), global shortwave (S↑, S↓) and photosynthetically-active radiative flux (PAR↑, PAR↓) upwelling and downwelling; 3D wind velocity U = (u, v, w), air temperature T, CO₂ molar fraction c_a, water-vapour mixing ratio q

At the HF site, only data from the foliated growing season (from early May until late September) are considered (Moore et al. 1996; Sakai et al. 1997; Fitzjarrald et al. 2001). At the LBA site, about 75% of data were recorded during the dry season, (July through January). We examined seasonal differences at the LBA site with the available data and found little difference, likely because we concentrate only on partly-cloudy days. At both sites, the event selection period is 0800–1600 local standard time (LST) when convective conditions are expected.

The two towers have slightly different instrument configurations above the canopy, which leads to different upwelling longwave radiative-flux (L↑) footprints (used here to estimate the blackbody temperature at canopy top) and eddy-flux footprints (see Table 1). The HF site features flux measurements only at the 29-m level (velocity U = (u, v, w) and temperature T fluctuations) associated with CO₂ molar fraction (c_a) and water-vapour mixing ratio (q) fluctuations at the 27.5-m level, about 7 m above the canopy. At the LBA site there were two flux-measurement levels, at 57.7-m and 46.7-m heights; the upper level located about 14 m above the canopy was used for analysis. This choice of levels allows us to obtain measurements at similar relative positions within the respective roughness sublayers, about one to five times the height of the typical roughness element (see Finnigan 1985; Katul et al. 1999), here estimated by the standard deviation of the canopy height, to be about 5 m for the HF site and 10 m for the LBA site (Parker 1995; G. G. Parker, pers. comm., 2015).

We take the eddy-flux footprints (Hsieh et al. 2000; Chen et al. 2009) to be the representative areas of the whole-canopy fluxes, noting that the turbulent mixing effectively averages response characteristics of different plant species. Radiative-flux footprints were estimated following procedures similar to those outlined in Marcolla and Cescatti (2018). We recognize that the turbulent-flux footprints can vary seasonally and over smaller time periods (Baldocchi 2003; Xu et al. 2017), but for our two forest cases, the canopies are sufficiently horizontally uniform that this likely will not significantly alter our conclusions. Moreover, the larger LBA footprints may partly compensate for the greater canopy inhomogeneity related to the larger variability of the outer-canopy depth at that site.

We took instantaneous grab samples at 1 Hz (from 10-Hz eddy-system data, HF site) or 0.2 Hz (from 5-Hz eddy-system data, LBA site) to match the acquisition frequency of radiative-flux data. Note that the 1-Hz grab sampling occurs approximately at flux-integral-scale intervals (Kaimal and Finnigan 1994, Fig. 7.1; Finkelstein and Sims 2001), making for independent estimates at both HF and LBA sites.

The raw data were calibrated with reference to the hourly-averaged U, T, c_a, and q data provided by the Harvard University team (Munger and Wofsy 1999; Saleska et al. 2013). The c_a and q signals acquired from the Licor-6262 close-path gas analyzer were synchronized with the vertical-velocity output from the 3D sonic anemometer by applying a delay correction, found on a daily basis by seeking the time-averaged covariance maxima in $\overline{{w^{\prime}c_{a}^{\prime} }}$ and $\overline{{w^{\prime}q^{\prime}}}$ at each site. For the HF site, the delay was 6–7 s and the fluxes were corrected by shifting the c_a and q time series forward relative to the vertical-velocity signal. This signal synchronization also helped to minimize possible measurement discrepancies, which may have occurred due to the 1.5-m height difference in the CO₂/H₂O sensor and sonic anemometer locations. For the LBA site, the delay was about 1–2 s (Hutyra et al. 2007), and because the delay was smaller than the time resolution of the data (5-s interval sampling) no correction was applied. Our tests indicated that this approximation leads to, at most, a 2% flux underestimate.

2.2 Ensembles of Shadow-to-Light and Light-to-Shadow Conditional Events

The ensemble-event selection from the HF and LBA time series follows methodology detailed in Kivalov and Fitzjarrald (2018):

• 1. After Freedman et al. (2001), we identified the forced-cumulus-cloud conditions by calculating the lifting condensation level and comparing it with operational meteorological aerodrome reports (HF site) or ceilometer estimates (LBA site) of cloud-base heights no more than 300 m above the lifting condensation level.

• 2. We applied a 0.8 threshold to the downwelling global shortwave radiative flux (S↓) (as normalized by the estimated clear-sky radiative flux) during the identified forced-cumulus-cloud conditions. Continuous periods above the threshold were associated with the single clear periods (light events), and the continuous periods below the threshold were associated with a single cloud passage (shadow events).

• 3. Adjacent shadow and light events were combined into the single shadow-to-light and light-to-shadow conditional events. These conditional events (hereafter, simply events) were aligned in ensembles by the strong radiative-flux jump accompanying the light transition occasioned by the cumulus-cloud passages (time t = 0, Fig. 1). To satisfy radiative-flux similarity conditions (see below), each ensemble consists of the events binned by median clear-period S↓ values: [500, 700], [600, 800], [700, 900], [800, 1000], [900, 1000], [900, 1100], [1000>], [1100>] W m⁻².

  Fig. 1

  

  Incident shortwave radiative-flux ensembles of the short light events that include an initial 300 s of shadow period (a, b) and long (20 + min) light events following the shadow periods (c, d) grey—individual-event ensemble; red—ensemble mean; blue—equilibrium-curve fit into light; green—number of the events available for the time specified; [900,1100] W m⁻² clear-period S↓ bin. e, f the short-event (left graph) and long-event (right table) counts for different radiative fluxes (200 W m⁻² radiative-flux–bin size); g, h the hour of day distribution of the short events by clear-period S↓ levels (see legend) over the investigation periods; for the HF (left column) and the LBA (right column) sites

To examine the sensitivity to light change, we constructed two sets of shadow-to-light and light-to-shadow ensembles: short-duration events that begin with an adjustment period of at least 300 s (Fig. 1a, b for the shadow-to-light transition) and long-duration events, which include all possible events > 20 min in duration when the adjustment-period duration is not restricted to ≥ 300 s (Figs. 1c, d for the shadow-to-light transition). The long adjustment period before transition in the ensembles of short-duration events ensures that the plants equilibrate to the initial conditions. Though the long-duration events lack a sufficient initial time to ensure equilibrium conditions before the radiative-flux transition, they are adequate to finding the eventual equilibrium (the approximately steady-state flux) to which the system returns, the initial condition having lost importance.

Some observational evidence supports this choice of the short-event adjustment period. From Knapp and Smith (1988, Fig. 5), we note that the 5-min shadow period is sufficient to completely recover intercellular CO₂ concentration in the leaf tissues in some alpine species, recharging them before the next light period. Also, since the time constant for water-potential changes in leaf tissues is about 1.2 min (Powell and Thorpe 1977) and that for RuBP (rubilose bisphosphate, see Gross et al. 1991) regeneration is about 1.5 min (Way and Pearcy 2012), a 5-min shadow period should be sufficient for the photosynthetic systems to recover (see Wyber et al. 2017) and to replenish the water lost, so the leaves are ready to operate again when high radiative fluxes return.

Duration distributions for the short-duration events are quite non-uniform, with maxima at very short durations (about 100 s) and a rapid drop in the number of events at longer durations (details in Kivalov and Fitzjarrald 2018). For example, for the clear periods with S↓ ∈ [900,1100] W m⁻² the number of events falls from 300 to 50 within the first 500 s following transition (green lines, Fig. 1a, b and graphs Fig. 1e, f). The initial fall in the number of the 300-s shadow events (green lines, Fig. 1a, b) is due to the high variability of the data series that occurred during the first 100 s of the short cloudy intervals, which was filtered by the outlier-removal procedure. The number of long-duration events remains about the same for the entire 20-min period following transition (green lines, Fig. 1c, d and tables Fig. 1e, f).

At each site, selected short events are distributed more or less uniformly during the event selection period of 0800–1600 LST (Fig. 1g, h), with the low-S↓ over the clear-period events shifted toward afternoon after 1400 LST. We note that the most common convective-cloud onset time, just after 1000 LST, occurs as the top of the convective boundary layer approaches the lifting condensation level (Freedman et al. 2001; Berg and Kassianov 2007; Czikowsky and Fitzjarrald 2009).

Figure 1a, b indicate a large radiative-forcing jump of ≈ 600 W m⁻², starting from the shadow-period S↓ of 300–400 W m⁻² approaching values of 1000 W m⁻² in sunlit conditions, at a rate of about 30–50 W m⁻² s⁻¹ for the average 10–20-s jump. This represents a restriction determining the data rate adequate in identifying correctly the radiative-flux jumps for proper alignment in the ensembles. Kivalov and Fitzjarrald (2018) showed that both 1-s and 5-s grab samples are sufficient to correctly identify cumulus-cloud events used here.

2.3 Assessing the Degree of Similarity Among Members of the Event Ensembles

An important issue for the ensemble approach as applied to real-world situations is to ensure that for a given environmental variable the events within an ensemble are sufficiently similar to justify the calculation. Following Pearcy (1987) and Givnish (1988), we argue that the main driver of the vegetation response at this time scale is the abrupt change in sunlight, while the most notable consequent short-term environmental change is in temperature. We hypothesize that the local environmental responses to the abrupt light change dominate over the many other possible influences that distinguish events. Our radiative-flux based event-similarity condition ensures similar heating rates and consequently leads to similar photosynthetic responses. Owing to the strong relationship between the variables PAR and S↓, little difference in event definition results from the radiative-flux variable chosen.

The choice of the study period 0800–1600 LST and the degree of similarity among distinct events can be assessed in part by examining the stability of the atmospheric surface layer, denoted by (Wyngaard 2010),

$$ \frac{z - d}{L} = \frac{{\left( {z - d} \right)kg\overline{{w^{\prime}\theta^{\prime}}} }}{{ - u_{*}^{3} \theta_{v} }}, $$

(3)

where L is the Obukhov length (Monin and Obukhov 1954), d is the displacement height, z is the height, k = 0.4 is the von Karman constant, g = 9.8 m s⁻² is the acceleration due to gravity, $\overline{{w^{\prime}\theta^{\prime}}}$ is the kinematic sensible heat flux, $u_{*} = \sqrt { - \overline{{u^{\prime}w^{\prime}}} }$ is the friction velocity, and $\theta_{v}$ is a reference virtual potential temperature (≈ 300 K). During daytime convective conditions (0800–1600 LST, indicated by the vertical lines in Fig. 2a, b), the observed change in $\left( {z - d} \right)/L$ (blue lines) is associated with only a 10–20% change in the mean dimensionless gradients and variances in the roughness sublayer (see Shapkalijevski et al. 2016). This provides evidence that the atmospheric stability conditions that control turbulent mixing during all investigation periods are sufficiently similar to allow us to combine events of different times of the day and seasons into ensembles.

Fig. 2

a, b Stability index ($\frac{z - d}{L}$, blue lines) and photosynthetically-active radiative flux (PAR, red lines) hourly averages; c, d blackbody leaf (T_leaf, red) and air (T_air, blue) temperatures at 5 min after the shadow-to-light transition; e, f vapour pressure deficit (VPD) at transition (blue), 5 min (green), and 10 min (red) after the shadow-to-light transition; 25% and 75% quantile bars (dashed lines); for the HF (left column), the LBA (right column) sites

Figure 2c, d show that for both sites blackbody leaf ($T_{leaf} \approx \left( {L \uparrow /\sigma } \right)^{1/4}$) and air (T_air) temperatures, measured 5 min after the shadow-to-light transition, change little across the radiative-flux bin ranges. The same is true for the vapour pressure deficit ($V\!P\!D \approx e_{{v_{sat} }} \left( {T_{leaf} } \right) - e_{{v_{sat} }} \left( {T_{air} } \right)RH/100$) (Fig. 2e, f), where $e_{{v_{sat} }}$ is the saturation vapour pressure and RH is the relative humidity of air. Estimating vapour pressure deficit using data from two levels rather than at a single canopy level amounts to only about a 10% enhancement, an unimportant difference for the purposes of assessing the potential for water stress in the canopy. Both temperatures and vapour pressure deficits are located mostly in the ‘comfort ranges’ in which plants thrive (e.g., El-Sharkawy et al. 1986; Rylski and Spigelman 1986; Gates et al. 1998; Tan et al. 2017). Thus, we expect similar responses from the plant physiological perspective as well.

These observations are consistent with typical midday changes in conditions (Fig. 1g, h) for partly cloudy, non-precipitating days. No strong heating occurs in the middle of these cloudy days owing to the moderating effect of the radiative-flux intermittency. The [500,700] S↓-bin outlier for the HF site (Fig. 2c, e) is due to a small number of events with high relative humidity present very late in the day and should be disregarded.

Rather than seek similarity in mean scalars, we examine similarity in the time series after their environmental offsets are removed (see below). To summarize, our event-selection procedure guarantees that the selected events are uncorrelated:

• Firstly, the selected events are uniformly distributed over the observation periods (Fig. 1g, h), and they are relatively infrequent (see the event distributions in Kivalov and Fitzjarrald 2018). Hence, they can be treated as the independent ensemble realizations, representing the entire observation periods.

• Secondly, the selected adjustment periods are long enough (≥ 300 s) to allow for recovery of the leaf physiological parameters to their CO₂- and water-saturated steady-state conditions in shadow and approach the light-saturated steady-state conditions during light periods. Hence, we expect that at similar times of day characterized by similar radiative fluxes, the resulting physiological responses should also be similar. A possible clue to the validity of this approach is that the midday CO₂ flux on clear days, the so-called ‘midday rest’ (e.g., Kostychev et al. 1926, 1930; Kostychev and Kardo-Sisoyeva 1930; Kostychev and Berg 1930; Rocha et al. 2004), is about 20% smaller than the CO₂ flux on those partly-cloudy days when forced cumulus clouds prevail (Freedman et al. 2001, Fig 1.3, for the HF site), being closer to full-capacity photosynthesis response throughout the midday period.

These findings allow us to describe turbulent flows as random fields with existing probability distributions with the expected values given by the continuous function μ(t) and the variances given by the, not necessarily, continuous finite function $\sigma^{2} \left( t \right)$ (Kolmogorov 1936; Monin and Yaglom 1971). We define ensembles of independent realizations (events, time series) $\left\{ {u_{i} \left( t \right), i = 1..N} \right\}$ and their statistical properties: ensemble mean ${\langle}{u(t)}{\rangle}$ (Eq. 4), fluctuations $\left\{ {u_{i}^{\prime} \left( t \right), i = 1..N} \right\}$ (Eq. 5), variance $\sigma^{2} \left( t \right)$ (Eq. 6), and convergence (Eq. 7) in the same manner as in conventional fluid dynamics (see Kolmogorov 1936; Monin and Yaglom 1971, Wyngaard 2010),

$$ \left\langle {u\left( t \right)} \right\rangle = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} u_{i} \left( t \right), $$

(4)

$$ u^{\prime } (t) = u(t) - \left\langle {u(t)} \right\rangle , $$

(5)

$$ \sigma^{2} \left( t \right) = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} u^{\prime 2} \left( t \right), $$

(6)

$$ \left\langle {u\left( t \right)} \right\rangle \mathop \to \limits^{d} {\mathbf{N}}\left( {\mu \left( t \right),{\raise0.7ex\hbox{${\sigma^{2} \left( t \right)}$} \!\mathord{\left/ {\vphantom {{\sigma^{2} \left( t \right)} N}}\right.\kern-0pt} \!\lower0.7ex\hbox{$N$}}} \right). $$

(7)

From the central limit theorem, the ensemble mean is seen as a realization drawn from the accompanying normal distribution (Eq. 7). For $N \to \infty ,\sigma^{2} \left( t \right)/N \to 0$, hence, almost everywhere $\left\langle {u\left( t \right)} \right\rangle \to \mu \left( t \right)$, so that the ensemble mean (Eq. 4) can be represented by a continuous function: $\left\langle {u\left( t \right)} \right\rangle = \mu \left( t \right)$. The classical formulation of the central limit theorem requires the identical distribution for the realizations to use the accompanying normal distribution, defining exact similarity in realizations with the same mean and standard deviation as in Eq. 7.

We acknowledge that in field measurements, even when strict constraints are given, different realizations $\left\{ {X_{i} \left( t \right), i = 1..N} \right\}$ inevitably come with somewhat different environmental mean conditions $\left\{ {\mu_{i} = \frac{1}{{t_{2} - t_{1} }}\mathop \int \nolimits_{{t_{1} }}^{{t_{2} }} X_{i} \left( t \right)dt, i = 1..N} \right\}$, characteristic of interval $\left[ {t_{1} ,t_{2} } \right]$, which may offset them from the ensemble-mean condition Eq. 4: ($\mu_{e} = \frac{1}{N}\mathop \sum \nolimits_{i = 1}^{N} X_{i} = \frac{1}{N}\mathop \sum \nolimits_{i = 1}^{N} \mu_{i}$). To focus on the sensitivity, we remove these bulk realization offsets ≡ ($\mu_{i} - \mu_{e}$), the mean differences, and leave only the transient tendencies to be accounted for $\left\{ {X_{i}^{{\left[ {t_{1} ,t_{2} } \right]}} \left( t \right) = X_{i} \left( t \right) - \left( {\mu_{i} - \mu_{e} } \right), i = 1..N} \right\}$, where $\left[ {t_{1} ,t_{2} } \right]$ is the offset-removal time interval. This new set of the realizations has the same ensemble average as does the original set (see Appendix 3 for derivation).

Are there important consequences of the offset removal? The transient tendencies or responses for individual realizations after their offsets are removed (Fig. 3, grey lines) and their ensemble averages (green lines) show that the only strong scalar responses to the radiative-flux jump appear in T_leaf and T_air changes (Fig. 3a–d). Mean-quantity ensembles indicate that c_a and q (Fig. 3e–h) and the mean horizontal wind speed $U_{u,v} = \sqrt {u^{2} + v^{2} }$ (see Online Resource 1) vary little across the light transition, indicating that our ensembles properly isolate thermal and light effects. We expect that by combining different events into the ensembles and applying the radiative-flux–similarity constraints, we obtain valid transient responses characteristic of the light change.

Fig. 3

Ensembles of the leaf (canopy surface) blackbody (radiative) temperature (T_leaf) (a, b), air temperature (T_air) (c, d), c_a (e, f), and q (g, h) after offsets of the individual realizations were removed (grey sets of lines) and their ensemble averages weighted by the number of realizations (green lines) with the fitted smooth lines (blue); dotted black lines—their ensemble standard deviations; left axis—fluctuations around ensembles after offset removal; right axis—absolute values of respective scalars after adding the ensemble averages; solid red lines—radiative-flux jumps (see Fig. 1); [900,1100] W m⁻² clear-period S↓ bin; for the HF (left column) and the LBA (right column) sites

These features mean that we can use simple linear fits to the c_a and q ensemble averages (blue lines, Fig. 3e–h) because the mean-ensemble (green lines) changes are very small compared to the individual-realization (grey lines) fluctuations. We use the simple linear fits to the ensemble averages for the velocity as well (see Online Resource 1). For the cases of T_leaf and T_air, large ensemble excursions after the transition (t = 0, green lines, Fig. 3a–d) being comparable to the individual fluctuations (grey lines), we use the following approach-to-the-equilibrium curves (blue curves, Fig. 3a–d; Eqs. 8b and 9b) to fit the temperature ensemble averages.

Since the temperature, but not the wind speed or specific humidity, is the environmental scalar that varies most immediately following transition, a transient thermal IBL appears, into which the dynamic heat flux converges. Knowing fluxes and mean state changes, we can assess the dimensions and other properties of this layer.

The first-order thermal budget equation and its solution for T_air (derivation in Appendix 4) is

$$ \tau_{air} \frac{{dT_{air} }}{dt} = \left( {T_{e} - T_{air} } \right), $$

(8a)

$$ T_{air} \left( t \right) = T_{e} - \left( {T_{e} - T_{0}^{air} } \right){ \exp }\left\{ { - {\raise0.7ex\hbox{$t$} \!\mathord{\left/ {\vphantom {t {\tau_{air} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${\tau_{air} }$}}} \right\}, $$

(8b)

where $\tau_{air} = r_{HR} \Delta z/2$ is the thermal time constant for the air and Δz is the thickness of a transient thermal IBL that forms above the canopy. The system of equations is closed once we know the convective and radiative heat-transfer resistance of air r_HR or its conductance: $g_{HR} = 1/r_{HR}$.

The blackbody thermal budget equation and its solution for T_leaf (derivation in Appendix 4) is

$$ \tau_{leaf} \frac{{dT_{leaf} }}{dt} = \left( {\frac{{r_{HR} }}{{\rho c_{p} }}R_{neti} + T_{air} \left( t \right)} \right) - T_{leaf} , $$

(9a)

$$ \begin{aligned} T_{leaf} \left( t \right) & = \frac{{r_{HR} }}{{\rho c_{p} }}R_{neti} + T_{eq}^{air} - \left( {T_{eq}^{air} - T_{0}^{air} } \right)\frac{{\tau_{air} }}{{\tau_{air} - \tau_{leaf} }}\exp \left\{ { - {\raise0.7ex\hbox{$t$} \!\mathord{\left/ {\vphantom {t {\tau_{air} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${\tau_{air} }$}}} \right\} \\ & \quad - \left( {\frac{{r_{HR} }}{{\rho c_{p} }}R_{neti} + T_{eq}^{air} - \left( {T_{eq}^{air} - T_{0}^{air} } \right)\frac{{\tau_{air} }}{{\tau_{air} - \tau_{leaf} }} - T_{l0} } \right)\exp \left\{ { - {\raise0.7ex\hbox{$t$} \!\mathord{\left/ {\vphantom {t {\tau_{leaf} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${\tau_{leaf} }$}}} \right\}, \\ \end{aligned} $$

(9b)

where $\tau_{leaf} = \frac{{\rho^{*} c_{p}^{*} l^{*} }}{{\rho c_{p} \left\{ {\frac{1}{{r_{HR} }} + \frac{s}{{\gamma \left( {r_{q} + r_{s} } \right)}}} \right\}}}$ is the big-leaf thermal time constant, $\rho^{*}$ is the average density of a leaf, and $c_{p}^{*}$ is the average leaf heat capacity. In this case, l^*, the total leaf thickness, is interpreted to be an estimate of the thickness of the thermally-active big-leaf layer. From Fig. 3, we see that T_leaf increases more rapidly than does T_air, indicating $\tau_{leaf} \ll \tau_{air}$, so that no singularity should occur in Eq. 9b. We found that the ensemble fluxes vary only about 2% or less if the Reynolds decomposition (Eq. 5) is made using the fits (Eqs. 8b and 9b) rather than using the smoothed ensemble averages with weighting proportional to the number of realizations.

A small ensemble size is unwelcome, as it introduces larger errors (dotted black lines, Fig. 3), while the outcome can also be strongly affected by large isolated outliers. However, given a sufficient ensemble size, larger asymmetry can be tolerated as such solitary fluctuations become increasingly less important. This is illustrated in Fig. 3, in which the trace is quite smooth near the initial point, but the ensemble averages (green lines) become more erratic as the number of realizations decreases with increasing time since transition.

There are statistical approaches to assess validity that complement these observational considerations. We can evaluate the acceptable order of dissimilarity in realizations to be members of the same ensemble through the “if” Lyapunov (1954) condition: the sequence of independent random variables $\left\{ {X_{i} , i = 1..N} \right\}$, each with finite expected value $\mu_{i}$ and variance $\sigma_{i}^{2}$, will converge to the same normal distribution if for some $\delta > 0$, the following is satisfied,

$$ {\text{Lyapunov}}: \mathop {\lim }\limits_{N \to \infty } \frac{1}{{\left( {\mathop \sum \nolimits_{1}^{N} \sigma_{i}^{2} } \right)^{{{\raise0.7ex\hbox{${\left( {2 + \delta } \right)}$} \!\mathord{\left/ {\vphantom {{\left( {2 + \delta } \right)} 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}}} }}\mathop \sum \limits_{1}^{N} \left( {X_{i} - \mu_{i} } \right)^{{\left( {2 + \delta } \right)}} \to 0. $$

(10)

For the offset removal for the individual realizations: $X_{i}^{off} = X_{i} - \left( {\mu_{i} - \mu_{e} } \right)$ and some small $\alpha > 0$, we obtain the following expression in the third moment of tendencies (see Appendix 2 for derivations),

$$ \mu_{{X^{off} \_ens}}^{3} \propto N^{{\left( {{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}} - \alpha } \right)}} \sigma_{ens}^{3} , $$

(11)

where $\sigma_{ens} = \sqrt {\frac{1}{N}\mathop \sum \nolimits_{1}^{N} \sigma_{i}^{2} }$. Equation 11 is the constraint on the acceptable asymmetry of the tendencies around their ensemble (skewness: $Sk = \frac{{\mu_{{X^{off} \_ens}}^{3} }}{{\sigma_{ens}^{3} }}$) to be considered members of the same ensemble in terms of convergence with using normal distribution (“if” Lyapunov (1954) condition, Eq. 10), which is the indicator of the similarity. We see that if $\alpha > {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}$ the skewness may even reduce when the size of the ensemble increases. Moreover, from Eq. 11 we see that larger standard deviations accommodate a larger spread of the individual realizations in the ensemble (larger third moment), allowing for larger acceptable fluctuations. That the near normality of the tendency fluctuation we identified is a valid assumption for turbulent fluctuations follows from the Millionchikov hypothesis (Monin and Yaglom 1971) from which, for zero skewness, the kurtosis ($Ku = \frac{{\mu_{{X^{off} \_ens}}^{4} }}{{\sigma_{ens}^{4} }}$) should be the Gaussian. This means that the convergence to normality in our case should be adequate for the distributions of fluctuations with small skewness.

We examine turbulent moments during a transient adjustment period, uncommon for field eddy-flux measurements, but routine in laboratory and, recently, in large-eddy modelling studies. Recall that the skewness of the turbulent fluctuations directly reflects a scalar transport asymmetry owing to the turbulent transport of variance (Wyngaard and Weil 1991).

2.4 Ensemble Flux

We refer to the sets of samples arranged by time following the light transition with offsets removed as tendencies. After assessing that such tendencies are indeed representative for single ensembles (see Sect. 2.3), we define an ensemble Reynolds flux (Monin and Yaglom 1971; Wyngaard 2010). For each environmental variable, for the ensemble of tendencies $\left\{ {X_{i} \left( t \right), i = 1..N} \right\}$, we obtain the statistical ensemble of their fluctuations $\left\{ {X_{i}^{\prime} \left( t \right), i = 1..N} \right\}$ (Eq. 13). Following the conventional approach to calculate Reynolds fluxes, we take the vertical-velocity fluctuation $w_{i}^{\prime} \left( t \right)$ and a second variable (X’), the fluctuation of a conserved scalar, [$T_{i}^{\prime} \left( t \right),c_{i}^{\prime} \left( t \right),q_{i}^{\prime} \left( t \right)$, or longitudinal velocity fluctuation $u_{i}^{\prime} \left( t \right)$], to obtain the eddy-covariance ensemble flux (Eq. 12) at a time t following transition. This approach, common in laboratory fluid experiments, aims to estimate flux covariances in unsteady conditions (e.g., He and Jackson 2000),

$$ F_{X} \left( t \right) = {\langle}w^{\prime}\left( t \right)X^{\prime}\left( t \right) {\rangle}= \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} w_{i}^{\prime} \left( t \right)X_{i}^{\prime} \left( t \right), $$

(12)

$$ X_{i}^{\prime} \left( t \right) = X_{i} \left( t \right) - {\langle}X\left( t \right){\rangle}, $$

(13)

$$ F_{X} \left( t \right)\mathop \to \limits^{d} {\mathbf{N}}\left( {\mu_{F} \left( t \right),{\raise0.7ex\hbox{${\sigma^{2}_{{w^{\prime}X^{\prime}}} }$} \!\mathord{\left/ {\vphantom {{\sigma^{2}_{{w^{\prime}X^{\prime}}} } N}}\right.\kern-0pt} \!\lower0.7ex\hbox{$N$}}} \right),\sigma^{2}_{{w^{\prime}X^{\prime}}} \propto \sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}} . $$

(14)

Equation 12 defines kinematic (${\langle} {w^{\prime } T^{\prime } }{\rangle}$) or sensible heat flux $(H = \rho c_{p} {\langle}w'T'{\rangle}$), net ecosystem exchange $(NEE = {\langle}w'c_{a}^{\prime}{\rangle}$), and evapotranspiration $(E = {\langle}w'q'{\rangle}$) as eddy-covariance ensemble fluxes. We define the convergence to the ensemble flux (Eq. 14), where $\sigma^{2}_{{w^{\prime}X^{\prime}}}$ is the ensemble-flux variance, which can be estimated through the product of the individual variances: $\sigma^{2}_{{w^{\prime}X^{\prime}}} \propto \sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}}$ (Appendix 1). For $N \to \infty ,\sigma^{2}_{{w^{\prime}X^{\prime}}} /N \to 0$, hence, almost everywhere $F_{X} \left( t \right) \to \mu_{F} \left( t \right)$ as above.

The standard deviation of the normal distribution in Eqs. 7 and 14, the standard error of the mean of the realizations, approaches zero with large ensemble size, and the mean should go to a stable continuous function (Monin and Yaglom 1971). For continuous functions, convergence “almost everywhere” is equivalent to $F_{X} \left( t \right) \equiv \mu_{F} \left( t \right)$.

Offset removal changes neither the ensemble average nor the ensemble flux (see Appendix 3 for the flux-covariance derivation),

$$ {\text{additive feature:}}\quad \left\langle {w^{\prime } X^{\prime } } \right\rangle \approx \left\langle {w^{{\left[ {t_{1} ,t_{2} } \right]\prime}} X^{{\left[ {t_{1} ,t_{2} } \right]\prime}} } \right\rangle , $$

(15)

$$ {\text{ergodicity feature:}}\quad \left\langle {w^{{\left[ {t_{1} ,t_{2} } \right]\prime}} X^{{\left[ {t_{1} ,t_{2} } \right]\prime}} } \right\rangle \approx \left\langle {w^{{T^{\prime}}} X^{{T^{\prime}}} } \right\rangle - \left\langle {w^{{T^{\prime}}} } \right\rangle \left\langle {X^{{T^{\prime}}} } \right\rangle . $$

(16)

The difference between the ensemble- and time-averaged flux estimates is the product of two terms $\left\langle {w^{{T^{\prime}}} } \right\rangle \left\langle {X^{{T^{\prime}}} } \right\rangle$, representing the ensemble average of the time-averaged fluctuations. For steady-state conditions, we can substitute the offset with the same length bulk mean removal $\frac{1}{N}\mathop \sum \nolimits_{i = 1}^{N} \left( {X_{i} \left( {t,x} \right) - \overline{{X_{i} }}^{T} } \right) = \frac{1}{N}\mathop \sum \nolimits_{i = 1}^{N} X_{i} \left( {t,x} \right) - \frac{1}{N}\mathop \sum \nolimits_{i = 1}^{N} \overline{{X_{i} }}^{T} \to \left\langle X \right\rangle - \left\langle X \right\rangle = 0$, which leads to the flux equality $\left\langle {w^{{\left[ {t_{1} ,t_{2} } \right]\prime}} X^{{\left[ {t_{1} ,t_{2} } \right]\prime}} } \right\rangle \approx \left\langle {w^{T^{\prime}} X^{T^{\prime}} } \right\rangle$ in Eq. 16. If the duration of the mean removal is small (not equal to) comparing to the offset-removal time interval $\left[ {t_{1} ,t_{2} } \right]$ during changing conditions, $\left\langle {w^{{T^{\prime}}} } \right\rangle \left\langle {X^{{T^{\prime}}} } \right\rangle \ne 0$ in Eq. 16. This indicates a time-flux deficiency for shorter time averaging, corresponding to the well-known low-frequency loss issue addressed in Sect. 1.2 and in the following Sect. 3.3.

Equations 15 and 16 indicate that the ensemble flux can be seen as a proper non-stationary extension of the standard time-averaged flux, while the offset-removal time interval $\left[ {t_{1} ,t_{2} } \right]$ is the direct dynamic equivalent to the mean-removal averaging interval used to obtain time-averaged fluxes. Moreover, Eq. 15 introduces the additive feature of the ensemble flux: if the environmental conditions of the underlying events are not identical, the ensemble-averaged flux is a representative average flux among these events in the same way as is the time-averaged flux, giving us an ability to estimate typical transient responses.

2.5 Practical Procedure to Satisfy Event-Similarity Conditions for the Ensemble-Flux Calculation

We give the procedure used to calculate the ensemble fluxes for the period following a step change in incident radiative flux:

• Following the Kivalov and Fitzjarrald (2018) procedure outlined in Sect. 2.1, we identified shadow-to-light and light-to-shadow conditional events. These are combined into small clear-period based radiative-flux increment bins, aligning by the time of light transition, thus generating similar-by-radiative-flux ensembles of realizations.

• We apply median filters (Velleman and Hoaglin 1981; Press 2007) to filter large single sets of fluctuations when they occur near the beginning and the end of the conditional-event time series, to eliminate possible contamination from cloud onsets and cloud-regime changes.

• We remove mean offsets of the individual realizations by the averaging in the [− 200,700]-s interval to organize them into their respective ensembles and to minimize possible flux-calculation errors within this offset-removal interval (see Rannik et al. 2016). This large offset-removal interval accounts for all fluctuations up to 15 min.

• We apply the Reynolds decomposition to the prepared subsets by removing the ensemble-mean tendencies from the individual realizations, after which we obtain the ensembles of fluctuations (Eq. 13).

• We construct the “longitudinal” (corresponding to each time point following transition) eddy-covariance ensemble-flux calculation (Eq. 12).

• To compensate for the decrease in the number of the realizations with time following transition, we bin 1-s data to form ensemble fluxes calculated to ensure that each resulting flux is represented by a similar number of ensemble-time averaged cohorts, in this case containing ≈ 3000 ensemble-flux covariances, consistent with the number of flux integral time scales typically used in the conventional time-averaged approach (e.g., George 2017). For 1-s grab sampled data, just after the transition there are ≈ 300 realizations available (Sect. 2.2), corresponding to the flux being estimated at 10-s intervals. Eight minutes after the transition, when only ≈ 50 realizations are available, each bin corresponds to flux estimated at a 60-s interval.

3 Results

3.1 Validation of Similarity: Examining Ensemble-Member Distributions by Sampling Above-Canopy Turbulence

To assess whether or not members assigned to an ensemble cohort belong together (see Fig. 3), we examined the histograms of the distributions of the tendencies around their ensemble means obtained at selected times following the shadow-to-light transition (red, green and blue histograms; at 0, 150, and 300 s, Fig. 4). To be consistent with our flux-calculation procedure (Sect. 2.5), each distribution is taken from the set of ≈ 3000 grab samples (corresponding to a 10-s time bin at the shadow-to-light transition, t = 0). The near-normal histograms (Fig. 4), describing the changes in T, w, c_a, and q, were checked using a Kolmogorov–Smirnov test. Results support our argument that these cohorts can be treated as members of their respective ensembles. The distributions are approximated by the normal-distribution probability density function (p.d.f.) curves (red line, t = 0; green line, t = 150 s; and blue line, t = 300 s) with single parametrizations describing any given time bin following transition. We argue that this convergence of the ensembles toward normality confirms the similarity of realizations (see Sect. 2.3).

Fig. 4

Distributions of the tendencies around their ensemble means (from Fig. 3) of T_air (a, b), w (c, d), c_a (e, f), and q (g, h) taken and labelled at 0 (red), 150 (green), and 300 (blue) seconds after the shadow-to-light transition; lines—respective normal-distribution p.d.f. fits; [900, 1100] W m⁻² clear-period S↓ bin; for the HF (left column) and the LBA (right column) sites

The distributions (red) at the time of the shadow-to-light transition show the smallest standard deviations when compared to those taken at 150 and 300 s after transition (green and blue respectively). The increase in the standard deviations is the result of the increasing size of the time bins, having fewer available realizations as the time since transition increases. Nature does not provide us with a uniform number of cases to allow bins equally spaced in time. Note that these distributions are very much like those that might have been obtained using the conventional time-averaging approach, in which conditions are assumed to be steadier. However, these results describe the turbulent state of the layer at a particular moment during rapidly-changing conditions.

Consider the vertical-velocity skewness ($Sk_{w} = \frac{{\mu_{{w^{off} \_ens}}^{3} }}{{\sigma_{ens}^{3} }}$) and kurtosis ($Ku_{w} = \frac{{\mu_{{w^{off} \_ens}}^{4} }}{{\sigma_{ens}^{4} }}$). Over a canopy during convective conditions, Sk_w changes sign from negative to positive just above the canopy–atmosphere interface (e.g., Baldocchi and Meyers 1988; Sakai 2000, Fig. 3.10). Negative above-canopy Sk_w during convective conditions is often observed (Raupach and Thom 1981; Fitzjarrald et al. 1988, 1990; Poggi et al. 2004), a fact sometimes associated with the turbulent transfer of stable air from inside the canopy to canopy top. Further into the roughness sublayer, Sk_w becomes positive, as would be found in the convective plane surface layer. The change in the sign of the skewness implies that it passes through zero where Ku_w should be approximately Gaussian (see the Millionchikov hypothesis, Monin and Yaglom 1971), meaning that the fluctuations can be distributed normally, consistent with Fig. 4.

Non-zero skewness in the T and w distributions (Sk_T and Sk_w, Table 2) is evidenced by the distribution means being slightly shifted from the normal-curve fits (Fig. 4a–d). We found that Sk_T and Sk_w change sign from negative to positive depending on the time following transition. The eventual positive Sk_T and Sk_w are about half the values observed during convective conditions in a plane atmospheric surface layer (e.g., Chiba 1978; Andreas et al. 1998; Liu et al. 2011). At 300 s following transition, Sk_w is comparable to that found in a stationary roughness sublayer (e.g., Rannik et al. 2003). Tampieri et al. (2000) and Maurizi (2006) suggest that the kurtosis should gradually increase from its Gaussian values in a parabolic Ku–Sk relation when the skewness departs from zero, and we can trace the similar kurtosis-increasing behaviour with the skewness significantly shifted from zero for most of the scalar distributions (* in Table 2).

Table 2 Observed standard deviation (σ), skewness (Sk), and kurtosis (Ku) of tendency distributions in w, T, c_a, q from Fig. 4 at 0, 150, and 300 s after the shadow-to-light transition; for the HF and LBA sites

Further evidence of the success of our description of the evolving turbulent layer following the transition is given by the time series of correlation coefficients for the relevant flux variables ($r_{Xw}$, see Appendix 1). We found that the correlation coefficients approach those expected for an equilibrium roughness sublayer. Within 200 s of the rapid onset of light, $r_{Tw} ,$ for example, approaches ≈ 0.5, a value expected for very convective conventional surface layer (e.g., Haugen et al. 1971), while $\left| {r_{uw} } \right|$ remains between 0.3 and 0.4 during the entire transition, values typically seen in the roughness sublayer above a forest. These results give us confidence that, while the abrupt onset of convective conditions stimulated by the light impulse changes the character of turbulence just above the canopy, approximate turbulent equilibrium is established relatively quickly, well before the fluxes approach their equilibrium time-averaged values.

We conclude that, after the Reynolds decomposition is applied to enforce zero means, the resulting fluctuations are the ensemble equivalents to the conventional time-series fluctuations used in time-averaged eddy-covariance flux, except that we can take our results to be valid during the transient conditions following the light transition.

3.2 Canopy (Bulk Leaf)- and Air-Temperature Thermal Time Constants

As the sunlit period begins, the enhanced radiative absorption warms the leaves, and the consequent heating of air over the surface promotes an enhanced buoyancy flux. As turbulent mixing increases, this perturbation heat flux converges into a transient thermal IBL. Both processes are non-stationary during the first few minutes, but eventually reach approximate equilibria. To quantify these changes, we seek the time constants for this adjustment and the equilibria for the temperature and ensemble heat flux.

Applying Eq. 8b to fit the air-temperature ensembles (Fig. 3c, d) and then Eq. 9b to fit the leaf-temperature ensembles (Fig. 3a, b) using the known $T_{eq}^{air} ,T_{0}^{air}$, τ_air for air, we find good fits (blue lines) in both the non-stationary and stationary parts of the ensembles (green lines). Note that in Fig. 3a, b the equivalent blackbody temperature that we refer to as a ‘leaf temperature’ includes a mixture of sunlit and shadowed canopy in the radiometer footprint.

The air temperature (Fig. 3c, d) eventually approximates equilibrium conditions $\frac{1}{{\rho c_{p} }}\left( {\frac{{\partial R_{net} }}{\partial z} + \frac{\partial H}{\partial z}} \right) \approx 0$, a steady state that is reached in the HF-site ensemble within 10 min and in the LBA-site ensemble within 20 min following the radiative-flux jump. This indicates that, after 10–20 min, the flux divergence becomes negligibly small and we return to approximate constant-flux conditions of a conventional surface layer.

The leaf-temperature results feature an initial 30-s sharp increase corresponding to a small time constant τ_leaf, followed by a 10–20 min slow increase as the air temperature approaches equilibrium with a large time constant τ_air, then remaining approximately stationary when the leaf-temperature final state is related to the air-temperature equilibrium with the applied leaf heating rate: $\frac{{r_{HR} }}{{\rho c_{p} }}R_{neti} + T_{eq}^{air}$ (Fig. 3a, b). This yields the time constants τ_air and τ_leaf (Fig. 5). Linear trends in τ_air are associated with a larger transient thermal IBL thickness for larger radiative-flux jumps.

Fig. 5

Air-temperature time constants τ_air (a) and leaf (canopy) blackbody-temperature time constants τ_leaf (b); points—measured values, lines—linear fits; see legend for the colour scheme; for the LBA (two levels: dry season—D; rainy season—R) and the HF growing season—G

It is remarkable that the estimated whole-canopy blackbody-temperature time constants (Fig. 5b), bulk descriptors of an equivalent ‘big leaf’, are strikingly similar to time constants observed for some kinds of individual leaves (Wiegand and Swanson 1972; Schymanski et al. 2013) and those reported by Jones (1992) calculated for simple geometric forms after Monteith (1981): beech (d = 60 mm; l^* = 1 mm) 20–60 s, where d is an effective leaf diameter and l^* is an effective leaf thickness. Jones (1992) reported that “respiring leaves” (r_l = 50 s m⁻¹) have a 30% smaller time constant than that for inactive leaves. However, there is no way a priori to determine which of individual-species’ properties the whole canopy might emulate. This finding is also largely consistent with the early direct measurements of transpiring leaves (Clements 1934; Shull 1936; Ansari and Loomis 1959), who found that such leaves exhibit slightly smaller thermal time constants than do non-transpiring leaves, and their cooling time is correspondingly smaller. Both the Ansari and Loomis (1959) and Jones (1992) results agree with those shown in Fig. 5b, in which τ_leaf values are larger for lower radiative-flux values when transpiration is reduced.

3.3 Dynamic and Equilibrium Ensemble Fluxes

Time- and ensemble-flux estimates converge in the steady state (Sect. 2.4), and we identify the start of these conditions with the point at which the dynamic (short-time event, red circles, Figs. 6, 7) and steady-state (long-time event, blue squares) fluxes converge. Ensemble fluxes were calculated and time averaged by the procedure described in Sect. 2.5. Just after the transition, when the maximum number of realizations are available (dashed black lines), the dynamic ensemble fluxes (red circles) undergo rapid steady growth while in direct sunlight (Figs. 6a, c, e, 7a, c, e) and a reduction when the shadow returns (Figs. 6b, d, f, 7b, d, f), demonstrating the true flux sensitivity to light change. These fluxes are adequately approximated by exponential approach-to-equilibrium curves (green lines) with the characteristic time constants τ₊ for light and τ₋ for shadow/long-term equilibrium respectively (see green marks on Figs. 6, 7).

Fig. 6

For the HF site: ensemble fluxes: (a, b) kinematic heat flux ($\left\langle {w^{\prime } T^{\prime } } \right\rangle$), (c, d) net ecosystem exchange (NEE), and (e, f) evapotranspiration (E) for the short-duration dynamic events (red) and long-duration static events (blue); dashed purple lines—5-min time-averaged fluxes (as in Doughty and Goulden 2008); green lines—exponential parametrizations; standard errors are shown as vertical bars; for the shadow-to-sunlit transition (left plot) and for the sunlit-to-shadow transition (right plot); solid red square—outliers possibly due to after-cloud downdraft flux suppression are seen in all the fluxes; dashed red box—outliers after the offset-removal time period; the number of events is shown by the dashed black line; for the [900, 1100] W m⁻² clear-period S↓ bin

Fig. 7

For the LBA site: ensemble fluxes: a, b kinematic heat flux ($\left\langle {w^{\prime } T^{\prime } } \right\rangle$), c, d net ecosystem exchange (NEE), and e, f evapotranspiration (E) for the short-duration dynamic events (red) and long-duration static events (blue); dashed purple lines—5-min time-averaged fluxes (as in Doughty and Goulden 2008); green lines—exponential parametrizations; standard errors are shown as vertical bars; for the shadow-to-sunlit transition (left plot) and for the sunlit-to-shadow transition (right plot); solid red square—outliers possibly due to after-cloud downdraft flux suppression seen in all the fluxes; the number of events is shown by the dashed black line; for the [900, 1100] W m⁻² clear-period S↓ bin

The initial local maximum observed during the sunlit period in the kinematic heat flux ($\left\langle {w^{\prime } T^{\prime } } \right\rangle$, Figs. 6a, b, 7a, b) is consistent with the flux representation according to Fick’s law, given the difference in time constants for leaf heating (≈ 30 s) and air heating (≈ 200 s) (see Sect. 3.2). Thus, following the shadow-to-light transition, kinematic heat flux rapidly increases along with leaf heating, stabilizes as the air warms, and then decreases somewhat after the air reaches equilibrium and the leaf–air-temperature difference decreases.

For the shadow-to-light transitions, for net ecosystem exchange and evapotranspiration (NEE and E, Figs. 6c, e, 7c, e), transient-flux maxima within first 500 s following transition (the transient overshoot) are clearly present for the dynamic ensemble-averaged fluxes (red). The dynamic- (red) and steady-state (blue) fluxes converge after 700 s, signalling the approach to steady-state conditions. For light-to-shadow transitions, for net ecosystem exchange and evapotranspiration (Figs. 6d, f, 7d, f), a monotonic reduction with convergence of the dynamic- (red) and steady-state (blue) fluxes occurs after 500 s. The dynamic fluxes deviate at longer durations, likely a result of the small number of realizations then available.

Note that the 5-min time-averaged fluxes (dashed purple lines, Figs. 6, 7) underestimate the ensemble fluxes by up to 40% during the first few minutes following transition. They become comparable to the ensemble-averaged fluxes only after steady state is reached, 15–20 min following transition, for both the shadow-to-light and light-to-shadow transitions. These 5-min time-averaged fluxes are comparable to values reported by Doughty and Goulden (2008), who analyzed a much smaller dataset obtained at the LBA km-67 and km-83 sites in the Tapajós National Forest, the former site being one source of our data. Clearly, during the transient period, averaging at 5-min intervals is inadequate to the present task.

When the number of realizations exceeds 100, all the ensemble fluxes are well approximated by the exponential relaxation curves (green lines, see Eq. 8b but here for fluxes). However, as the number of realizations falls below 100 during the first 200 s after the transition, all fluxes simultaneously decrease below these exponential fits. The decrease during the light period (left panels) appears from 230 to 350 s following shadow and could be a response provoked by a downdraft that occurs following the cumulus-cloud passage (Gronemeier et al. 2017). The decrease during the shadow periods (right panels) appears from 200 to 300 s following cloud onset and could be a response on the possible outflows from the developed cumulus clouds (see Garstang and Fitzjarrald 1999, Fig. 6.6c). This makes all these flux-suppression points to be outliers. Future studies are needed to clarify these points.

The successful description of the sensitivity following the light transitions by the exponential approach-to-equilibrium fits (green lines) represents an observational validation of the use of the big-leaf hypothesis to describe the dynamic behaviour of the outer canopy during fluctuating-light conditions. This is the region of the canopy in which most of the radiative flux is absorbed and the bulk of fluxes originate (see Introduction).

3.3.1 Ensemble-Flux Time Constants

By applying the exponential relaxation fits to the best available set of data in the S↓ range for clear periods between 900 and 1100 W m⁻², we found the ensemble-flux time constants as shown in Table 3.

Table 3 Exponential-fit time constants for light-onset (shadow-to-light, Figs. 6 and 7, left) and shadow-onset (light-to-shadow, Figs. 6 and 7, right) ensemble-flux change by the Eq. 8b; where “τ₊” represents the initial dynamic-flux change approaching to the transient-overshoot equilibrium, and “τ₋” refers to the approach to a steady-state long-time equilibrium; clear-period S↓ in the range [900, 1100] W m⁻²

The evapotranspiration bulk-opening time constant, τ₊, exceeds that of the typical stomatal-opening time constant for dominant species—oak and maple, determined from single-leaf studies at the HF site (≈ 3 min, Woods and Turner 1971). This suggests that the slow increase in evapotranspiration results from stomatal opening after the radiative flux intensifies, rather than a result of the rapid change in vapour pressure deficit that closely follows T_leaf. The evapotranspiration bulk-decay time constant, τ₋, after the local maximum for the HF site is shorter than is typical for the closing value for these species (≈ 5 min, Woods and Turner 1971).

Net ecosystem exchange and evapotranspiration are linked through the stomatal conductance. However, the net-ecosystem-exchange bulk-opening time constant is ≈ 40% smaller than is the corresponding evapotranspiration time constant, perhaps a result of dynamic decoupling of net ecosystem exchange and evapotranspiration (and hence, stomatal conductance) in fluctuating-light conditions, as suggested by McAusland et al. (2015, 2016).

We emphasise that all time constants in Table 3 are about five times larger than characteristic turbulent-mixing time constants (30–40 s; see Fig. 15, Appendix 1), indicating that the observed flux dynamics are largely independent of transient changes in turbulence intensity. This is also shown by the fact that c_a and q stay largely unchanged during and after the transition (see Fig. 3e–h), signalling an absence of the significant changes in flux convergence of these quantities in the lowest layers.

3.3.2 Durations of Transient Overshoots (Sect. 2.2) and Cumulative Radiative Dosages

The transient overshoots (red circles, Figs. 6, 7), when the dynamic fluxes exceed steady-state equilibrium values, exhibit consistent behaviour, generally increasing in duration when the radiative flux decreases (red, green, and blue lines, Fig. 8). The transient-overshoot duration (the regime-change point on the green line, Figs. 6, 7, left panels) is obtained by identifying the intersection of the fitted transient exponential curve with the respective flux curve as the flux begins consistently to relax toward the long-term equilibrium.

Fig. 8

Duration of the maximal transient responses (transient overshoots) for $\left\langle {w^{\prime } T^{\prime } } \right\rangle$ or H (red line), NEE (green line), and E (blue line) and the cumulative radiative dosage received under the transient-overshoot periods for NEE (green area) and E (blue area); for the HF (left column) and the LBA (right column) sites

In Sect. 3.3, we found that the overshoots occur after a sufficiently long relaxation period in shadow, during which radiative-flux levels remain below the light-saturation threshold (PAR ≈ 1000 μmol m⁻² s⁻¹; Critchley and Russell 1994) and can be interpreted as the full capacity physiological response of plants after this recovery occurs (see Sect. 2.2 for discussion).

This hypothesized link between the transient overshoots and a physiological response becomes more plausible when we estimate the cumulative dosage of radiative flux (green for NEE and blue for E, Fig. 8), found by integrating over the transient-overshoot duration. This is possible because the radiative-flux jump occurs almost instantaneously compared to the slower responses of net ecosystem exchange or evapotranspiration (see Appendix 4 for details). The levels of cumulative dosages stay about the same for different S↓ values, decreasing only in the afternoon. The dosage net-ecosystem-exchange maxima practically coincide with the most prevalent cloudy-day S↓ values at both HF (≈ 900–950 W m⁻²) and LBA (≈ 1000 W m⁻²) sites (see Fig. 1e, f for maximal event counts). We speculate that this results from a plant physiological adaptation to preferred radiative-flux levels during fluctuating-light regimes.

3.4 Transient Forest Response to the Step Change in Incident Radiative Flux

3.4.1 Sensible (Kinematic) Heat Flux, Net Ecosystem Exchange, Evapotranspiration, and Water-Use Efficiency

Combining the dynamic sunlit and shadow fluxes (Sect. 3.2, Figs. 6, 7) for the kinematic heat flux ($\left\langle {w^{\prime } T^{\prime } } \right\rangle$), net ecosystem exchange (NEE), and evapotranspiration (E) for different radiative-flux bins with their time constants (Sect. 3.2.1, Table 3), and the transient-overshoot durations (Sect. 3.2.2, Fig. 8), we visualize the transient ecosystem response to radiative-flux change for the HF growing season (Fig. 9a, c, e) and the combined LBA wet and dry seasons (Fig. 9b, d, f).

Fig. 9

Composite dynamic ensemble fluxes following transition at different values of photosynthetically-active radiative flux (PAR): a, b kinematic heat flux ($\left\langle {w^{\prime } T^{\prime } } \right\rangle$), c, d net ecosystem exchange (NEE), and e, f evapotranspiration (E). The periods representing time since the shadow-to-light transition are in the top portion of each panel, marked “in sun”. These were preceded by about 300-s long shadow periods. The responses during the light-to-shadow transition are in the lower portion of each panel, marked “in shadow”. These were preceded by an approximately 300-s long sunlit period. Dashed lines are the most common clear-period light intensity (red; vertical axis) and duration (purple; horizontal axis) on partly cloudy days for these forests (Kivalov and Fitzjarrald 2018); for the HF (left column) and the LBA (right column) sites

The sensible (kinematic) heat flux for both HF and LBA sites (Fig. 9a, b) shows behaviour consistent with a big-leaf model, reducing synchronously with the radiative flux during both the sunlit (“in sun”) and “in shadow” periods.

Net ecosystem exchange during the sunlit period for the HF site (“in sun”, Fig. 9c) shows maximum assimilation at about 500 s after transition and PAR ≈ 1800 μmol m⁻² s⁻¹, with the subsequent reduction for PAR > 2000 μmol m⁻² s⁻¹. Net ecosystem exchange during the sunlit period at the LBA site (“in sun”, Fig. 9d) shows a consistently high assimilation band at about 500 s after transition in the PAR > 1800 μmol m⁻² s⁻¹ region.

Evapotranspiration during the sunlit period at both HF and LBA sites (“in sun”, Fig. 9e, f) features a less distinct, more broadly spread, peak about 500–600 s after transition. This could be associated with the fact that the water vapour flux includes both forest transpiration and direct evaporation. Evapotranspiration maxima occur at the larger photosynthetically-active radiative flux values than do the respective net-ecosystem-exchange maxima. For cloud-shadow events (“in shadow”, Fig. 9), there is a consistent decrease in all fluxes when the photosynthetically-active radiative flux decreases.

The instantaneous water-use efficiency is usually defined as the ratio between assimilation and transpiration (e.g., Medrano et al. 2015). Assimilation exceeds net ecosystem exchange, because the latter includes respiration (see Kirschbaum et al. 2001). Transpiration is lower than evapotranspiration, which includes other evaporation terms, which are not easy to separate. Here we approximated the whole-canopy water-use efficiency as $WUE = {\raise0.7ex\hbox{${NEE}$} \!\mathord{\left/ {\vphantom {{NEE} E}}\right.\kern-0pt} \!\lower0.7ex\hbox{$E$}}$. Though this definition differs somewhat from that used in Medrano et al. (2015), we argue that the near-constant forest-respiration rates (≈ 8 µmol m⁻² s⁻¹ for the HF site; ≈ 12 µmol m⁻² s⁻¹ for the LBA site) and the small change in evaporation over short time scales allow for substantially correct estimates of changes in the instantaneous water-use efficiency.

At both the HF and LBA sites, water-use efficiency exhibits large variations for a given photosynthetically-active radiative flux range (Fig. 10a, b), generally larger in the periods of high net ecosystem exchange during the transient response to the light change (see Fig. 9c, d). That the sunlit water-use efficiency at the HF site (Fig. 10a) is larger than the water-use efficiency at the LBA site (Fig. 10b) is consistent with the higher net ecosystem exchange at the HF site. The shadow water-use efficiency (low sub-panel) is about twice that of the sunlit water-use efficiency (upper sub-panel) for both HF and LBA sites, due to the low evapotranspiration during shadow. The increase of shadow water-use efficiency with the photosynthetically-active radiative flux (Fig. 10a, b) is consistent with the shadow net ecosystem exchange increase with photosynthetically-active radiative flux (Fig. 9c, d), and can be explained by its location on the light-limited branch of the photosynthesis curve, while the shadow evapotranspiration remains relatively low (Fig. 9e, f) at all times.

Fig. 10

Composite water-use efficiency ($WUE = {\raise0.7ex\hbox{${NEE}$} \!\mathord{\left/ {\vphantom {{NEE} E}}\right.\kern-0pt} \!\lower0.7ex\hbox{$E$}}$) following transition at different photosynthetically-active radiative flux (PAR) values. Dashed lines are the most common clear-period light intensity (red; vertical axis) and duration (purple; horizontal axis) on partly-cloudy days for these forests (Kivalov and Fitzjarrald 2018); for the HF (left column) and the LBA (right column) sites

We found a preferred period at approximately 500 s following the shadow-to-light transition for which there are joint maxima of net ecosystem exchange and evapotranspiration. This feature appears more strongly at the HF site than at the LBA site, but it could be a robust feature of the outer-canopy system response. The vertical dashed lines in Fig. 9 and 10 indicate the average duration of sunlit intervals on partly-cloudy days at each site; the horizontal dashed line shows the average final photosynthetically-active radiative flux after the jump (Kivalov and Fitzjarrald 2018). Clearly, in the partly-cloudy conditions we studied, carbon uptake does not often reach its potential maximum, and therefore light usage cannot quite be optimal during dynamic photosynthesis.

3.4.2 Dynamics of Derived Bulk Aerodynamic and Canopy Conductances and Active Thermal-Layer Thicknesses (the Inverse Problem Solutions)

In Sects. 2.2 and 3.3, we noted that the average scalar concentration and flux behaviour can be approximated by exponential relaxation functions, featuring equivalent time constants. Thus, the time constants obtained from the outer-canopy temperature (Sect. 3.2) and from the flux ensembles (Sect. 3.3) directly correspond to a big-leaf formulation being adequate to describe the transient events. To define the thermally-active big-leaf layer parameters, we apply the bulk approach to the ensemble fluxes (Eqs. 17 and 18) from which we find resistances r_H and $\left( {r_{q} + r_{c} } \right)$ to kinematic heat flux ($\left\langle {w^{\prime } T^{\prime } } \right\rangle$) and evapotranspiration ($\left\langle {w^{\prime } q^{\prime } } \right\rangle$), respectively.

Examining the tower profile measurements, we note that the specific humidity of air (q_a) is nearly constant, not significantly varying with height above the canopy (see also Fitzjarrald et al. 1990, Fig. 1), implying that $q_{a} \approx q_{air}$. On the contrary, T_a, the average temperature within a transient thermal IBL with thickness Δz above the canopy, does not necessarily equal the tower-measured air temperature T_air. Knowing the vertical gradient of the tower − surface temperature, $\varGamma = \frac{{T_{leaf} - T_{air} }}{h}$, where h is the distance from the outer canopy to the tower measurement-site location (taking h = 7 m for the HF site, and h = 14 m for the LBA site), we can re-evaluate the flux equations for the variable-thickness (Δz) of the transient thermal IBL as follows,

$$ {\langle}w'T'{\rangle} = \frac{1}{{r_{H} }}\left( {T_{leaf} - T_{a} } \right) = \frac{1}{{r_{H} }}\varGamma \frac{\Delta z}{2}, $$

(17)

$$ {\langle}w'q^{\prime}{\rangle} = \frac{1}{{\left( {r_{q} + r_{c} } \right)}}\left( {q_{leaf} - q_{a} } \right). $$

(18)

To estimate the thermally-active big-leaf layer thicknesses (l^*) and the transient thermal IBL thickness (Δz), we use expressions for τ_leaf and τ_air for the leaf-(canopy) and air-temperature time constants from Sect. 3.2. Because $\rho^{*} c_{p}^{*}$ ≈ 2.7 × 10⁶ J K⁻¹ m⁻³ for a leaf (Jones 1992) and $\rho c_{p}$ ≈ 1.231 × 10³ J K⁻¹ m⁻³ for air, $\rho^{*} c_{p}^{*} > 2 \times 10^{3} \rho c_{p}$, and we can avoid including the air heat capacity within the canopy in the heat budget, unless the air volume within the canopy is ≈ 10³ times larger than the leaf volume. The thermally-active outer-canopy layer can be estimated to be 3–5 m thick (70–90% of radiative flux is absorbed, see Introduction). However, because Δ_z ≫ 5 m, we can also consider this outer-canopy air as part of the transient thermal IBL.

Combining expressions for τ_leaf and τ_air from Sect. 2.2 (see Appendix 4) with Eqs. 17 and 18 evaluated for the whole canopy, we obtain Δz and l^* as,

$$ \Delta z = 2\frac{{\tau_{air} }}{{r_{HR} }} = 2\tau_{air} \left( {\left( {\frac{{4\sigma T_{air}^{3} }}{{\rho c_{p} }}} \right) + \frac{{{\langle}w^{\prime}T^{\prime}{\rangle}}}{{\varGamma \frac{\Delta z}{2}}}} \right), $$

(19a)

$$ l^{*} = \tau_{leaf} \frac{{\rho c_{p} \left\{ {\frac{1}{{r_{HR} }} + \frac{s}{{\gamma \left( {r_{w} + r_{c} } \right)}}} \right\}}}{{\rho^{*} c_{p}^{*} }} = \tau_{leaf} \frac{{\rho c_{p} }}{{\rho^{*} c_{p}^{*} }}\left( {\frac{\Delta z}{{2\tau_{air} }} + \frac{s}{\gamma }\frac{{{\langle}w^{\prime}q^{\prime}{\rangle}}}{{\left( {q_{leaf} - q_{a} } \right)}}} \right). $$

(19b)

Equation 19a leads to a quadratic equation for Δz with the single positive solution,

$$ \Delta z = \tau_{air} \left( {\frac{{4\sigma T_{air}^{3} }}{{\rho c_{p} }}} \right) + \sqrt {\left( {\tau_{air} \left( {\frac{{4\sigma T_{air}^{3} }}{{\rho c_{p} }}} \right)} \right)^{2} + 4\tau_{air} \frac{{{\langle}w^{\prime}T^{\prime}{\rangle}}}{\varGamma }} . $$

(20)

Taking $r_{q} \approx r_{H}$, a common assumption for turbulent conditions (Monteith and Unsworth 1990), the canopy resistance $r_{c} = \left( {r_{q} + r_{c} } \right) - r_{H}$, noting (Fig. 11) that:$g_{H} = 1/r_{H} ,g_{q + c} = 1/\left( {r_{q} + r_{c} } \right),g_{c} = 1/r_{c}$.

Fig. 11

Evolution of conductances to heat (g_H) and to water-vapour transport (g_q) and calculated from them canopy conductance (g_c); see legend for photosynthetically-active radiative flux (PAR) values for each line; for the HF (left column) and the LBA (right column) sites

That the thickness of the transient thermal IBL at the HF site is about two times smaller than that at the LBA site (Fig. 12a, b) is the result of the difference in τ_air for these sites (Fig. 5). The smaller air heat resistance at the LBA site also leads to the higher heat exchange with larger Δz.

Fig. 12

Temporal evolution of transient thermal IBL thickness (Δz) above the canopy and thickness of the thermally-active big-leaf layer (l*); see legend for photosynthetically-active radiative flux (PAR) values for each line; for the HF (left column) and the LBA (right column) sites

One may question whether or not the blackbody temperature represents the true leaf temperature—the field of view contains leaves in full sun and those in shadow. However, the inferred properties of our ‘big leaf’ could plausibly describe real leaves. Under the Jones’ (1992) assumptions of the average leaf density ρ^* ≈ 700 kg m⁻³ and 90% water content for fresh leaves, the thicknesses of the thermally-active big-leaf layer (l^*) is about the same at both HF and LBA sites, around 1 mm most of the time. Using l^* = 1 mm (Jones 1992) yields the thermally-active outer-canopy layer to be equivalent to one LAI unit. However, if we apply the Vile et al. (2005, Table 2) estimates for shrub and tree leaf thicknesses l^* ≈ 0.42 mm, we obtain more realistic l^* values equivalent to a depth into the canopy of approximately two LAI units during the sunlit periods.

The maximal stomatal conductances $g_{s} \approx g_{c} /LAI \approx g_{c} /2$ (Fig. 11e, f) are close to the corresponding leaf-scale measurements when stomata are open (see Jones 1992). After 600–700 s there is some reduction in g_s. A similar bump and subsequent reduction can be seen in individual leaf measurements (e.g., Fay and Knapp 1996; Kaiser et al. 2017). However, this effect accounts for only about 10% of the stomatal closure and cannot explain the 20–25% dynamic reduction in net ecosystem exchange (Fig. 9c, d) that occurs a little earlier, after 10 min. This agrees with the hypothesis that the net ecosystem exchange and stomatal conductance are uncoupled in the transient regime (McAusland et al. 2015; 2016). Comparing Fig. 9 for NEE with Fig. 11e, f for g_c, we see that the stomata respond by partially closing, and this effect is slightly lagged after the net-ecosystem-exchange reduction to reduce evapotranspiration with its corresponding water loss.

4 Discussion: Photosynthesis Responds Directly to the Light Change Versus Diffuse Radiative Flux Increases Leaf Illumination

In this section we estimate the effect of abrupt light change on net ecosystem exchange and water-use efficiency during realistic partly-cloudy conditions (Kivalov and Fitzjarrald 2018) given sufficient long shadow periods to allow equilibrium conditions to obtain before the light returns (Sect. 4.2). The basis for comparison is provided by defining the steady-state whole-canopy photosynthesis curves (Sect. 4.1). This effort aims to assess the hypothesis that the enhancement in net ecosystem exchange and water-use efficiency due to light change could be a viable alternative to the diffuse-radiative-flux enhancement hypothesis and also helps further to validate the use of a single-layer big-leaf model in fluctuating light conditions.

4.1 Single Steady-State Whole-Canopy Photosynthesis Curves Compared to Time-Averaged Fluxes

If one hypothesizes that the forest responds directly to the incident radiative-flux (S↓) or to the photosynthetically-active radiative flux (PAR) rather than to an increase in the diffuse radiative flux, we suggest that a single whole-canopy photosynthesis curve (solid blue lines, Fig. 13a, b) describes the steady-state whole-canopy NEE—S↓ (or PAR) dependence. In Fig. 13, the combined events were separated by “Day”—shadow periods from 0800 to 1600, “Noon”—light periods from 1000 to 1400, and “Afternoon”—light periods mostly after 1400–1600 given the paucity of morning events (see Fig. 1g, h for event distributions). Since “Day” and “Noon” are the mostly uniformly representative, they are naturally connected by the solid blue line, representing the single daytime whole-canopy photosynthesis curve. The “Afternoon” data (red squares) show the NEE dip below the single photosynthesis curve (blue line), outlining the well-known afternoon reduction in photosynthesis.

Fig. 13

a, b Net ecosystem exchange (NEE) and c, d evapotranspiration (E) during full light (red) and under shadow (blue) for equilibrium, i.e., 15–20 min after transition (squares), and at maximal transient response (circles); dashed blue lines—shadow NEE and E values shifted to the corresponding light-period radiative flux for references; solid blue line (a, b)—single photosynthesis curve combining shadow and light NEE values; “Day”—shadow events 0800–1600, “Noon”—light midday events 1000–1400, “Afternoon”—light events mostly after 1400 since there are few events present before 1000 in the morning (see Fig. 1g, h); for the HF (left column) and the LBA (right column) sites

Net ecosystem exchange and evapotranspiration steady-state results (squares, Fig. 13) are similar to those found using time-averaged fluxes reported for the same sites or for similar regions (Amthor et al. 1994, Fig. 1a; Goulden et al. 2004, Fig. 4; Rocha et al. 2004, Fig. 1; Hutyra 2007; Hutyra et al. 2007, Fig. 3.9).

Goulden et al. (2004) and Hutyra (2007), Hutyra et al. (2007) reported temporally-averaged eddy-flux observations from the LBA km-83 and km-67 tower sites, respectively. Their net ecosystem exchange, along with the morning–afternoon difference discussed earlier, is close to our net ecosystem exchange for the LBA site (red squares, Fig. 13b) between 7 and 15 μmol m⁻² s⁻¹.

Amthor et al. (1994) reported measurements from the same HF site. They found the similar upward NEE—PAR trend we see during the clear high radiative-flux periods (both red squares and blue line, Fig. 13a) with similar values between 10 and 20 μmol m⁻² s⁻¹ for the steady state.

Rocha et al. (2004) report measurements from northern hardwood forest in upper Great Lakes Region, USA, a forest similar to the HF site. In agreement with Rocha et al., we observe net ecosystem exchange increasing rapidly with radiative flux during the cloudy periods (light-dependent branch of the photosynthesis curve) and net ecosystem exchange increasing more slowly during the higher radiative-flux clear periods (light-saturated branch of the photosynthesis curve). Under partly-cloudy skies, the net ecosystem exchange lies between the result for clear periods and cloudy periods increasing with radiative flux. However, Rocha et al. (2004) observed a strong upward tendency in net ecosystem exchange during partly-cloudy conditions that cannot be explained by reference to the familiar rectangular hyperbola, the Michaelis–Menten fit (Thornley 1976). Because such partly-cloudy periods are often associated with cumulus clouds (see Introduction; Kivalov and Fitzjarrald 2018), this result can be explained by accounting for the impact of the transient increase in net ecosystem exchange (red circles, Fig. 13a, b) in the time-averaged fluxes. At the HF site, the transient increase in net ecosystem exchange exactly coincides with PAR ≈ 1200–1800 μmol m⁻² s⁻¹ (or S↓ ≈ 600–900 W m⁻²) and with the increase in net ecosystem exchange during partly-cloudy conditions reported by Rocha et al. (2004). The transient increase in net ecosystem exchange at the HF site also follows the ascending branch pattern, similar to results for partly-cloudy conditions reported by Rocha et al. (2004), with similar values of photosynthesis, increasing from > 10 to < 30 μmol m⁻² s⁻¹.

4.2 Diagnosing the Perturbation Forest Net Ecosystem Exchange and Water-Use Efficiency from the Transient Forest Responses

Can the response to a variable-light environment explain the increase in net ecosystem exchange observed during partly-cloudy conditions? Following Tomson (2014, Fig. 2), we note that during shadow periods largely diffuse radiative flux occurs ($S\!\!\downarrow_{d\_s}$), but the sunlit periods include both diffuse and direct radiative fluxes ($S\!\!\downarrow_{g\_l}$) with about the same diffuse radiative flux as did the preceding shadow period. This allows us to estimate the magnitude of the diffuse radiative flux during the light periods of the conditional events. To be consistent with our event selection, we impose the following constraints:

• The radiative-flux jump between shadow and light periods ($S\!\!\downarrow_{g\_l} - S\!\!\downarrow_{d\_s}$) is at least 300 W m⁻².

• The shadow-period radiative flux naturally occurs and selected within the diffuse radiative-flux related ascending branch of the photosynthesis curves (solid blue lines, Fig. 13a, b): $S\!\!\downarrow_{d\_s}$ < 500 W m⁻².

• The light-period radiative flux naturally occurs and selected within the global (direct and diffuse) radiative-flux related light-saturating branch of the photosynthesis curves (solid blue lines, Fig. 13a, b): $S\!\!\downarrow_{g\_l}$ > 500 W m⁻².

Taking different durations of light (t₊) and shadow (t₋) periods during radiative-flux fluctuations, we estimate different average radiative forcing ($\overline{S\!\!\downarrow }$, Eq. 21a) affecting ecosystems as well as the value of cloudiness (Cld, Eq. 21b) and diffuse fraction (f_d, Eq. 21c),

$$ \overline{S\!\!\downarrow } = \frac{{S\!\!\downarrow_{g\_l} t_{ + } + S\!\!\downarrow_{d\_s} t_{ - } }}{{\left( {t_{ + } + t_{ - } } \right)}}, $$

(21a)

$$ Cld = \frac{{t_{ - } }}{{\left( {t_{ + } + t_{ - } } \right)}}, $$

(21b)

$$ f_{d} = \frac{{S\!\!\downarrow_{d\_s} }}{{\overline{S\!\!\downarrow } }}. $$

(21c)

To diagnose the consequences of cloudiness on the fluxes, we examined a larger set of conditional events (t₋ > 200 s and t₊ < 600 s), which largely overlaps with the initial set used in this work. This choice is justified by looking at the dynamic forest fluxes. We see that

• About 70% of the rapid flux changes occur during the first 200 s following transition (Figs. 6, 7), and these can be approximated by exponential fits (green lines, Figs. 6, 7) using the diagnosed characteristic time constants (Table 3);

• The transient response time for net ecosystem exchange and evapotranspiration is limited to 600–700 s for the HF site and 600 s for the LBA site (Fig. 8).

In conclusion, we note that if we select conditional events with the shadow period > 200 s and the light period < 600 s:

• The shadow-period fluxes are close to the shadow equilibria at the end of the shadow periods for any selected event, which defines the equilibrium starting point for the light-period fluxes.

• The light-period fluxes are within the transient-response ascending branches of the fluxes reaching toward transient maxima for any selected event.

The first conclusion ensures that for these selected conditional events we can examine idealized periodic solutions to estimate the real fluxes. The second conclusion allows us to estimate the natural transient forest response by using the simple Ooba and Takahashi (2003) periodic solution for the following first-order system (written for NEE),

$$ NEE_{{t_{ + } }}^{ + } = NEE_{eq}^{ + } - \left( {NEE_{eq}^{ + } - NEE_{{t_{ - } }}^{ - } } \right){ \exp }\left\{ { - \frac{{t_{ + } }}{{\tau_{ + } }}} \right\}, $$

(22a)

$$ NEE_{{t_{ - } }}^{ - } = NEE_{eq}^{ - } - \left( {NEE_{eq}^{ - } - NEE_{{t_{ + } }}^{ + } } \right){ \exp }\left\{ { - \frac{{t_{ - } }}{{\tau_{ - } }}} \right\}, $$

(22b)

where Eq. 22a describes the light-period dynamic fluxes with equilibrium $NEE_{eq}^{ + }$ and associated time constant τ₊, and Eq. 22b describes the shadow-period dynamic fluxes with equilibrium $NEE_{eq}^{ - }$ and associated time constant τ₋. Following Ooba and Takahashi (2003, Eq. 2), the total cumulative solution evaluated for NEE (or H or E) over the event is

$$ \begin{aligned} NEE_{total} & = \left[ {NEE_{eq}^{ + } t_{ + } + NEE_{eq}^{ - } t_{ - } } \right] \\ & \quad + \left( {\tau_{ - } - \tau_{ + } } \right)\left( {NEE_{eq}^{ + } - NEE_{eq}^{ - } } \right)\frac{{\left( {1 - \exp \left\{ { - \frac{{t_{ + } }}{{\tau_{ + } }}} \right\}} \right)\left( {1 - \exp \left\{ { - \frac{{t_{ - } }}{{\tau_{ - } }}} \right\}} \right)}}{{1 - \exp \left\{ { - \frac{{t_{ + } }}{{\tau_{ + } }}} \right\}\exp \left\{ { - \frac{{t_{ - } }}{{\tau_{ - } }}} \right\}}}, \\ \end{aligned} $$

(23)

while the average value $\overline{NEE} = NEE_{total} /\left( {t_{ + } + t_{ - } } \right)$. To evaluate the absolute increase in the dynamic $\overline{NEE}$ values above the clear-period steady-state values, we use the following ratio with the values from the photosynthesis curves $NEE_{eq} = NEE_{\text{eq}} \left( {S\!\!\downarrow_{g\_l} } \right)$ (solid blue lines, Fig. 13a, b),

$$ \frac{{\overline{NEE} }}{{NEE_{eq} }} = \frac{{NEE_{total} }}{{\left( {t_{ + } + t_{ - } } \right){NEE}_{\text{eq}} \left( {S\!\!\downarrow_{g\_l} } \right)}}. $$

(24)

To evaluate average water-use efficiency ($\overline{WUE}$) over the event, we write,

$$ \overline{WUE} = \frac{{NEE_{total} }}{{E_{total} }}. $$

(25)

By separating the $\overline{NEE} /NEE_{eq}$ ratios (Eq. 24) into average radiative-flux bins ($\overline{S\!\!\downarrow }$, Eq. 21a): [1000 +], [900, 1000], [800, 900], [700, 800], [600, 700], [500, 600], we trace the maxima for different radiative-flux values (Fig. 14a–d). For the HF site, the $\overline{NEE} /NEE_{eq}$ maximum lies in the [700, 800] $\overline{S\!\!\downarrow }$ bin and is about 20% above the clear-period equilibrium values. For the LBA site, the $\overline{NEE} /NEE_{eq}$ maximum lies in the [800, 900] $\overline{S\!\!\downarrow }$ bin and is about 10% above the clear-period equilibrium values.

Fig. 14

Separation of the absolute increase in $\overline{NEE} /NEE_{eq}$ (Eq. 24) and $\overline{WUE}$ (Eq. 25) by the average radiative-flux levels ($\overline{S\!\!\downarrow }$, see legend): a, b$\overline{NEE} /NEE_{eq}$ by cloudiness; c, d$\overline{NEE} /NEE_{eq}$ by diffuse fraction; e, f$\overline{WUE}$ by cloudiness; g, h$\overline{WUE}$ by diffuse fraction; i, j$\overline{NEE}$ (from Eq. 23); k, l$\bar{E}$ (from Eq. 23 evaluated for E); dots—individual events; coloured lines—fitted third-degree polynomials to illustrate trends (following Gu et al. 1999); for the HF (left column) and the LBA (right column) sites

For water-use efficiency (Fig. 14e–h), we see a gradual increase in water use efficiency with increasing Cld and f_d values. This is consistent with:

• The observed large difference in evapotranspiration values, with under-cloud shadow values about four times smaller than those during clear periods (Fig. 13c, d);

• The small difference in net ecosystem exchange, with under-cloud values located on the ascending branch of the photosynthesis curve while the clear-period values located on its light-saturated branch (Fig. 13a, b).

The two forests show distinct separation of water-use efficiency maxima by Cld and f_d. At the HF site, water-use efficiency increases with the reduction in $\overline{S\!\!\downarrow }$, which for the same Cld and f_d implies that the HF forest is more water efficient at a lower incident radiative flux. However, at the LBA site, there are water-use efficiency maxima in the [800, 900] and [700, 800] $\overline{S\!\!\downarrow }$ bins, suggesting that the LBA forest has higher water-use efficiency at high incident radiative flux, perhaps an adaptation in the tropics. Due to the larger net ecosystem exchange (Fig. 14i, j) accompanied by the similar evapotranspiration (Fig. 14k, l), the HF site exhibits about 50% more water-use efficiency than does the LBA site, consistent with the water-use efficiency sensitivity (Fig. 10).

Net ecosystem exchange and water-use efficiency maxima are evident for both Cld (Fig. 14a, b, e, f, i, j) and f_d (Fig. 14c, d, g, h). The cloudiness in the selected events is greater than 25% due to $t_{ - } \ge 200 s$, which may obscure possible further net-ecosystem-exchange and water-use efficiency maxima for smaller Cld and shorter t₋.

5 Summary and Conclusions

To determine the sensitivity of the forest response in a variable-light environment often imposed by cumulus clouds, we presented a practical methodology using an ensemble-flux approach to obtain valid eddy-covariance flux estimates at short time intervals. We found (Sect. 3.3) dynamic ensemble values of kinematic ($\left\langle {w^{\prime } T^{\prime } } \right\rangle$) or sensible heat flux ($H = \rho c_{p} \left\langle {w^{\prime } T^{\prime } } \right\rangle$), net ecosystem exchange ($NEE = \left\langle {w^{\prime } c_{a}^{\prime } } \right\rangle$), and evapotranspiration ($E = \left\langle {w^{\prime } q^{\prime } } \right\rangle$) or latent heat flux (λE). With 300-realization ensembles grab sampled at 1-s intervals (HF site) and 500-realization ensembles grab sampled a 5-s intervals (LBA site) available near transition (Sect. 2.2), we resolved dynamic ensemble fluxes on 10-s and 20-s time scales respectively with ≈ 5% standard error (Sect. 3.3, Fig. 6, 7). These resolution scales are about 15–20 times (HF site) and 6–10 times (LBA site) smaller than the ensemble-flux time constants obtained (150–200 s), ensuring adequate representativeness of the flux results obtained. This also shows that 1-s sub-sampling is preferable to the 5-s sub-sampling for determining dynamic ensemble fluxes and their time-dependent features with higher resolutions.

We justified the convergence of the ensemble flux to a valid mean in terms of the properties of a normal distribution given similar temporal changes in the realizations following similar light changes. We defined sufficient conditions for this convergence by evaluating distribution skewness (Sect. 2.3, Eq. 11 and Sect. 3.1). We showed that the offset removal does not affect the ensemble-flux convergence, deriving the ensemble-flux additive property ensuring that the typical average fluxes can be obtained for non-identical conditions (Sect. 2.4, Eqs. 15 and 16). We found that sufficiently long after the transitions (steady-state conditions), the ensemble-flux converges to the standard time-averaged flux, confirming both the method applicability and the ergodic hypothesis (Sect. 3.3).

Applying the ensemble-average technique to the changes in leaf and air temperatures, we obtained the leaf- and air-warming time constants (τ_leaf and τ_air and their dependence on the radiative-flux jump (Sect. 3.2). These τ_leaf estimates, characteristic of the whole canopy, are close to the time constants obtained from the individual leaf measurements (Shull 1936; Wiegand and Swanson 1972; Jones 1992; Schymanski et al. 2013).

Applying the ensemble-flux technique to the time series of T, c_a, q, and w fluctuations:

• We distinguished the transient responses of the ecosystems by the clear-period defined values of S↓ or PAR and showed that there are previously unidentified transient maxima (overshoots) in net ecosystem exchange above the steady-state values at about 200–600 s following the light onset (Sect. 3.3). These maxima then relax toward steady-state values. As the light alternates with shadow, the evapotranspiration time constants (Sect. 3.3.1, Table 3) resemble stomatal opening and closing time constants previously observed for individual leaves for certain species, with the evapotranspiration time constant 30–40% larger than that for net ecosystem exchange.

• The big-leaf heat and water-transport conductances (g_H and g_q) obtained for changes at short time intervals following the shadow-to-light transition show similar sensitivity to abrupt light changes at both forests, with slightly higher g_H and lower g_q at the LBA site. The derived stomatal conductances (g_s) (Sect. 3.4.2) are time dependent, showing lagged behaviour behind the initial transient increase and following reduction in net ecosystem exchange during light periods. However, the g_s reduction of about 10% cannot explain the 20–25% NEE reduction as it relaxes to steady-state values. This is consistent with McAusland et al. (2015, 2016), who argued that there should be a fluctuating-light induced decoupling between NEE and g_s at the leaf level.

• Using g_H, g_q, τ_leaf, and τ_air we calculated dynamic changes in the thickness of the transient thermal IBL (Δz) and the thickness of the thermally-active big-leaf layer (l*) (Sect. 3.4.2). The transient thermal IBL thickness above the forest canopy reaches ≈ 20 m at the HF site and ≈ 40–45 m at the LBA site, consistent with the larger LBA temperature time constant for air (Sect. 3.2), larger heat-flux time constant (Sect. 3.3.1), and larger transient duration for the heat flux (Sect. 3.3.2). The thermally-active big-leaf layer thickness, at about 1 mm, however, is about the same for both forests, corresponding to a penetration into the canopy, equivalent to approximately two LAI units. This value is comparable to previous estimates of the outer-canopy thickness over which up to 90% of the incident photosynthetically-active radiative flux is being absorbed (Yoda 1974, Fig. 6; Ellsworth and Reich 1993, Fig. 1; Parker 1995, Fig. 9; Lalic et al. 2013, Figs. 2 and 3).

At both HF (midlatitude deciduous forest) and LBA (tropical rainforest) sites, the cumulative transient net-ecosystem-exchange responses (Fig. 8) are consistent with there being adaptive behaviour, with maxima close to, but slightly below the respective average sunlit duration conditions (Table 4), with the HF maxima better defined than are the LBA maxima (Fig. 9). Even though the HF and LBA ecosystem responses to light change have similar shapes (Figs. 6, 7), they are distinct in terms of both average and flux time constants (Fig. 5; Table 3) and net-ecosystem-exchange and evapotranspiration maxima and magnitudes (Figs. 9, 13). The biggest difference between the two forests is seen in the transient-response times and cumulative transient net-ecosystem-exchange responses in high photosynthetically-active radiative flux values (Fig. 8). The pronounced cumulative transient reduction especially seen for the HF site, but only moderate reduction observed for the LBA site. To reiterate, one possible physiological cause to this in terms of the stomatal control to water stress, perhaps, is due to possible lower water availability at the HF site than at the LBA site, where the trees are known to have very deep roots (Nepstad et al. 1994).

Table 4 Summary of transient forest-response properties for the most common fluctuating-light conditions (in light, S↓_l, and in shadow, S↓_s; first rows in the Table)

The similarity in the sensitivity following light transition between net ecosystem exchange and evapotranspiration and their successful approximations by exponential relaxation curves (Sect. 3.3) represents an observational vindication of the use of the big-leaf hypothesis in describing dynamic outer-canopy behaviour in fluctuating light. This result represents an observational validation of the utility of this parametrization even during dynamic photosynthesis. For this phenomenon, the model is more than an analogy.

Though our realizations cover an entire midday period, we found no 20% net-ecosystem-exchange decrease associated with high radiative flux, the so-called ‘midday rest’ seen on clear days (e.g., Kostychev et al. 1926, 1930; Kostychev and Kardo-Sisoyeva 1930; Kostychev and Berg 1930; Freedman et al. 2001; Rocha et al. 2004). This can be explained by noting that we concentrate on partly-cloudy days, when there are likely to be less thermal- and water-stress–related photoinhibition (Geiger 1928; Maskell 1928a, b; Stalfelt 1935; Friend et al. 1962; Sawada 1978; Powles 1984; Long et al. 1994).

Examining the ensemble-flux change during the fluctuating-light events for real conditional-event distributions using selected cumulus-cloud events with a sharp change in radiative flux (see Sect. 4), we found an absolute increase in net ecosystem exchange over the clear-period steady-state maxima: 20% for the HF site and 10% for the LBA site that occurs near the preferred average light conditions (Fig. 14). For the lower average radiative flux, the net-ecosystem-exchange maxima move to the higher cloudiness values, and this appears to support the diffuse-radiative-flux carbon-uptake enhancement hypothesis (e.g., Gu et al. 2002). However, we note that the corresponding conditional events have higher radiative flux during shorter light periods, and this increase in net ecosystem exchange is actually associated with this high-radiative-flux period and not necessarily with the diffuse-fraction increase. Even though the selected long-shadow-period events ($t_{ - } \ge 200 s$) are not typical for the fair-weather conditions with the average shadow durations less than 100 s (see Kivalov and Fitzjarrald 2018), they cover about 20% of all cumulus-cloud events for both HF and LBA sites, which could make them significant from the plant-adaptation perspective.

We suggest that the twin features of enhanced carbon uptake and increased water-use efficiency on partly-cloudy days can be explained by noting that the uptake enhancement is realized during the sunlit periods, and the resting shadow periods account for the overall increases in the water-use efficiency (see Sect. 4). We believe that the correlation of net ecosystem exchange with diffuse radiative flux, while perhaps robust, does not indicate a causal connection. Rather, it is an emergent property that becomes apparent upon averaging. One benefit of understanding how the whole canopy responds to realistic skies will be in justifying how leaf-scale responses can be scaled-up to whole canopies or ecosystems, what is now done largely by appeal to analogy. Note that the many accumulating tower eddy-flux datasets can be more effectively used to understand of whole-canopy processes by parsing the raw data in novel ways. This different understanding should inform our understanding of whole-canopy photosynthetic mechanisms. In subsequent work, we aim to perform this analysis on data from other ecosystems and to examine simple models of whole-canopy dynamic photosynthesis.

Appendices

Appendix 1: Addressing Ensemble-Flux Convergence Though the Product of the Individual Standard Deviation Representation and Accessing Changes in Correlation Coefficients

When estimating the variance of the product in the convergence Eq. 14 through the individual variances: $\sigma^{2}_{{w^{\prime}X^{\prime}}} \propto \sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}}$, we note that in Eq. 12 the variables should be correlated for flux estimations. Means of the variables are set to zero: $\mu_{{X^{\prime}}} = \mu_{{w^{\prime}}} = 0$, and we assume that they have finite variances: $\sigma^{2}_{{X^{\prime}}}$ and $\sigma^{2}_{{w^{\prime}}}$. Also we accept a general assumption of continuity of processes in the atmosphere, so that the fluctuations should be bounded, within some intervals around a zero mean, with the size of fluctuating eddies within the inertial subrange (Kolmogorov 1941).

For our specific cases, we estimate the degree of dependence of the ensemble variables and their product ensemble variance $\sigma^{2}_{{w^{\prime}X^{\prime}}}$ through the covariances (Goodman 1960),

$$ \begin{aligned} \sigma^{2}_{{w^{\prime}X^{\prime}}} & = \sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}} + \sigma^{2}_{{X^{\prime}}} \left( {\mu_{{w^{\prime}}} } \right)^{2} + \sigma^{2}_{{w^{\prime}}} \left( {\mu_{{X^{\prime}}} } \right)^{2} + {\text{Cov}}\left( {X^{\prime2} ,w^{\prime2} } \right) \\ & \quad - {\text{Cov}}^{2} \left( {X^{\prime},w^{\prime}} \right) - 2{\text{Cov}}\left( {X^{\prime},w^{\prime}} \right)\mu_{{X^{\prime}}} \mu_{{w^{\prime}}} . \\ \end{aligned} $$

(26)

Taking into account that $\mu_{{X^{\prime}}} = \mu_{{w^{\prime}}} = 0$, Eq. 26 can be simplified to

$$ \frac{{\sigma^{2}_{{w^{\prime}X^{\prime}}} }}{{\sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}} }} = 1 + \left[ {\frac{{{\text{Cov}}\left( {X^{\prime2} ,w^{\prime2} } \right)}}{{\sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}} }} - \frac{{{\text{Cov}}^{2} \left( {X^{\prime},w^{\prime}} \right)}}{{\sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}} }}} \right] \Rightarrow \sigma^{2}_{{w^{\prime}X^{\prime}}} \propto \sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}} . $$

(27)

Basing on the mean values (Table 5), we estimate the product ensemble variance through the individual ensemble variances of the variables in Eq. 27 for the HF site (Eq. 28a) and the LBA site (Eq. 28b),

$$ \sigma^{2}_{{w^{\prime}X^{\prime}}} \le 1.02\sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}} , $$

(28a)

$$ \sigma^{2}_{{w^{\prime}X^{\prime}}} \le 1.25\sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}} . $$

(28b)

Then from the central limit theorem applied to the similar and uncorrelated realizations, the ensemble flux is seen as the realization drawn from the following normal distribution,

$$ F_{X} \left( {t,x} \right) = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} w_{i}^{\prime} \left( {t,x} \right)X_{i}^{\prime} \left( {t,x} \right)\mathop \to \limits^{d} {\mathbf{N}}\left( {\mu_{F} \left( {t,x} \right),{\raise0.7ex\hbox{${\alpha \sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}} }$} \!\mathord{\left/ {\vphantom {{\alpha \sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}} } N}}\right.\kern-0pt} \!\lower0.7ex\hbox{$N$}}} \right), $$

(29)

where for $N \to \infty ,\alpha \sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}} /N \to 0$. Hence, noting bounded conditions of fluctuations, the $\sigma^{2}_{{w^{\prime}}} \sigma^{2}_{{X^{\prime}}}$ should be bounded too, and there is ensemble flux convergence $F_{X} \left( {t,x} \right) \to \mu_{F} \left( {t,x} \right)$ to some continuous function almost everywhere.

Table 5 Mean values of the individual terms and the result of Eq. 27

One way to estimate fluxes and transport efficiency is to represent them as the product of individual-scalar standard deviations and correlation coefficients ($r_{Xw}$): $F_{X} = r_{Xw} \sigma_{{X^{\prime}}} \sigma_{{w^{\prime}}}$, where $r_{Xw} = \frac{{F_{X} }}{{\sigma_{{X^{\prime}}} \sigma_{{w^{\prime}}} }}$ can be readily estimated through the obtained ensemble fluxes (Sect. 3.3). Figure 15 shows that the correlation coefficients for the friction velocity ($\left| {r_{uw} } \right|$, purple points and lines) are practically not affected by the shadow-to-light transition. The correlation coefficients for net ecosystem exchange ($\left| {r_{{c_{a} w}} } \right|$, green points and lines), and evapotranspiration ($r_{qw}$, blue points and lines) are seen to increase together and keep similar values after the shadow-to-light transition. The most affected are the correlation coefficients for the sensible heat flux ($r_{Tw}$, red points and lines), exhibiting rapid regime changes. Note that, in the course of the adjustment to the rapid onset of light, $r_{Tw}$ moves toward a value (≈ 0.5) characteristic of a very convective conventional surface layer (e.g., Haugen et al. 1971), while $\left| {r_{uw} } \right|$ remains between 0.3 and 0.4 during the entire transition, values typically seen in the roughness sublayer above a forest canopy (Kaimal and Finnigan 1994, p. 80).

Fig. 15

Correlation coefficients ($r_{Xw}$) obtained from the ensemble-flux estimates: $F_{X} = r_{Xw} \sigma_{{X^{\prime}}} \sigma_{{w^{\prime}}}$ near the shadow-to-light transition (t = 0), points—measurements, lines–fits; [900,1100] W m⁻² clear-period S↓ bin; a the HF site, and b the LBA site

Figure 15 shows that the turbulent heat transfer due to buoyancy (red, in light following transition, time > 0) is twice as effective as transport by forced convection (red, in shadow before transition, time < 0). As the shadow-to-light transition occurs, there is the rapid change (within the first 150 s following transition) in surface-layer stability toward more convective conditions. These results clearly show the evolution of the more convective surface layer, ending with correlation coefficients similar to those normally found in an equilibrium surface layer.

Appendix 2: Conditions for Ensemble-Flux Convergence

Following the “if” Lyapunov (1954) condition, the sequence of independent random variables X_i, each with finite expected value $\mu_{i}$ and variance $\sigma_{i}^{2}$, converges to the normal distribution if for some δ > 0 the following is satisfied

$$ {\text{Lyapunov}}: \mathop {\lim }\limits_{N \to \infty } \frac{1}{{\left( {\mathop \sum \nolimits_{1}^{N} \sigma_{i}^{2} } \right)^{{{\raise0.7ex\hbox{${\left( {2 + \delta } \right)}$} \!\mathord{\left/ {\vphantom {{\left( {2 + \delta } \right)} 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}}}} }}\mathop \sum \limits_{1}^{N} \left( {X_{i} - \mu_{i} } \right)^{{\left( {2 + \delta } \right)}} \to 0. $$

(30)

Taking $X_{i}^{\prime} = X_{i} - \mu_{i}$ as a fluctuating part of the time series with regard to its own mean μ_i, we see that for this convergence to occur the average of moments higher than the second (with some δ > 0 in Eq. 30) of the individual distributions should be smaller than the average of their second moments. This will be trivial for the third moment of a symmetrical distribution. This allows us to compare similarity and convergence for the normal distribution to the symmetry (the third moment: δ = 1) of the realizations around their ensemble means. Evaluation for the third moment δ = 1 yields

$$ \frac{1}{{\left( {\mathop \sum \nolimits_{1}^{N} \sigma_{i}^{2} } \right)^{3/2} }}\mathop \sum \limits_{1}^{N} \left( {X_{i} - \mu_{i} } \right)^{3} = \frac{1}{\sqrt N }\frac{1}{{\left( {\frac{1}{N}\mathop \sum \nolimits_{1}^{N} \sigma_{i}^{2} } \right)^{3/2} }}\left[ {\frac{1}{N}\mathop \sum \limits_{1}^{N} \left( {X_{i} - \mu_{i} } \right)^{3} } \right]. $$

(31)

To ensure convergence to zero (Eq. 30), it follows that the following expression holds

$$ \frac{1}{N}\mathop \sum \limits_{1}^{N} \left( {X_{i} - \mu_{i} } \right)^{3} \propto N^{{{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}} - \alpha }} \sigma_{ens}^{3} ; \quad \forall \alpha > 0, $$

(32)

in which $\sigma_{ens} = \sqrt {\frac{1}{N}\mathop \sum \limits_{1}^{N} \sigma_{i}^{2} }$. Now we can evaluate the condition of symmetry of X_i realizations around the ensemble mean: μ_a. Modifying Eq. 32 with μ_a,

$$ \begin{aligned} \frac{1}{N}\mathop \sum \limits_{1}^{N} \left( {X_{i} - \mu_{i} } \right)^{3} & = \frac{1}{N}\mathop \sum \limits_{1}^{N} \left( {\left( {X_{i} - \mu_{a} } \right) - \left( {\mu_{i} - \mu_{a} } \right)} \right)^{3} = \frac{1}{N}\mathop \sum \limits_{1}^{N} \left( {X_{i} - \mu_{a} } \right)^{3} \\ & \quad - \frac{3}{N}\mathop \sum \limits_{1}^{N} X_{i}^{'2} \mu_{i}^{\prime} + \frac{3}{N}\mathop \sum \limits_{1}^{N} \mu_{i}^{'2} X_{i}^{\prime} - \frac{1}{N}\mathop \sum \limits_{1}^{N} \left( {\mu_{i} - \mu_{a} } \right)^{3} = \mu_{X\_ens}^{3} - \mu_{\mu \_ens}^{3} \\ & \quad - \frac{3}{N}\mathop \sum \limits_{1}^{N} X_{i}^{\prime} \mu_{i}^{\prime} \left( {X_{i}^{\prime} - \mu_{i}^{\prime} } \right) \approx \mu_{X\_ens}^{3} - \mu_{\mu \_ens}^{3} . \\ \end{aligned} $$

(33)

From Eqs. 32 and 33, we obtain the following expression for the third moment

$$ \mu_{X\_ens}^{3} - \mu_{\mu \_ens}^{3} \propto N^{{{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}} - \alpha }} \sigma_{ens}^{3} . $$

(34)

Another version of Eq. 32 is for the case with offset removal from the individual realizations: $X_{i}^{off} = X_{i} - \left( {\mu_{i} - \mu_{a} } \right)$. Modifying this with μ_a and taking the offset removal, we directly obtain the following expression for the third moment,

$$ \begin{aligned} \frac{1}{N}\mathop \sum \limits_{1}^{N} \left( {X_{i} - \mu_{i} } \right)^{3} & = \frac{1}{N}\mathop \sum \limits_{1}^{N} \left( {\left( {X_{i} - \left( {\mu_{i} - \mu_{a} } \right)} \right) - \mu_{a} } \right)^{3} = \frac{1}{N}\mathop \sum \limits_{1}^{N} \left( {X_{i}^{off} - \mu_{a} } \right)^{3} \\ & = \mu_{{X^{off} \_ens}}^{3} \propto N^{{{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 2}}\right.\kern-0pt} \!\lower0.7ex\hbox{$2$}} - \alpha }} \sigma_{ens}^{3} . \\ \end{aligned} $$

(35)

Appendix 3: Ensemble-Flux Additive Feature, Its Ergodicity Relative to the Time-Averaged Flux, and Its Invariance to Offset Removal

We define the offsets of realizations through the difference between means (time integrals) of the realization and the ensemble: $X_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right) = \frac{1}{{t_{2} - t_{1} }}\mathop \int \limits_{{t_{1} }}^{{t_{2} }} X_{i} \left( {t,x} \right)dt - \frac{1}{{t_{2} - t_{1} }}\mathop \int \limits_{{t_{1} }}^{{t_{2} }} X\left( {t,x} \right)dt$, where $t \in \left[ {t_{1} ,t_{2} } \right]$ is the base for the offset removal and x is any space-point in the domain.

We derive the new set of the realizations with offset removal,

$$ \begin{aligned} X_{i}^{{\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right) & = X_{i} \left( {t,x} \right) - X_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right) = X_{i} \left( {t,x} \right) - \frac{1}{{t_{2} - t_{1} }}\mathop \int \limits_{{t_{1} }}^{{t_{2} }} X_{i} \left( {t,x} \right)dt \\ & \quad + \frac{1}{{t_{2} - t_{1} }}\mathop \int \limits_{{t_{1} }}^{{t_{2} }} \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i} \left( {t,x} \right)dt. \\ \end{aligned} $$

(36)

The new set of the realizations has the same ensemble as the original set of the realizations,

$$ \begin{aligned} {\langle}X^{{\left[ {t_{1} ,t_{2} } \right]}}{\rangle} \left( {t,x} \right) & = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i} \left( {t,x} \right) - \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i} \left( {t,x} \right) \\ & \quad - \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \frac{1}{{t_{2} - t_{1} }}\mathop \int \limits_{{t_{1} }}^{{t_{2} }} X_{i} \left( {t,x} \right)dt + \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \frac{1}{{t_{2} - t_{1} }}\mathop \int \limits_{{t_{1} }}^{{t_{2} }} \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i} \left( {t,x} \right)dt \\ & = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i} \left( {t,x} \right) + \frac{1}{{t_{2} - t_{1} }}\mathop \int \limits_{{t_{1} }}^{{t_{2} }} \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i} \left( {t,x} \right)dt - \frac{1}{{t_{2} - t_{1} }}\frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \mathop \int \limits_{{t_{1} }}^{{t_{2} }} X_{i} \left( {t,x} \right)dt \\ & = \left\langle X \right\rangle (t,x). \\ \end{aligned} $$

(37)

Applying the Reynolds decomposition, we calculate the fluctuations from the new set of realizations,

$$ \begin{aligned} X_{i}^{{\left[ {t_{1} ,t_{2} } \right]\prime}} \left( {t,x} \right) & = X_{i}^{{\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right) - \left\langle {X^{{\left[ {t_{1} ,t_{2} } \right]}} } \right\rangle \left( {t,x} \right) = \left( {X_{i} \left( {t,x} \right) - X_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right)} \right) - \left\langle X \right\rangle \left( {t,x} \right) \\ & = X_{i} \left( {t,x} \right) - \frac{1}{{t_{2} - t_{1} }}\mathop \int \limits_{{t_{1} }}^{{t_{2} }} X_{i} \left( {t,x} \right)dt + \frac{1}{{t_{2} - t_{1} }}\mathop \int \limits_{{t_{1} }}^{{t_{2} }} \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i} \left( {t,x} \right)dt - \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i} \left( {t,x} \right) \\ & = \left[ {X_{i} \left( {t,x} \right) - \frac{1}{{t_{2} - t_{1} }}\mathop \int \limits_{{t_{1} }}^{{t_{2} }} X_{i} \left( {t,x} \right)dt} \right]\\ &- \left[ {\frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \left( {X_{i} \left( {t,x} \right) - \frac{1}{{t_{2} - t_{1} }}\mathop \int \limits_{{t_{1} }}^{{t_{2} }} \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i} \left( {t,x} \right)dt} \right)} \right] \\ & \approx \left[ {X_{i} \left( {t,x} \right) - \overline{{X_{i} }}^{T} } \right] - \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \left( {X_{i} \left( {t,x} \right) - \overline{{X_{i} }}^{T} } \right) = X_{i}^{{T^{\prime}}} \left( {t,x} \right) - \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i}^{{T^{\prime}}} \left( {t,x} \right). \\ \end{aligned} $$

(38)

For the fluxes we find,

$$ \begin{aligned} {\langle}w^{{\left[ {t_{1} ,t_{2} } \right]^\prime}} X^{{\left[ {t_{1} ,t_{2} } \right]^\prime}} {\rangle}& = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \left( {w_{i}^{{T^{\prime}}} \left( {t,x} \right) - \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} w_{i}^{{T^{\prime}}} \left( {t,x} \right)} \right)\left( {X_{i}^{{T^{\prime}}} \left( {t,x} \right) - \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} x_{i}^{{T^{\prime}}} \left( {t,x} \right)} \right) \\ & \approx {\langle}w^{T'} X^{T'}{\rangle} - \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \left( {\frac{1}{N}\mathop \sum \limits_{i = 1}^{N} w_{i}^{{T^{\prime}}} \left( {t,x} \right)} \right)\left( {X_{i}^{{T^{\prime}}} \left( {t,x} \right)} \right) \\ & \quad - \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \left( {w_{i}^{{T^{\prime}}} \left( {t,x} \right)} \right)\left( {\frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i}^{{T^{\prime}}} \left( {t,x} \right)} \right)\\ &+ \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \left( {\frac{1}{N}\mathop \sum \limits_{i = 1}^{N} w_{i}^{{T^{\prime}}} \left( {t,x} \right)} \right)\left( {\frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i}^{{T^{\prime}}} \left( {t,x} \right)} \right) \\ & = \left\langle {w^{T'} X^{T'} } \right\rangle - \left\langle {w^{T'} } \right\rangle \left\langle {X^{T'} } \right\rangle . \\ \end{aligned} $$

(39)

Comparing the ensemble fluxes without $\left\langle {w^{\prime } X^{\prime } } \right\rangle$ and with $\left\langle {w^{{\left[ {t_{1} ,t_{2} } \right]^{\prime } }} X^{{\left[ {t_{1} ,t_{2} } \right]^{\prime } }} } \right\rangle$ the offset removed,

$$ \begin{aligned} {\langle}w^{{\left[ {t_{1} ,t_{2} } \right]^\prime}} X^{{\left[ {t_{1} ,t_{2} } \right]^\prime}}{\rangle} & = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} \left( {w_{i}^{\prime} \left( {t,x} \right) - w_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right)} \right)\left( {X_{i}^{\prime} \left( {t,x} \right) - X_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right)} \right) \\ & = \left\langle {w^{\prime}X^{\prime}} \right\rangle - \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} w_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right)X_{i}^{\prime} \left( {t,x} \right) \\ & \quad - \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} w_{i}^{\prime} \left( {t,x} \right)X_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right) + \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} w_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right)X_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right) \\ & = \left\langle {w^{\prime}X^{\prime}} \right\rangle - \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} w_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right)X_{i}^{\prime} \left( {t,x} \right)\\ &- \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} w^{{\left[ {t_{1} ,t_{2} } \right]\prime}} \left( {t,x} \right)X_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right), \\ \end{aligned} $$

(40)

we see that there are two terms distinguishing our initial ensemble flux $\left\langle {w^{\prime } X^{\prime } } \right\rangle$ from the ensemble flux with the offsets removed $\left\langle {w^{{\left[ {t_{1} ,t_{2} } \right]\prime}} X^{{\left[ {t_{1} ,t_{2} } \right]\prime}} } \right\rangle$.

The first term is $\frac{1}{N}\mathop \sum \nolimits_{i = 1}^{N} w_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right)X_{i}^{\prime} \left( {t,x} \right)$ and represents the correlation between the initial fluctuation X′ and the offset size $w_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right)$. Because w always fluctuates around zero, $w_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right) = \frac{1}{{t_{2} - t_{1} }}\mathop \int \nolimits_{{t_{1} }}^{{t_{2} }} w_{i} \left( {t,x} \right)dt - w\left( {t,x} \right) \approx 0 - 0 = 0$, which effectively cancels the first term.

The second term is $\frac{1}{N}\mathop \sum \nolimits_{i = 1}^{N} w^{{\left[ {t_{1} ,t_{2} } \right]\prime}} \left( {t,x} \right)X_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right)$ and represents the correlation between the fluctuation $w^{{\left[ {t_{1} ,t_{2} } \right]\prime}} \left( {t,x} \right)$ after the offsets are removed and the offset size $X_{i}^{{off\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right)$. However, after the offsets are removed, the $w^{{\left[ {t_{1} ,t_{2} } \right]}} \left( {t,x} \right)$ terms are uniformly placed around the ensemble mean, implying that any offset-related correlations in $w^{{\left[ {t_{1} ,t_{2} } \right]\prime}} \left( {t,x} \right)$ are removed, effectively cancelling the second term.

Appendix 4: Dynamic Air and Leaf (Canopy) Energy Balances

The thermal balance equation in a layer of air is

$$ \frac{{dT_{a} }}{dt} = \frac{{\partial T_{a} }}{\partial t} + u_{i} \frac{{\partial T_{a} }}{{\partial x_{i} }} = \vartheta \frac{{\partial^{2} T_{a} }}{{\partial x_{i}^{2} }} - \frac{1}{{\rho c_{p} }}\left( {\frac{\partial Q}{{\partial x_{i} }}} \right), $$

(41)

$\vartheta = \frac{k}{{\rho c_{p} }}$ is the thermal diffusivity, $Q = R_{net} + H + \lambda E$ is the sum of the net radiative flux and sensible and latent heat fluxes (positive upward). If we neglect advection $u_{i} \frac{{\partial T_{a} }}{{\partial x_{i} }} = 0$ and phase change $\frac{\partial \lambda E}{\partial z} = 0$, and note that the thermal diffusivity in air is much smaller than is the convective and radiative heat transfer $k\frac{{\partial^{2} T_{a} }}{{\partial x_{i}^{2} }} \ll \frac{{\partial R_{net}^{air} }}{\partial z} + \frac{\partial H}{\partial z}$, we approximate the budget as,

$$ \frac{{dT_{a} }}{dt} = - \frac{1}{{\rho c_{p} }}\left( {\frac{{\partial R_{net}^{air} }}{\partial z} + \frac{\partial H}{\partial z}} \right). $$

(42)

Assuming a constant-flux layer, we write the budget in finite difference form, with Δz representing the thickness of the transient thermal internal boundary layer,

$$ \rho c_{p} \Delta z\frac{{dT_{a} }}{dt} \approx \left( {R_{net + } - R_{net - } } \right) + \left( {H_{ + } - H_{ - } } \right). $$

(43)

Absorption of the shortwave radiative flux across the Δz layer of air is approximated through the Bouguer law of extinction,

$$ R_{net + } - R_{net - } = R_{0} \left( {{\text{e}}^{{ - \beta z_{ + } }} - {\text{e}}^{{ - \beta z_{ - } }} } \right) \approx R_{0} \beta \Delta z, $$

(44)

while the turbulent and radiative-flux divergence is written,

$$ H_{ + } - H_{ - } = \frac{{\rho c_{p} }}{{r_{HR} }}\left( {T_{a + } - T_{a} } \right) - \frac{{\rho c_{p} }}{{r_{HR} }}\left( {T_{a} - T_{a - } } \right) = 2\frac{{\rho c_{p} }}{{r_{HR} }}\left( {\frac{{T_{a + } + T_{a - } }}{2} - T_{a} } \right), $$

(45)

where $T_{e} = \frac{{T_{a + } + T_{a - } }}{2}$ represents an average environmental temperature.With the symmetry of T_a+ and T_a− around T_e, this satisfies constant-flux-layer conditions $\frac{{dT_{a} }}{dt} = \vartheta \frac{{\partial^{2} T_{a} }}{{\partial z^{2} }} = 0,$ making $\frac{{\partial T_{a} }}{\partial z}$ constant. Combining Eqs. 43, 44, 45,

$$ \rho c_{p} \Delta z\frac{{dT_{a} }}{dt} \approx R_{0} \beta \Delta z + 2\frac{{\rho c_{p} }}{{r_{HR} }}\left( {T_{e} - T_{a} } \right). $$

(46)

Taking clear-air conditions $\beta \ll 1$, and the air time constant $\tau_{air} = r_{HR} \Delta z/2$, we find the final air-layer temperature approach-to-equilibrium equation,

$$ \tau_{air} \frac{{dT_{a} }}{dt} = \left( {T_{e} - T_{a} } \right), $$

(47a)

$$ T_{a} \left( t \right) = T_{e} - \left( {T_{e} - T_{0}^{air} } \right){ \exp }\left\{ { - {\raise0.7ex\hbox{$t$} \!\mathord{\left/ {\vphantom {t {\tau_{air} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${\tau_{air} }$}}} \right\}. $$

(47b)

If we take the air layer centred around the site of the measurements, T_a can be replaced by T_air, and Eqs. 47a, b describe the change in temperature at the measurement site as well.

For the big-leaf thermal balance, we follow Jones (1992) for absolute flux values,

$$ \rho^{*} c_{p}^{*} l^{*} \frac{{dT_{l} }}{dt} = R_{net} - C - \lambda E. $$

(48)

In this case l^* is interpreted as the thickness of the thermally-active big-leaf layer, and $\rho^{*} c_{p}^{*}$ is the leaf heat capacity in terms of leaf water content. For the single-leaf case, $\rho^{*} = 700 {\text{kg m}}^{ - 3}$ and $c_{p}^{*} = 3800 {\text{J kg}}^{ - 1} {\text{K}}^{ - 1}$, giving $\rho^{*} c_{p}^{*} = 2.7 {\text{MJ m}}^{ - 3} {\text{K}}^{ - 1}$ (Jones 1992).

The equilibrium leaf temperature ($T_{eq}^{leaf}$) can be described through the Penman–Monteith equation (Monteith and Unsworth 1990; Jones 1992),

$$ \frac{{\rho c_{p} }}{{r_{HR} }}\left( {T_{eq}^{leaf} - T_{air} } \right) = C = \frac{{\gamma \frac{{\left( {r_{w} + r_{s} } \right)}}{{r_{HR} }}R_{neti}^{leaf} - \rho c_{p} \frac{\delta e}{{r_{HR} }}}}{{\left( {\gamma \frac{{\left( {r_{w} + r_{s} } \right)}}{{r_{HR} }} + s} \right)}}, $$

(49)

where $R_{net}^{leaf} = R_{neti}^{leaf} - \left( {\frac{{4\sigma T_{a}^{3} }}{{\rho c_{p} }}} \right)\rho c_{p} \left( {T_{l} - T_{a} } \right) = R_{neti}^{leaf} - \frac{{\rho c_{p} }}{{r_{R} }}\left( {T_{l} - T_{a} } \right)$ and $\frac{1}{{r_{HR} }} = \frac{1}{{r_{R} }} + \frac{1}{{r_{H} }}$.

Jones (1992) finds the leaf-temperature approach-to-equilibrium equation,

$$ \tau_{leaf} \frac{{dT_{l} }}{dt} = \left( {T_{eq}^{leaf} - T_{l} } \right), $$

(50)

in which $\tau_{leaf} = \frac{{\rho^{*} c_{p}^{*} l^{*} }}{{\rho c_{p} \left\{ {\frac{1}{{r_{HR} }} + \frac{s}{{\gamma \left( {r_{w} + r_{s} } \right)}}} \right\}}}$ is the leaf-temperature time constant. Note that to estimate l^* we still require all the resistances describing the sensible and latent heat fluxes.

The time constant for the radiative flux (τ_light) is 3–4 s (LBA site, see Fig. 1b) and 6–7 s (HF site, see Fig. 1a) so we can assume that the radiative flux reaches its equilibrium value almost instantaneously when compared to the leaf or air temperature or the stomatal conductance and $\tau_{light} \ll \tau_{{leaf_{T} }} \ll \tau_{stomata} < \tau_{air\_T}$. Assuming small stomatal conductance at the end of the shadow and beginning of the light period, $r_{s} \to \infty$ in Eq. 49, we obtain the first estimation for the equilibrium leaf-temperature change,

$$ T_{eq}^{leaf} = \frac{{r_{HR} }}{{\rho c_{p} }}R_{neti} + T_{air} . $$

(51)

This estimate holds for the small initial time period because $\tau_{leaf\_T} \ll \tau_{stomata}$, and we can extract adequate leaf-temperature time constants from it. The justification for assuming a small stomatal conductance arises from noting that the ensemble evapotranspiration estimates during the shadow periods are about four times smaller than are those during the light periods (see Figs. 6, 7).

Because by definition R_netti does not depend on the leaf temperature, Eq. 51 yields a linear relationship of equilibrium leaf temperature with air temperature. Using the approach-to- equilibrium air temperature Eq. 47b and combining Eqs. 50 and 51,

$$ \tau_{leaf} \frac{{dT_{l} }}{dt} = \left( {\frac{{r_{HR} }}{{\rho c_{p} }}R_{neti} + T_{air} \left( t \right)} \right) - T_{l} , $$

(52a)

whose exact solution is

$$ \begin{aligned} T_{l} \left( t \right) & = \frac{{r_{HR} }}{{\rho c_{p} }}R_{neti} + T_{eq}^{air} - \left( {T_{eq}^{air} - T_{0}^{air} } \right)\frac{{\tau_{air} }}{{\tau_{air} - \tau_{leaf} }}\exp \left\{ { - {\raise0.7ex\hbox{$t$} \!\mathord{\left/ {\vphantom {t {\tau_{air} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${\tau_{air} }$}}} \right\} \\ & \quad - \left( {\frac{{r_{HR} }}{{\rho c_{p} }}R_{neti} + T_{eq}^{air} - \left( {T_{eq}^{air} - T_{0}^{air} } \right)\frac{{\tau_{air} }}{{\tau_{air} - \tau_{leaf} }} - T_{l0} } \right)\exp \left\{ { - {\raise0.7ex\hbox{$t$} \!\mathord{\left/ {\vphantom {t {\tau_{leaf} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${\tau_{leaf} }$}}} \right\}. \\ \end{aligned} $$

(52b)

Following the Penman–Monteith-equation approach, the same Eq. 52b can be applied to the whole-canopy scale by replacing the stomatal conductance (g_s) with the canopy conductance (g_c) or canopy resistance ($r_{c} = 1/g_{c}$) through $g_{s} \approx g_{c} /LAI$. Then τ_leaf is replaced by $\tau_{canopy} = \frac{{\rho^{*} c_{p}^{*} l^{*} }}{{\rho c_{p} \left\{ {\frac{1}{{r_{HR} }} + \frac{s}{{\gamma \left( {r_{w} + r_{c} } \right)}}} \right\}}}$.