The importance of radiation for semi-empirical water-use efficiency models

Water-use efficiency (WUE) is a fundamental property for the coupling of carbon and water cycles in plants and ecosystems. Existing model formulations predicting this variable differ in the type of response of WUE to the atmospheric vapor pressure deficit of water (VPD). We tested a representative WUE model on ecosystem scale at 110 eddy-covariance sites of the FLUXNET initiative by predicting evapotranspiration (ET) based on gross primary productivity (GPP) and VPD. We found 10 that introducing an intercept term in the formulation increases model performance considerably, indicating that an additional factor needs to be considered. We demonstrate that this intercept term varies seasonally and we subsequently associate it with radiation. Replacing the constant intercept term with a linear function of global radiation was found to further improve model predictions of ET. Our new semi-empirical ecosystem WUE formulation indicates that, averaged over all sites, this radiation term accounts for up to half (40–49 %) of transpiration. These empirical findings challenge the current 15 understanding of water-use efficiency on ecosystem-scale.

obtain data of the gross-primary productivity (GPP) and the evapotranspiration (ET). On ecosystem scale, water-use efficiency is then defined as WUE = GPP ET, while the intrinsic water-use-efficiency is accordingly formulated as WUE ' = GPP ) , where G s denotes the surface conductance of the ecosystem.
Analytical models that predict WUE on the leaf-scale (Katul et al., 2010;Medlyn et al., 2011) were derived from theoretical 5 optimality considerations. Corresponding variants were evaluated with ecosystem-scale flux measurements gathered by the FLUXNET in empirical studies (Beer et al., 2009;Zhou et al., 2014;Zhou et al., 2015). The central difference between the existing models is the response of g s or G s to VPD, resulting in different dependencies of WUE on VPD. The concept by Beer et al. (2009) assumed an inverse dependency between WUE and VPD due to an increasing diffusion of water vapor with higher values of VPD. On the other hand, physiological models (e.g. Katul et al., 2010), imply stomatal contraction 10 with increasing VPD, as plants aim to minimize water loss. This was found to be an important factor on ecosystem scale as shown by Zhou et al. (2014) for half-hourly and by Zhou et al. (2015) for daily observations. In both studies, the response of WUE could be approximated to be proportional to VPD -0.5 . Despite their discrepancy, both the models of Beer et al. and Zhou et al. imply that, at ecosystem-scale, WUE is strictly an inverse function of VPD or VPD 0.5 , respectively.

15
By contrast, physiological WUE models that were derived on plant scale include an additional conductance term g 0 that is invariable to changing levels of photosynthesis (Ball et al., 1987;Medlyn et al., 2011). Any transpiration resulting from this part of stomatal conductance should be expected to be proportional to the gradient of the partial pressure of water vapor, quantified by the atmospheric VPD observations. This presents a process that current ecosystem-scale WUE models do not accommodate. 20 Incoming solar radiation is a driving factor for both photosynthesis and transpiration, yet in the existing WUE model formulations, radiation is not explicitly considered. The models implicitly assume that, at ecosystem-scale, GPP and ET respond equally to changes in radiation and that, therefore, the ratio of both is constant with regard to this factor. However, photosynthesis saturates at high radiation levels (Farquhar et al., 1980), even under well-watered conditions. Models of 25 Potential Evapotranspiration (PET), by contrast, do not include a similar limiting behavior (Leuning et al., 2008).
Mechanistically, the process of equilibrium transpiration (Jarvis and McNaughton, 1986) implies that sizeable transpiration can occur even when the leaf is fully decoupled from the atmosphere, i.e. when VPD is very low. This is a second process that current ecosystem WUE models cannot accommodate.

30
In this study, we address these unresolved inconsistencies regarding the importance of additional model terms for predictions of ecosystem-scale transpiration. We do this empirically on ecosystem scale by optimizing and assessing different WUE models with FLUXNET observations from 110 globally distributed towers. In our approach, ET is selected as target variable, while the different WUE models utilize GPP as one of multiple explanatory variables. The substantial degree of Biogeosciences Discuss., doi: 10.5194/bg-2016-524, 2017 Manuscript under review for journal Biogeosciences Published: 5 January 2017 c Author(s) 2017. CC-BY 3.0 License. correlation between GPP and ET is thus harnessed for the predictions of ET. In a first step, we identify existing biases in ecosystem-scale WUE models. In the next step, these biases are tested for their dependency on VPD and radiation. Lastly, we infer a tentative partitioning of transpiration according to its association with radiation and discuss the substantial magnitude of this metric.

Data
The daily day-time integrals of gross primary productivity (GPP) and evapotranspiration (ET) were taken from the La Thuile FLUXNET (open and fair use data policy sites) collection (Baldocchi et al., 2001). The aggregation to day-time values was 10 based on values of potential radiation larger than 10 W m -2 . Additionally, we used global radiation (Rg) and the day-time vapor pressure deficit (VPD) measured at the same eddy-covarance (EC) sites. For some analyses we further used long-term data of precipitation (P) and potential evapotranspiration (PET) based on a bias-corrected ERA-Interim reanalysis. 1 The EC data were processed according to the standard methods (Papale et al., 2006;Reichstein et al., 2005) to assure 15 consistent quality of the observations. Eddy-covariance GPP results were based on the flux partitioning method of Reichstein (2005). We used only data with GPP > 0.1 gC d -1 m -2 , ET > 0.05 mm d -1 and VPD > 0.001 kPa) to reduce the relatively large impact of random measurement errors under low flux conditions. Following the procedure of Beer et al. (2009) we further used only data after three consecutive rain free days. This reduces contributions by evaporation to the measured evapotranspiration as physical evaporation typically declines rapidly after rain events due to the depletion of water stored in 20 the topmost soil layer (Wythers et al., 1999). That assumption similarly applies to precipitation that is intercepted on leafsurfaces and other plant parts in the canopy. As Miralles et al. (2010) summarize, the interception storage for forest ecosystems reported in different studies amounted to a mean of 1.2 mm (± 0.4 mm). With the mean interception evaporation rates reported as 0.3 mm h -1 (± 0.1 mm h -1 ), this storage can be concluded to be typically depleted within the first days after a precipitation event. Therefore, the measured evapotranspiration after three consecutive rain free days is expected to 25 approximate transpiration, and we additionally verified that our results are robust when considering longer rain free periods (see below). For the main analysis, we included sites with at least 25 data-points fulfilling the requirements noted above. A list of these 110 sites used for the parameter estimation can be found in the Supplement (S2).
The presented analyses presume that the observed evapotranspiration is dominated by transpiration after three consecutive rain-free days. To test the robustness of our findings against this assumption, we varied the number of consecutive rain free days from 0 to 14. For each considered step, the data were filtered accordingly and parameters of the WUE models were estimated for each site. We then calculated the mean of the parameter estimates over all sites and the associated uncertainty (95% confidence interval) of the mean via bootstrapping. We excluded some humid sites from this sensitivity analysis for 5 which too few data points were available when filtering for longer rain free periods. This procedure can thus ensure that When analyzing the collinearity of parameters, we tested to what degree our results were affected by water limitation. The criterion used to classify sites into arid and non-arid was a cumulative difference between PET and P of more than 300 mm annually. 10

Concepts & models
For our analysis, we started out with the WUE model of Zhou et al. (2015), which we converted for an inversion against ET Hence, any significant intercept would indicate transpiration that cannot be explained by GPP • VPD =.> . The relative contribution of ET res to the mean predicted flux, C ETres , was calculated as: 25 where the denominator contains the mean predicted daily ET of the model.
To further test whether ET res relates to atmospheric variables, we postulated three different alternative models. The residual transpiration could be driven by an additional VPD term that is independent from photosynthesis, as some stomatal conductance models (Medlyn et al., 2011) We also considered the possibility that ET res is related to global radiation (Rg):

ET =
/00 •203 4.6 789: To test for possible interactions between the two additional variables, we considered a third alternative where ET res was modelled by included both Rg and VPD as independent factors: ET = /00 •203 4.6 789: (6) 15 In the following, we refer to the reference WUE definition (Eq. 1) as "Zhou". We abbreviate models with additional covariates by omitting the reference to the GPP-VPD-term ( /00 • 203 4.6 789: ) of Zhou, which is always used unless denoted otherwise. The model with an additional VPD term, for example, is thus designated "+VPD". 20

Parameter estimation & statistics
In the following, we refer to Eq. 1, 2, 4, 5, 6 as models, as we optimize their fit to the EC data by estimating free parameters.
The estimation was conducted with the Levenberg-Marquardt technique, minimizing the sum of squares of the model 25 residuals. The residuals were calculated as the difference between observed and predicted ET. We used the nlsLM package in R (Elzhov et al., 2015). The uncertainties and correlations of the parameters were calculated with the variance-covariance matrix returned by the fitting function (Omlin & Reichert, 1999).
where `aB are observations and b@c are predictions by a model.
A MEF of 1 implies a perfect fit of the model to the data. A MEF below 0 implies that the mean of the observations 5 outperforms the fit of the model. All MEFs were calculated in a leave-one-out cross-validation to account for the problem of over-fitting. Thus, the cross-validated MEFs can be used to compare models with differing numbers of free parameters.

2.4
Partitioning of linear models

10
To assess the contribution of driving variables to the predicted fluxes, we performed an analysis of attribution to the individual model terms. Consider a simple multiple linear regression model,

15
where is the dependent variable, e and g are independent variables and a 1 and a 2 denote the model parameters. Due to the additive character of the model, the contribution of one variable (e.g. e ) to the total flux is given by its product with the slope ( • e ). In our analysis, we tested the contribution of the linear radiation term (Eq. 5) to the total modelled evapotranspiration. Thus we defined the fraction of evapotranspiration that was attributed to the radiation term as: We considered two variants to estimate this metric: parallel and hierarchical. In the first case, both parameters, uWUE and r were estimated in a standard parameter estimation, i.e. concurrently. In the second case, uWUE was first estimated in the +ET res model. We then defined a term M as: 25 where denotes the residuals and ET res denotes the intercept parameter of the +ET res model. The parameter r was then estimated in a linear regression of the form: 30 Biogeosciences Discuss., doi: 10.5194/bg-2016-524, 2017 Manuscript under review for journal Biogeosciences Published: 5 January 2017 c Author(s) 2017. CC-BY 3.0 License.
where c denotes a constant intercept. By giving precedence to the uWUE parameter in this approach, we expect to get an estimate of a reasonable lower bound for ET frac . All parameters were constrained to positive values.

5
To further assess uncertainties due to problems of parameter identifiability among uWUE and r in the parallel variant, we sampled 200 parameter vectors from the posterior parameter uncertainty distribution for each site (See Supplement S1 for a full description of the applied procedure).

Results
In our analysis we tested different water-use efficiency (WUE) models that predicted evapotranspiration (ET) using the product of gross-primary productivity (GPP) and the water vapor pressure deficit, GPP • VPD =.> as predictor variable. When plotting ET as a function of STO, we observed significant intercepts, e.g. for the Mediterranean FLUXNET site IT-BCi (Fig.  15 1a). In these cases, significant ET was observed when the driving force of the established models, the GPP-VPD-term, was small or zero. When we explicitly included this term in the model (+ETres), the cross-validated MEF increased notably (Fig.   1b). As Table 1 shows, the +ETres variant outperformed the Zhou model at 91% of the sites. The respective mean difference in MEF between the two variants was 0.08 (Table 2)

5
The site specifically estimated intercept values range from -0.43 to 2.88 mm (1.25 mm when excluding the highest value) with a mean of 0.34mm (Fig. 2a). The mean relative intercept, C ETres , was 0.25 when averaged over all sites, implying that a quarter of transpiration is was not attributed to GPP • VPD =.> in this model formulation (Fig. 2b). The importance of the intercept for the prediction of daily ET and its diverging values raise the question whether the intercept compensates for the absence of a physical or biological process in the model or whether there is a problem with the observations. 10 Our first hypothesis was that the ET res intercept was due to remaining contributions of soil and interception evaporation to measured ET after three consecutive rain-free days. To test this, we estimated the parameter ET res for periods of successively longer consecutive rain free days. If our hypothesis was right, we would expect a trend of declining ET res with increasing consecutive rain free days. However, no reduction of ET res beyond the exclusion of the three days after precipitation 15 proposed by Beer et al. (dotted line) could be observed (Fig. 3a). We therefore concluded that potential contributions of soil and interception evaporation to ET cannot explain the existence of the intercept term in the WUE model.
We then hypothesized that a missing process representation in the model would be discernible as a temporal pattern in the ET res estimates. For that, we estimated the intercept for each month and site separately. The monthly means of the intercept 20 for all sites varied in a clear seasonal pattern (Fig. 3b).

5
The seasonality of the parameter suggests a relationship to meteorological variables such as VPD or radiation that vary seasonally too. A relationship with VPD could represent the g 0 term of canopy conductance models (Ball et al., 1987). It was therefore introduced in an additional linear term (Eq. 4). A relationship with global radiation (Rg) (Eq. 5) could represent "equilibrium transpiration" (Jarvis and McNaughton, 1986), where the energy surplus of incoming radiation forces a 10 transpirational flux independent from the vapor pressure gradient. We tested modeling the ETres intercept by including both variables separately (Eq. 4, 5) and jointly (Eq. 6).
We found only a small performance increase of the VPD-variant with regard to the Zhou model with intercept term ("+ETres", Fig. 4a). In fact, the MEF diff suggested that the ET res -model was still superior to the VPD-variant at 42% of the 15 sites (Table 1). By comparison, the Rg-model showed a substantial increase in MEF; compared to the ET res -model. The increased model complexity (two free parameters) was justified at 81% of the sites. The model variant with both VPD and Rg terms ("+VPD+Rg") was only slightly better than considering only the Rg term: Superior at 62% of the sites (Table 1), but with a small mean MEF diff of 0.01 (Table 2).

Figure 5. Correlation of the parameter estimates for r and uWUE across sites (a). Distribution of the parameter correlations split according to aridity state of the site (b).
The difference between the +Rg-and +VPD-model is equally distinct, when the presence of a remaining model bias was 10 tested. For this, the intercept was estimated with all four model variants (Zhou, +VPD, +Rg, +VPD+Rg). Only the models with an Rg-term had considerably reduced residual intercepts (Fig. 4b).
All models compared to the original definition Zhou had two or more parameters. Of those models with two parameters, the +Rg emerged as the best model after cross-validation, indicating that the additional model complexity was in fact justified.
The presence of two parameters raised the question whether they could be identified independently. In fact, we found a high degree of correlation between the parameters for all sites (Fig. 5a). It is likely that the correlation originates from the 5 correlation between GPP and Rg. This is confirmed by noting that the correlation was lower at sites classified as arid (Fig.   5b), where the water limitation of GPP reduces correlations with radiation.
As transpiration of the +Rg model is a linear combination of stomatal and radiation-driven components, it is possible to calculate the relative contribution of each component to daily ET fluxes. We compared two approaches to calculate this 10 quantitiy: Parallel for estimating both parameters concurrently, hierarchical for estimating r only after uWUE has been calibrated (Eq. 10, 11). In both approaches, we observed that a sizeable fraction of mean daily ET could be attributed to the radiation term (Fig. 6). For the hierarchical approach, the mean global ET frac was 24% (95% confidence interval 21…28%); for the parallel approach, the mean global ET frac was 45% (95% confidence interval 40…49%).

15
Notably, ET frac estimated with the parallel approach varied widely between the sites (Fig. 7a). We assessed to which degree this variability can be interpreted as between-site variability of the expected value of ET frac or whether it is due to poorly constrained and correlated parameters (Fig. 5a). The 'within-site' uncertainty of ET frac caused by parameter uncertainty was quantified as the range between the 97.5 and the 2.5 percentiles (95% confidence interval, CI) of ET frac estimates from 200 parameter vectors sampled from their respective posterior distributions. The majority of sites had an CI lower than 0.3 (Fig.  20 7b). This suggests that the large variability of ET frac was not a result of parameter uncertainties. The conducted ANOVA (see Supplement S1) supported this conclusion, as it revealed that 95% of the global ET frac variability could be attributed to the variability between sites.
Biogeosciences Discuss., doi: 10.5194/bg-2016-524, 2017 Manuscript under review for journal Biogeosciences Published: 5 January 2017 c Author(s) 2017. CC-BY 3.0 License.     Finally, we assessed the impact of ET res and Rg as additional covariates on the global variability of uWUE estimates (Fig. 8).
We calculated Kendall's t rank correlation coefficient (Kendall, 1938) between the site-level estimates of uWUEs derived from the Zhou model and two different variants: +ET res and +Rg. The degree of correlation of the uWUE estimates quantifies whether changes in the model structure permutes the ordering of the estimates between sites. The extent to which 10 parameter estimates are affected by the model structure is crucial because any explanation or prediction of parameter values between sites would be highly desirable. While a moderate correlation between the uWUE estimates of the Zhou model with the +ET res -variant can be seen (t = 0.75), the correlation of the uWUE estimates of the Zhou model with the respective values of the +Rg-variant is low (t = 0.49). The between-site variability of uWUE explained by the estimates of the Zhou variant was 77% for the +ET res -variant and 34% for the +Rg-variant, as quantified by the R 2 . 15

20
In this study, we identified radiation as an important variable for ecosystem-scale transpiration and water-use efficiency.
Depending on the approach used, we attributed between a quarter and half of mean daily transpiration of all included Biogeosciences Discuss., doi: 10.5194/bg-2016-524, 2017 Manuscript under review for journal Biogeosciences Published: 5 January 2017 c Author(s) 2017. CC-BY 3.0 License.
FLUXNET sites to a linear radiation term. These findings raise the question which biophysical or ecophysiological processes can account for the estimated magnitudes of this attributed fraction.
The influence of radiation on stomatal conductance has been noted and discussed in the literature (Whitehead et al., 1981;Jarvis and McNaughton, 1986) and is also reflected in existing transpiration models (Leuning et al., 2008). By contrast, we 5 detected a substantial transpiration component that was statistically independent from the product of GPP and VPD 0.5 , tentatively suggesting an insensitivity to stomatal conductance. Our results suggest that this additional flux could not be associated with a g 0 conductance term, as this would imply the dependency of the additional ET on a linear VPD term (Ball et al., 1987). By contrast, the intercept was shown to be more consistent with radiation. Consequently, models that integrated radiation-driven transpiration with such an additive linear response had a superior predictive performance across flux towers. 10 The observed effect of radiation indicates that equilibrium transpiration could play an important role in ecosystem-scale transpiration. Equilibrium transpiration (or equilibrium evaporation rate) is the transpiration occurring if the leaf is completely decoupled from the atmosphere (Jarvis and McNaughton, 1986). Therefore, equilibrium transpiration is independent from stomatal conductance and driven solely by the dissipation of energy, provided by incoming solar radiation. 15 This can be contrasted with a competing explanation that the radiation dependency of transpiration does not reflect an additional process but rather a systematic problem of the VPD observations. VPD is measured together with the fluxes above the canopy. While the recorded water and carbon fluxes do in fact represent the net fluxes of the tower footprint, the same cannot be said for VPD. Its measurement above the canopy may differ substantially from the relevant magnitude of the variable at the leaf-scale. However, leaf temperature and thus the VPD of the leaf boundary layer is dependent on solar 20 radiation (Tenhunen et al., 1990). Adding solar radiation to the equation could therefore be seen as compensating for the lack of the aforementioned leaf-scale VPD observations. One limitation stems from the selection of rain-free periods for the parameter estimation. As previously described, this is a necessary step to justifiably assume that the observed latent heat fluxes constitute mostly a transpiration flux, rather than 25 evaporation from bare-soil and leaf surfaces in addition to transpiration. It also makes our work comparable with the study of Beer et al. (2009) that used this method to derive their estimates. For the observed residual evapotranspiration ET res we could show that it is not an artefact of insufficient exclusion of days after precipitation events. However, the environmental conditions during and after rain events generally represent some of these specific conditions: Low VPD due to the moisture available for evaporation, a higher share of diffuse radiation and a more or less sudden increase in soil moisture among 30 others. All mentioned variables could plausibly be assumed to have an influence on the stomatal opening of plants.
Therefore, the presented WUE model must not necessarily predict flux relationships during or immediately after precipitation events due to the underrepresentation of similar conditions in the sample of observations used for parameter estimation.
The effect of atmospheric CO2 concentrations on WUE (Morison, 1985;Conley et al., 2001) was not considered in this analysis. One could suspect that the seasonal variability of ET res is affected by seasonal CO2 variability if both were in phase. However, ET res showed a global maximum during June-July, while northern hemisphere CO2 concentrations are at its minimum in September-October (Keeling et al., 1976), implying that CO2 concentrations are unlikely to cause the seasonal 5 variation of ET res .
The model structures we tested in this study were evaluated according to their global empirical adequacy. Thus, despite the possible identification of probable mechanisms responsible for the observed patterns, the model structure selected in the end likely does not reflect the exact physical mechanisms by which the ecosystem operates. Furthermore, this limitation can be 10 specifically important in the case of water limitation. As the Rg-term is completely independent from any variable reflecting vegetation activity, our model would predict transpiration scaling with radiation during periods of severe drought. By principle, the same problem could affect periods of low temperatures before leaf flushing, although these periods are generally also associated with low radiation levels and hence may not be as problematic.

Implications
Our empirical analysis suggests that ecosystem-scale transpiration depends to a sizeable degree on variables other, which seem to be linked to radiation. This implies that photosynthesis and transpiration might be less strongly coupled on 20 ecosystem scale than commonly assumed. We speculate that the additional effect of radiation could be due to "equilibrium transpiration" where radiation drives transpiration even when the canopy is fully decoupled from the atmosphere. However, we cannot disentangle direct physiological effects of radiation on transpiration from other radiation effects that are relevant on ecosystem scale at this point. For the latter, leaf to canopy scaling, micrometeorological conditions, and boundary layer dynamics might contribute to the observed relationship between ecosystem scale water use efficiency and radiation. Thus, 25 further research is needed to reconcile our empirical findings with detailed ecosystem-scale modeling and theory on the one hand and plant-physiological research under controlled conditions on the other hand.
Finally, we caution against prematurely interpreting the between-site variability of underlying water-use efficiency (uWUE).
We showed that estimates of uWUE derived with the Zhou model explained only a third of the observed variability of 30 ecosystem uWUE derived from our empirically superior model formulation. The dependence of uWUE estimates on the chosen model formulation makes the interpretation of uWUE as an ecosystem property problematic. This concern would Biogeosciences Discuss., doi: 10.5194/bg-2016-524, 2017 Manuscript under review for journal Biogeosciences Published: 5 January 2017 c Author(s) 2017. CC-BY 3.0 License. also hold for assessing temporal dynamics of uWUE such as long-term trends. For example, the unexpectedly large global trend of WUE across FLUXNET sites (Keenan et al., 2013) would need to be tested for its omission of radiation in the model that was used. Overall, this study highlights the importance of model structure uncertainty for interpretations of parameter variability.

Data Availability
For this study, we used observations of the FLUXNET initiative from sites with an open and fair use data policy 2 . The data sets can be downloaded at http://fluxnet.fluxdata.org//data/download-data/