Methods

Study Site

The study was conducted at the Kakubari eddy covariance tower located in Santa Rosa National Park, Costa Rica (Figure 2). The site is representative of a tropical dry forest with marked seasonality, and strong intra-annual variability in atmospheric demand. This tower constantly measured the turbulent exchanges of carbon dioxide, water vapor, and energy between the ecosystem and the atmosphere. The workflow followed in this research is shown in Figure 3. 

Eddy Covariance Data Processing

Measurements every 30 minutes were processed using EddyPro software developed by LI-COR, from 2016 to 2024 with gaps produced by problems with sensors or source of energy. At this point, gaps-filling techniques were avoided to present original data from this ecosystem. Data were aggregated from 30 minutes to daily averages after filtering, and only complete daily observations were used to avoid artificial covariance structures through gaps-filling. 

Estimation of Ecosystem Carbon Flux Components

Gross Primary Productivity (GPP in µmol CO₂ · m-2 · s-1) was estimated subtracting the Net Ecosystem Exchange (NEE) from the ecosystem respiration (Reco):

GPP = Reco - NEE                                                                                               (1)

The ecosystem respiration was estimated using a temperature-dependent Q10 formulation (S. Chen et al., 2020):

Reco = Rref * Q10^(T-Tref)/10                                                                                       (2)

Where air temperature (T) was used as a driver of respiratory fluxes.

The water use efficiency (WUE in µmol CO₂ / mol H₂O) and the intrinsic water-use efficiency (iWUE in µmol CO₂ · kPa / mol H₂O) were approximated using equations (3) and (4) respectively (Liu et al., 2022):

WUE = GPP / ET                                                                                                 (3)

iWUE = (GPP * VPD) / ET                                                                                           (4)

Where ET is the evapotranspiration measured in mm/h, but it was transformed to mol H₂O · m-2 · s-1. 

We chose to use Vapor Pressure Deficit instead of relative humidity because relative humidity measures the percentage of water saturation in the air at a given temperature, while Vapor Pressure Deficit measures the air's actual "thirst," or the pressure difference (in Pa) between the current humidity and saturation. In forests, Vapor Pressure Deficit is a better indicator of plant transpiration.

Footprint Modeling 

Footprint modeling was made using the analytical parametrization formulated by Kljun et al. 2015, estimating the area of action of the tower through variables such as friction velocity (u*), wind direction and speed, canopy height, surface roughness length (z0), boundary layer height (h). The boundary layer height was extracted from the ERA5 Hourly reanalysis product, provided by the European Centre for Medium-Range Weather Forecasts. This variable is crucial in the Kljun model because it can affect the radial extent of the footprint model.

The footprint analysis was implemented to define the spatial representativeness of the tower measurements. It quantifies the specific area of the upwind forest canopy (the 'source area') that contributes to the observed gas exchanges. 

For this project, the surface roughness length was estimated as 10% of canopy height, while the displacement height (d) was 0.67 times canopy height. During the execution of Kljun model, a restriction of  zm/L >= -15.5  (where zm represent the measurement height above displacement height and L the Obukhov length) was applied to ensure stability conditions (Abdaki et al., 2024).

Remote Sensing Data Processing

Vegetation structural and spectral metrics within a one-kilometer buffer around the eddy covariance tower were obtained using the Google Earth Engine (GEE) cloud-computing platform (Figure 2). The buffer was selected based on an approximation of the tower’s footprint under unstable atmospheric conditions (a metric obtained with the Kljun model).

Canopy height average was obtained from the ETH Global Canopy Height Model (10 m resolution). The “reduceRegion” function was used to compute the zonal statistics, including mean, standard deviation, and percentiles.

Spectral vegetation indices (NDVI and EVI) were extracted using the same spatial buffer to ensure consistency between eddy covariance measurements and remotely sensed vegetation properties.

Homogeneity and Stability Assessment

A stability study was conducted to ensure that the tower measurements were spatially representative and functionally homogeneous. Ecosystem fluxes and water-use efficiency indicators were classified into three distance groups based on the daily X80​ footprint: short (<300 m), medium (300-600 m), and long (>600 m). Also, the influence of atmospheric transport was assessed by grouping data into wind speed classes (low, medium, and high) and directional sectors (North, East, West, and South). 

The Coefficient of Variation (CV%) was calculated for each group to ensure that the observed variability was driven by biological processes rather than sensor location or turbulent transport conditions. A low CV in Gross Primary Productivity across these gradients was used as the primary indicator of a functionally uniform canopy.

Statistical Analysis

Multivariate statistical analysis was applied to identify the dominant environmental driver and its influence on ecosystem fluxes. The Interquartile Range (IQR) method was applied to detect outliers, defined by 1.5 * IQR (below the first quartile or above the third quartile). These values falling outside the boundary were identified as potential instrumental noise or non-biological transients. To mitigate extreme values ​​and improve statistical calculations, logarithmic transformation and winsorization (2%–98%) were applied to the fluxes, stabilizing the variance to achieve normalization in the distribution. All variables were then standardized to zero mean and unit variance to perform principal component analysis. This analysis was used to obtain the dominant axes of variability and quantify variable loading through the correlation matrix. 

Canonical Correlation Analysis (CANCOR) was used to determine the direct relationship between two multidimensional sets: environmental predictors (Vapor Pressure Deficit, Air Temperature, Photosynthetic Active Radiation, and vegetation indices) and ecosystem responses (Gross Primary Productivity, Net Ecosystem Exchange, Latent Heat Flux, and water usage efficiency). This approach identified the highest correlation (Rc) between atmospheric forcing and biological reaction, resulting in a comprehensive understanding of the dry forest's functional dynamics.