Examining the link between vegetation leaf area and land-atmosphere exchange of water, energy, and carbon fluxes using FLUXNET data

Vegetation regulates the exchange of water, energy, and carbon fluxes between the land and the atmosphere. This regulation of surface fluxes differs with vegetation type and climate, but the effect of vegetation on surface fluxes is not well understood. A better knowledge of how and when vegetation influences surface fluxes could improve climate models and the extrapolation 15 of ground-based water, energy, and carbon fluxes. We aim to study the link between vegetation and surface fluxes by combining yearly average MODIS leaf area index (LAI) with flux tower measurements of water (latent heat), energy (sensible heat), and carbon (gross primary productivity and net ecosystem exchange). We show that the correlation between LAI and water and energy fluxes depends on vegetation type and aridity. In water-limited conditions, the link between LAI and water and energy fluxes is strong, which is in line with a strong stomatal or vegetation control found in earlier studies. In energy20 limited forest we found no link between LAI and water and energy fluxes. In contrast to water and energy fluxes, we found a strong spatial correlation between LAI and gross primary productivity that was independent of vegetation type and aridity. This study provides insight into the link between vegetation and surface fluxes. It indicates that for modelling or extrapolating surface fluxes, LAI can be useful in savanna and grassland, but LAI is only of limited use in deciduous broadleaf forest and evergreen needleleaf forest to model variability in water and energy fluxes. 25


Introduction
Vegetation and water, energy, and carbon fluxes are tightly coupled. Large-scale vegetation patterns are driven by the longterm memory of water and energy availability (Köppen, 1936;Prentice et al., 1992;Cramer et al., 2001). Recent climate change leads to shifts in the spatial distribution of vegetation, as well as shifts in the timing of the growing season (Jeong et al., 2011;Rosenzweig et al., 2008;Fei et al., 2017). Additionally, vegetation plays a crucial role in the exchange of water, energy, and 30 carbon between the land surface and the atmosphere, mainly through its effects on evapotranspiration, turbulence, redistribution of water, and surface heating (Shao et al., 2015;Jia et al., 2014;Esau and Lyons, 2002). Large-scale reforestation and afforestation increased evapotranspiration over most of Europe , and large-scale deforestation increased the air temperature in tropical regions and decreased air temperature in boreal regions (Perugini et al., 2017). This two-way interaction between vegetation and terrestrial surface fluxes has been known for a long time (e.g. Bates and Henry, 35 1928;Woodwell et al., 1978), but is still a very relevant research topic today (Forkel et al., 2019;Lu et al., 2019; Hoek van Dijke, 2020; Kirchner et al., 2020;Evaristo and McDonnell, 2019), given the importance of understanding the impacts of climate change on vegetation, as well as the effects of land cover change on climate.
Plants regulate the exchange of water, energy, and carbon with the atmosphere through their stomata. The stomatal regulation 40 of these fluxes depends on available energy, transpiration demand, and available soil moisture in the root zone. When both the available energy and soil moisture are abundant, stomata open and water and carbon can freely move in and out: the stomatal control on surface fluxes is low. When the available energy is high, but soil moisture is limiting, stomata tend to close and exert a large control on water and carbon fluxes (Mallick et al., 2016;O'Toole and Cruz, 1980). Zooming out from stomatal to canopy scale, there are several other ways in which vegetation influences surface fluxes. Soil and crown mutual shadowing 45 and deep ground water uptake by vegetation influence the latent heat flux whereas soil moisture influences ecosystem respiration and thereby carbon exchange (Chen et al., 2019;Schmitt et al., 2010). The vegetation control of ecosystem fluxes has been shown by different data or modelling studies and depends on climate and vegetation type (Williams et al., 2012;Xu et al., 2013;Wagle et al., 2015). Williams and Torn (2015) found a strong vegetation control on surface heat flux partitioning in both arid and humid grassland, cropland, and forest, but Padrón et al. (2017) concluded that globally, vegetation control on 50 evapotranspiration was low and even absent in the equatorial regions. Chen et al. (2019) showed that for wetland sites, temperature, precipitation and vegetation leaf area explained 91% of the mean annual variability in vegetation carbon uptake. Mallick et al. (2018) showed that vegetation control on evapotranspiration was stronger in arid ecosystems as compared to the mesic ecosystems. Similar results were found for dry and wet Amazonian forest (Costa et al., 2010;Mallick et al., 2016) and dry and wet grassland (De Kauwe et al., 2017). Ferguson et al. (2012) studied land-atmosphere coupling of fluxes, which 55 includes the effect of vegetation as well as other factors as soil wetness, soil texture, and surface temperature. From remote sensing data and model output, they concluded that transitional zones between arid and humid climates (shrublands, grasslands, and savannas) tend to have a strong land-atmosphere coupling, while in the energy-limited regions, land-atmosphere coupling is weak.
Vegetation is coupled to the atmosphere through its leaves. The leaf area index (LAI) is an important vegetation characteristic and is indicative of the total amount of foliage that intercepts light and assimilates carbon. Furthermore, both rainfall interception and canopy conductance increase with LAI (Van Heerwaarden and Teuling, 2014;Gómez et al., 2001). A high LAI is therefore related to high vegetation carbon uptake and high canopy evapotranspiration of water (Lindroth et al., 2008;Duursma et al., 2009). Highest mean yearly LAI is found in tropical and temperate forests, while a low LAI is found in cold 65 and in arid climate zones (Iio et al., 2014;Asner et al., 2003) (Figure 1). This global LAI pattern closely resembles large-scale patterns in estimates of water, energy, and carbon exchange (Miralles et al., 2011;Jung et al., 2011). With an increasing availability of remotely sensed LAI data, LAI − besides its usage in many remote sensing applications (e.g. Si et al., 2012;Zheng and Moskal, 2009) − became a frequently used variable to represent vegetation in land-surface models (Williams et al., 2016;Sellers et al., 1997;Lawrence and Chase, 2010 amongst many others) or to estimate or extrapolate regional or global 70 water and carbon fluxes (Beer et al., 2007;Yan et al., 2012;Turner et al., 2003;Xie et al., 2019). The algorithms to retrieve LAI from remotely sensed data improved during the past decades, increasing the accuracy of LAI products (Shabanov et al., 2005;Yan et al., 2016). Nevertheless, it is important to be aware of the product uncertainties, especially over dense forest, where saturated reflectance and canopy clumping can only provide limited information for LAI retrievals (Shabanov et al., 2005;Xu et al., 2018), and at high latitudes, where the solar zenith angle is low (Fang et al., 2019). 75 The interaction between vegetation LAI and surface fluxes on larger scale is not yet well understood and vegetation is not well represented in many land-atmosphere and climate models (Williams et al., 2016). A small scale study in temperate deciduous forest, for instance, revealed that the correlation between sap flow and the normalized difference vegetation index (NDVI) can change from positive to negative depending on the season and soil moisture availability . A 80 detailed knowledge of how and when vegetation LAI is linked to the surface fluxes is required to improve global climate modelling and extrapolation of water and carbon fluxes from canopy to ecosystems. The high availability of remote sensing LAI products, recent developments in cloud-based platforms for geospatial analysis (Mutanga and Kumar, 2019), and the availability of publicly available eddy covariance data from FLUXNET (Baldocchi et al., 2001) allows for an analysis of the link between vegetation characteristics and surface fluxes. The objective of our study is to get an insight about the link between 85 vegetation LAI and surface fluxes for different vegetation types along an aridity gradient. We address the following research questions: 1) What is the link between LAI versus water, energy, and carbon fluxes in different vegetation types? 2) How is the interaction between LAI versus water, energy, and carbon fluxes governed by climatological aridity? We hypothesise that the link between LAI and surface fluxes is strong in semi-arid and arid climates, owing to the strong stomatal control, while the link is weak in humid climates. 90 In our study we focus on five metrics of water, energy, and carbon fluxes measured by flux towers. Latent heat (LE), a measure for the evapotranspiration of water, and sensible heat (H), represent the exchange of water and energy between the Earth's surface and the atmosphere. LE and H are linked through the evaporative fraction (EF). The EF is the ratio of latent heat to the sum of LE and H and is a useful measure of the partitioning of total available energy between the evapotranspiration of water 95 and surface heating. Net Ecosystem Exchange (NEE) is the net exchange of carbon between the land and the atmosphere, which is directly measured by flux towers. Gross primary productivity (GPP) is derived from NEE and is the gross uptake of atmospheric carbon by the vegetation. Within the FLUXNET-2015 dataset (Baldocchi et al., 2001), we selected all Tier-1 sites (open and free for scientific purposes) within the five studied vegetation types. We completed the dataset with two sites from the OzFLUX network to increase the 110 number of sites in the EBF class (Liddell, 2013b, a). Two forest sites were excluded from the analyses because they were effected by a beetle outbreak that resulted in high tree mortality, and one heavily managed grassland site was excluded from the analysis. For each site, only years with good-quality data were selected, following the quality selection procedure that is explained below. This site selection procedure, in combination with the quality check, resulted in a dataset of 545 site-years spread over 93 sites ( Figure 2, Table 1). 115

Data averaging and aggregation
We studied yearly averaged LAI and surface fluxes for different vegetation types. In most vegetation types, LAI and surface fluxes showed seasonal variability, with high values during the growing season and lower or zero LAI and surface fluxes during the cold or dry season. The non-growing season might be non-relevant for finding the link between LAI and surface fluxes, however, selecting growing season values only lead to difficulties. The vegetation types differ in the timing, number, 120 and length of growing seasons, and for instance time-series analysis did not successfully select the growing seasons. To be consistent in the methodology, yearly averaged fluxes were used for all flux tower sites. Using yearly averaged values for every site (referred to as 'site-years') has few implications 1) we study both spatial (site-to-site) variability and temporal (year-to-year) variability simultaneously, and 2) averaged flux and meteorological measurements might not represent similar conditions. The latter is for example when a site-year receives plenty of precipitation in December, increasing the site-year's 125 aridity index, while this precipitation mainly impacts the next site-year's fluxes or LAI values. To test the effect of using siteyear data, we also studied spatial and temporal variability separately. For these analyses, the data was aggregated in three ways: 1) Site-year data, having one average value per site per year, 2) multi-year data, having one multi-year average LAI and flux value per site, to study spatial correlation, and 3) yearly average data for a few sites, to study the temporal correlation.
Sites were included in the multi-year data if at least three years of data were available. The three aggregation methods led to 130 similar conclusions for water and energy, but slightly different results for carbon, as is shown in the manuscript.

Flux measurements
Within the FLUXNET 2015 database, LE, H, NEE, and GPP measurements are gapfilled using the MDS (Marginal Distribution Sampling) method (Reichstein et al., 2005), and LE and H are corrected by an energy balance closure correction factor. The MDS method uses the correlation of fluxes with the driver variables (incoming radiation, temperature, and vapour 135 pressure deficit) to estimate flux values during gap periods. The energy balance closure corrects LE and H for the total incoming radiation, assuming that the Bowen ratio (the ratio of the sensible heat flux to the latent heat flux) is correct. A similar energy balance closure correction was applied to the LE and H measurements of the OzFLUX sites. Monthly averaged flux values were discarded if the percentage of measured and good quality gapfill data was below 50%. Yearly average fluxes were calculated if measurements for each month were available. The evaporative fraction (EF), the ratio between LE and the total 140 energy available at Earth's surface was calculated using Eq. (1) as follows: where LE is the latent heat flux and H is the sensible heat flux.

Meteorological measurements
Meteorological measurements are delivered with the flux tower data. Precipitation data is downscaled from the ERA-interim 145 reanalysis data (Vuichard and Papale, 2015). Net radiation and air temperature are measured at the flux tower and gap-filled using the MDS (Marginal Distribution Sampling) method (Reichstein et al., 2005). Yearly potential evaporation (Ep) was calculated from mean daily air temperature and net radiation using the Priestley-Taylor formulation (Priestley and Taylor, 1972). The Priestley-Taylor equation is a modification of the Penman equation and requires less measurements. The aridity index (AI), an indicator of dryness, was calculated according to Eq. (2) 150 where P is precipitation and Ep is the potential evaporation. An aridity value of one indicates that, on a yearly scale, precipitation equals potential evaporation, while values below one indicate site-years that received less precipitation than their potential evaporation. 155

Leaf Area Index
Leaf Area Index (LAI) is the ratio of green leaf area to ground area (in m 2 m -2 ). We used LAI derived from the MODIS data product MCD15A3H.006 (Myneni et al., 2015). This algorithm derives 4-day composite LAI values on 500 m spatial resolution from the Terra and Aqua satellites and is available for 2003 onwards. Within this 4-day period, the best pixel is selected from the MODIS sensors located on the Terra and Aqua satellite for the calculation of LAI. The LAI calculation 160 algorithm uses a Look-up-Table that was generated using a 3D radiative transfer equation (Myneni et al., 2015). Heinsch et al. (2006) compared the MODIS data product with ground measurements at FLUXNET sites and concluded that 62.5% of the MODIS LAI was well estimated, but that MODIS LAI overestimated ground measured LAI for the other sites. Despite this overestimation, MODIS LAI was used, because it has a long record length, good (and free) data availability, good spatial coverage, and high temporal resolution. The overestimation and saturation of the signal at high LAI could introduce noise in 165 the LAI data. We do however not expect this noise to change the conclusions of our analysis. The resolution of the LAI data product is 500 m, compared to a typical flux tower footprint length of 100 to 1000 m (Kim et al., 2006). The exact size and location of the footprint of flux towers however varies with among others wind direction and wind speed, surface roughness, and flux measurement height (Kim et al., 2006;Barcza et al., 2009). For our analyses, we selected the one nearest LAI pixel for each flux tower. Data were filtered to remove clouds, using the with the product delivered quality label. To smoothen 170 outliers, the moving mean LAI was calculated for three consecutive data points. Monthly mean values were calculated if at most one data point was missing. Site-year average LAI was calculated when no monthly data were missing.

Methodology
To study the link between LAI and surface fluxes, we performed a linear regression between LAI and the surface fluxes. We calculated the correlation coefficient for 1) site-year data, 2) multi-year average data (spatial variability) and 3) yearly data for 175 a few specific sites (temporal variability). Afterwards, to study if the link between LAI and fluxes changed with aridity, all site-years within one ecosystem type were ranked by aridity, from most arid to most humid. For each consecutive 30 site-years in this ranking, we performed a linear regression between LAI and the fluxes. For some site-years, part of the data was missing that was needed to calculate the regression. Within each window of 30 site-years, the slope of the regression was calculated if at least 15 complete site-years were available (Figure 3). 180

Figure 3 Illustration of the applied methodology. The correlation coefficient between leaf area index (LAI) and evaporative fraction (EF) is calculated for 30 site-years for grassland over a moving window of aridity index.
In the illustration, the correlation has a significant positive slope at p = 0.056 for the 30 most arid grassland sites, while for the 30 most humid grassland sites, the slope is nearly flat and not significant (p = 0.49).

The link between water, energy, and carbon fluxes versus LAI
LAI and LE were positively correlated in SAV, GRA, and EBF ( Figure 4, Table 2). The slope of the correlation between the different vegetation types is different; the slope was steepest for SAV (slope = 46.1 W m -2 ): a doubling in LAI (1 to 2) was associated with almost a doubling in LE (51 to 97 W m -2 ), compared to a flatter slope in GRA (9.80 W m -2 ) and EBF (13.0 W m -2 ). In ENF and DBF, LAI and LE were not significantly correlated. LAI and H were negatively correlated in SAV, GRA 190 and EBF, while there was no significant correlation in ENF and DBF. LAI and the EF were positively correlated in SAV, GRA and EBF, while no correlation was found in ENF and DBF. A positive slope indicates that, for a higher LAI, a higher fraction of the available energy is used for evapotranspiration of water, compared to surface heating. The slope between LAI and EF was steeper in SAV and GRA (slope = 0.27 for both) than in EBF (slope = 0.08). A positive correlation between LAI and GPP was found in all vegetation types (r = 0.47 -0.97), with a very strong correlation coefficient for SAV (r = 0.97). The correlation 195 followed a steep slope for SAV (slope = 3.37 gC m -2 d -1 ) and GRA (slope = 2.17 gC m -2 d -1 ), a similar slope in EBF (slope = 1.71 gC m -2 d -1 ) and ENF (slope = 1.81 gC m -2 d -1 ), and a less steep slope in DBF (slope = 0.76 gC m -2 d -1 ). The correlation between LAI and NEE is negative in SAV, EBF, and ENF. This indicates that net carbon uptake increases with LAI. Among the different fluxes, GPP showed the strongest correlation with LAI for all vegetation types. Comparing the different vegetation types, the correlation between LAI and fluxes was strongest in SAV. 200 Using multi-year average data reduced the number of data points to only 5 to 16 sites per vegetation type. Nevertheless, the spatial correlation (site-to-site variability) between LAI and surface fluxes is very similar to the spatio-temporal correlation ( Figure 5, Table 2). For SAV, GRA, and ENF, the slope and strength of the correlation were similar when compared with the site-year data. For the EBF, for the site-year data, the correlation with LE and EF was only significant at p ≤ 0.1 and the 205 correlation was not significant for H and NEE. Table 2 Strength and significance of the correlation between LAI versus surface fluxes for site-year and multi-year average data. The correlation coefficients are shown for significant correlations at p ≤ 0.05 (*) or at p ≤ 0.1 (·). A -indicates that the correlation was not significant. Temporal (year-to-year) variability in LAI and surface fluxes was smaller than spatial (site-to-site) variability ( Figure 6). For both SAV sites, and one of the two GRA, EBF, and DBF sites, LAI and LE were positively correlated in time. For H, one EBF site showed a significant negative correlation with LAI, and for EF, and one of the two SAV, GRA, EBF, and DBF sites showed a positive correlation with LAI (p ≤ 0.1 or p ≤ 0.05). For GPP and NEE, one of the SAV, GRA, EBF, and ENF sites 215 showed a positive correlation, and for NEE. Overall, the temporal correlations between LAI and surface fluxes was of similar direction as the spatio-temporal and spatial correlations. For more than half of the sites in Figure 6, however, year-to-year variability in LAI and surface fluxes was low and variability in fluxes was not significantly correlated with variability in LAI.

Figure 6 An illustration of the temporal correlation between yearly average surface fluxes and leaf area index (LAI).
For each land cover type, two sites were selected that had the highest number of available data. The colours of the symbols indicate the land cover type as in Fig 4 and Fig 5. Panels show (a) Figure 7 shows the steepness and significance of the correlation between LAI and surface fluxes for different aridity values.

The effect of climatological aridity on the link between LAI and surface fluxes 220
In dry vegetation types or regions, the correlation between LAI and fluxes was significant and had a steeper slope, while in the more humid vegetation types or regions, the slope was relatively horizontal and the correlation was often not significant. In SAV, GRA, and EBF, the correlation between LAI and LE was significant for the whole range of aridity values. In arid GRA, the correlation had a steeper slope, as compared to humid GRA. For LAI versus H and LAI versus EF, the slope was steep and 225 significant for SAV. For GRA, the correlation was strong and significant in the arid regions, and insignificant for the humid regions. For EBF, the slope and significance of the correlation did not change with aridity. For LAI and GPP, the slope and significance of the correlation did not change with aridity for SAV, GRA, EBF, and ENF. For DBF, the correlation between LAI and GPP was negative at higher aridity, but these results were strongly influenced by one site with an above average LAI for all the site-years. For LAI versus NEE, a steep slope with negative correlation was found in arid SAV and humid ENF. In 230 other humid regions, the correlation was less steep.

Figure 7 The effect of aridity on the relation between surface fluxes and leaf area index (LAI). The slope of the correlation between LAI and surface fluxes is shown for different aridity values for (a) the latent heat flux (LE), (b) the sensible heat flux (H), (c) the evaporative fraction (EF), (d) gross primary productivity (GPP), and (e) net ecosystem exchange (NEE). Each dot indicates the slope value for the 30 closest aridity values. The filled symbols indicate that the correlation was significant at p < 0.05, while the empty symbols indicate a non-significant correlation.
To study how the correlations varied with climatic drivers of surface fluxes, we calculated the correlation coefficient between the fluxes versus precipitation (P) and incoming shortwave radiation (Rg) (Figure 8). In SAV, GRA, and EBF, the water fluxes showed a strong correlation with P, indicating that water availability partly explained the spatio-temporal variability in surface fluxes. In ENF and DBF, there was a weak or no correlation between LE and P, but a strong correlation with Rg. This indicates 235 that available radiation was the primary driver of water and energy fluxes in these sites.

Discussion
The EBF site-years span a wide range of LAI values (LAI = 0.9 -6.1) and aridity conditions (AI = 0.3 -9.3), and both are a potential limitation of our analysis for the EBF vegetation type. The uncertainty of the LAI retrieval in dense vegetation is higher compared to other vegetation types due to saturation of the remotely sensed signal. The large range of climatic 240 conditions indicates that our EBF site-years range from arid, water-limited conditions to humid conditions. Despite this high variability in site-years, the sites fell within one vegetation type.
The correlation between LAI versus water and energy fluxes (LE, H, and EF) varied with vegetation type and aridity. For the spatio-temporal and spatial variability, we found 1) strong (positive or negative) correlations and (partly) steep slopes for SAV 245 and GRA, 2) a significant correlation, but less steep slope for EBF, and 3) no significant correlations for ENF and DBF. For the temporal variability, this pattern was similar for LE, but almost no significant correlations were found for LAI versus H and EF for SAV and GRA. Evapotranspiration is the sum of transpiration, soil evaporation and interception evaporation and the magnitude of each component depends on LAI. Transpiration increases with LAI at the cost of soil evaporation when there is sufficient moisture available (Gu et al., 2018;Wang et al., 2014). In arid climates, the transpiration component is higher 250 compared to wetter climates (Gu et al., 2018) and the link between transpiration and LAI is particularly strong in these arid climates (Sun et al., 2019). When soil moisture is deficient and vegetation encounters a high evaporative demand, stomatal control is stronger (Mallick et al., 2016). This accelerates a strong stomatal coupling between LAI and LE and could explain the strong correlation between LAI versus LE, H, and EF that was found in SAV and arid GRA. Soil water deficiency and high evaporative demand leads to a high increase in LE, for a small increase in LAI, which could explain the steep(er) slope 255 in arid GRA and SAV vegetation.
In forests, soil evaporation is low, while interception evaporation is large. The high interception evaporation is due to the large leaf area (both green leaves included in the LAI and brown leaves after leaf senescence) with a high canopy water storage capacity and a high turbulence, enhancing fast evaporation (De Jong and Jetten, 2007). In EBF, interception evaporation contributes to up to 30% of total evapotranspiration (Wei et al., 2017;Gu et al., 2018). This could explain the strong correlation 260 between LAI versus water and energy fluxes in EBF. A high interception evaporation was however also reported for temperate and boreal forest (Miralles et al., 2011), while for these forest types, we found no correlation between LAI and water and energy fluxes. The ENF and DBF sites were found in humid regions, and fluxes were in the first place energy-limited. In these energy-limited sites, LAI played no, or a weak role in controlling surface fluxes. This indicates a weak or no vegetation control on surface water and energy fluxes in energy-limited sites. This is in line with a low land-atmosphere coupling in energy-265 limited sites (Ferguson et al., 2012).
In contrast to the results for water and energy fluxes, the spatio-temporal and spatial correlation between GPP versus LAI was strong across all vegetation types and (almost) all aridity gradients. A strong link between LAI and carbon uptake on yearly timescale over all vegetation types is expected, as plants try to optimize carbon gain and would generally not display leaves 270 with a negative carbon balance. A strong link between LAI and mean yearly GPP was also shown by Hashimoto et al. (2012).
Other studies however found a weak link between LAI and GPP for annual time scales (Law et al., 2002). In contrast to the spatial variability, year-to-year variability in GPP was only in part of the sites correlated to LAI. Water availability is an important driver for temporal variability in GPP (Williams and Albertson, 2004;Kutsch et al., 2008), and GPP is strongly reduced under drought conditions (Vicca et al., 2016). The effect of drought is also visible in reduced LAI, but on a longer 275 time scale of one or two years in forest (Le Dantec et al., 2000;Kim et al., 2017). This different response time to water availability for forest LAI and GPP could partly explain the absence of a temporal correlation for part of the sites. The spatial correlation between LAI and NEE was less strong as compared to GPP, which is in agreement with results of Chen et al. (2019). NEE is the sum of carbon uptake by the vegetation (GPP) and carbon loss by ecosystem respiration. Ecosystem respiration varies with climate and soil carbon storage, which are not directly related with LAI. This could explain the absence 280 of a correlation between LAI and NEE.
The results partly confirmed our hypothesis. As hypothesised, the correlation between LAI and surface fluxes was strong in arid regions for water and energy fluxes, and the correlation was absent in humid ENF and DBF. For humid EBF, however, we found a strong correlation between LAI and water and energy fluxes, and for GPP, the correlation with LAI was strong 285 across all aridity gradients. While carbon uptake is the primary goal of vegetation, independent of the aridity gradient, ecosystem water loss comes inevitably with carbon uptake, but also depends on vapour pressure deficit, available radiation, and soil moisture, which are not directly linked to LAI.
Our statistical analysis cannot be used to study causality between LAI and surface fluxes, or to study vegetation control on the 290 surface fluxes. The correlation between LAI and water fluxes is confounded by the effect of soil moisture, especially in arid and semi-arid ecosystems, where both canopy development and LE increase with water availability (Kergoat, 1998;Mallick et al., 2018). Similarly, precipitation is the main controller for spatial variability in both vegetation and GPP (Koster et al., 2014). Furthermore, LAI is related to vegetation properties, but not a direct measure of canopy conductance. Despite, there are similarities with previous studies showing the stomatal or vegetation control on surface fluxes. A strong vegetation control on 295 water and energy fluxes in arid and semi-arid regions was shown on timescales of days or smaller (e.g. Mallick et al., 2016;Mallick et al., 2018) and also our study shows that, on large spatio-temporal scale, LAI versus water and energy fluxes show the strongest correlation in arid regions. For EBF however, we found a strong spatial correlation between vegetation versus water, and energy fluxes, while Padrón et al. (2017) showed that vegetation control in equatorial regions was absent. An interesting follow-up study would be to link stomatal control for different vegetation types (De Kauwe et al., 2017) to the 300 canopy-scale pattern investigated in this study.
Our analyses give insight in how and when vegetation LAI is related to surface fluxes. The results show that LAI is a good predictor for spatial variability in GPP across different vegetation types and aridity gradients. Furthermore, the analysis suggests that, in SAV, GRA, and EBF, LAI could be used to describe canopy-scale spatio-temporal variability water and 305 energy fluxes. LAI is however not a good predictor for water and energy fluxes in ENF and DBF and for NEE. It is important to be aware of these limitations when using LAI to describe or estimate water, energy, and carbon fluxes in climate models or extrapolation methods. This study provides insight in the link between surface fluxes and LAI and could be used to improve predictions of the effects of land cover change on surface fluxes.

Conclusions 310
The objective of this study was to get an insight about the link between vegetation LAI and land-atmosphere fluxes for different vegetation types along an aridity gradient. We studied this link at large spatio-temporal scale using flux tower measurements of water, energy, and carbon, combined with satellite derived LAI data. The data analysis led to the following conclusions: a) The link between LAI versus water and energy fluxes depends on vegetation type and aridity. The correlation between LAI versus water and energy fluxes is strong in SAV, GRA, and EBF. In DBF and ENF however, no significant 315 correlation was found. Contrary to water and energy fluxes, the spatial correlation between LAI versus GPP was strong, independent of vegetation type and aridity. This suggests that using LAI to model or extrapolate surface fluxes of water and energy is well possible in SAV, GRA, and EBF, but is limited in DBF and ENF. b) As hypothesised, the link between LAI and water and energy fluxes was strong in arid, water-limited conditions and absent or weak for humid, radiation-limited conditions. EBF, which was found over a high range of aridity conditions, 320 but mostly in humid environments, forms an exception: the spatial correlation between LAI versus water and energy fluxes was strong, despite the overall humid conditions. This researchfacilitated by the recent availability of large global datasets of remotely sensed LAI, flux tower data, and cloudcomputing platformshas added to the understanding of LAI interaction with surface fluxes and could help to improve modelling or extrapolating surface fluxes. 325

Author contribution
The data analyses were done by AJHvD in close consultation with KM, MS, MM, MH, and AJT. AJHvD prepared the draft manuscript and all authors contributed to the discussions and writing of the manuscript.

Acknowledgements
This study was supported by the Luxembourg National Research Fund (FNR) (PRIDE15/10623093/HYDROCSI). We also acknowledge Prof. Michael Liddell for providing the data of two OzFlux research sites. We further acknowledge the FLUXNET community, for acquiring and sharing the eddy covariance data, including these networks: AmeriFlux, AfriFlux, AsiaFlux, CarboAfrica, CarboEuropeIP, CarboItaly, CarboMont, ChinaFlux, Fluxnet-Canada, GreenGrass, ICOS, KoFlux, 335 LBA, NECC, OzFlux-TERN, TCOS-Siberia, and USCCC. The FLUXNET eddy covariance data processing and harmonization was carried out by the European Fluxes Database Cluster, AmeriFlux Management Project, and Fluxdata project of FLUXNET, with the support of CDIAC and ICOS Ecosystem Thematic Center, and the OzFlux, ChinaFlux and AsiaFlux offices. The ERA-Interim reanalysis data are provided by ECMWF and processed by LSCE.