Mature riparian alder forest acts as a strong and consistent carbon sink

Alder forests are widespread across the Northern Hemisphere, often occupying riparian zones and enhancing soil fertility through symbiosis with nitrogen-fixing bacteria. Despite their ecological importance, the ecosystem-level carbon and water exchange of alder forests remains poorly studied, particularly under contrasting hydroclimatic conditions. We studied ecosystem carbon and water fluxes over three contrasting years (“wet”, “drought”, “recovery”) in a mature riparian grey alder forest in Estonia. The forest was a strong and consistent net carbon sink with annual net ecosystem exchange (NEE) ranging from −496 to −663 g C m⁻² yr⁻¹, gross primary production (GPP) from −1258 to −1420 g C m⁻² yr⁻¹, ecosystem respiration (ER) from 595 to 923 g C m⁻² yr⁻¹ and evapotranspiration (ET) varied from 194 to 342 kg H₂O m⁻² yr⁻¹. Moderate soil water saturation (40 %–50 %) enhanced all ecosystem fluxes. In contrast, progressive drought reduced ER, ET, and to a much lesser extent GPP, with elevated EWUE and suppressed canopy conductance indicating strong stomatal regulation to limit water loss while maintaining carbon sequestration. While soil saturation affected canopy conductance, its effect was outweighed by vapour pressure deficit during the drought year, even after soil water availability recovered. We observed a full recovery in the following year, which was supported by favourable temperature and precipitation, although partially suppressed canopy conductance suggested some vulnerability to possible consecutive droughts in the future. Overall, the forest demonstrated drought resilience and high net carbon uptake across contrasting years, underscoring the capacity of riparian alder stands to sustain carbon sequestration under variable hydroclimatic conditions.

1 Introduction

Terrestrial ecosystems play an essential role in restraining the rise of atmospheric carbon dioxide (CO₂) concentrations and mitigating global warming (Pan et al., 2011; Piao et al., 2020). Over the preceding decades, they have sequestered approximately one-third of the total industrial carbon emissions (Friedlingstein et al., 2022). Forest ecosystems typically act as net carbon sinks, with the rate of photosynthetic uptake surpassing respiratory emissions on the annual scale (Harris et al., 2021). The strength of this carbon sink depends upon various factors, including, but not limited to, forest age, tree species composition, climatic conditions, soil properties, and management practices (Winkler et al., 2023). Moreover, a change in weather conditions or forest management decisions can turn a local carbon-sequestering forest stand into a net carbon source, thereby affecting ecosystem-atmosphere interactions at a regional scale (Hadden and Grelle, 2016; Lindroth et al., 1998). Thus, it is critical to evaluate the local forest carbon uptake in the face of varying climatic events (Allen et al., 2010; Bonan, 2008; Teskey et al., 2015).

Water availability plays a particularly critical role among the environmental factors affecting forest carbon uptake. Drought can reduce photosynthesis, increase tree mortality, and temporarily weaken or reverse a forest's sink function (Allen et al., 2010; Breshears et al., 2005; Cavin et al., 2013; Haberstroh et al., 2022; McDowell et al., 2008). The frequency and severity of extreme climate events, including droughts, have been growing in recent decades, a trend expected to continue (Fischer et al., 2021; Trenberth et al., 2014). The 2018 European drought was considered the most severe in the last 250 years (Gutierrez Lopez et al., 2021; Hari et al., 2020), resulting in a significant decline in forest carbon uptake and elevated tree mortality rates (Bastos et al., 2020; Buras et al., 2020; Haberstroh et al., 2022; Senf and Seidl, 2021; Smith et al., 2020). Improving our understanding of how forest carbon and water fluxes are modified during and after such events is crucial for assessing ecosystem resilience and informing adaptive forest management.

Riparian forests, located at the interface between terrestrial and aquatic systems, play a crucial role in mediating nutrient and carbon flows and are particularly sensitive to hydrological changes (Capon et al., 2013; Dybala et al., 2019; Naiman and Décamps, 1997). Grey alder (Alnus incana (L.) Moench.) is a typical pioneer species frequently occupying riparian zones and is widely spread in North America and Europe (Caudullo et al., 2017). Alder plantations can mitigate carbon losses in rewetted peatlands (Huth et al., 2018) and improve the soil structure of skid trails (Warlo et al., 2019). Their high adaptability also makes alders suitable for the afforestation of post-industrial sites (Krzaklewski et al., 2012). Owing to their symbiosis with atmospheric nitrogen-fixing bacteria (Benson, 1982; Rytter et al., 1989), alder trees play an essential role in forest soil nitrogen enrichment (Mander et al., 2008, 2021; Soosaar et al., 2011). Due to their rapid growth, alder species are frequently chosen for riparian buffer zones and short-rotation forestry (Aosaar et al., 2012; Rytter and Rytter, 2016; Uri et al., 2017).

However, ecosystem-level studies on carbon and water exchange in alder forests remain extremely limited. At the time of manuscript preparation, only two studies had reported ecosystem carbon exchange in grey alder forests. Uri et al. (2017) estimated net ecosystem production across an alder chronosequence using the carbon budgeting method. While informative, this traditional approach relies on discrete estimates of multiple carbon pools and fluxes and typically provides only an annual-scale assessment. In contrast, the eddy-covariance (EC) method provides continuous, high-frequency measurements of carbon and water fluxes between the ecosystem and the atmosphere, allowing for detecting intra-annual dynamics that can strongly influence the annual balance (Baldocchi, 2014). In our previous study (Krasnova et al., 2022), we conducted a two-year comparative analysis of EC carbon fluxes across four forested ecosystems, including the current site. That study focused on the effects of elevated air temperatures on carbon exchange and found that spring warming enhanced carbon uptake in the alder forest, indicating a positive temperature response during the early growing season. However, soil moisture variability, water fluxes, and post-drought recovery dynamics were beyond the scope of the analysis, leaving a substantial gap in our understanding of alder forest functioning.

In this study, we aim to investigate the ecosystem-level carbon and water exchange of a mature riparian alder forest stand in the hemiboreal zone in Estonia in relation to soil moisture variability. We utilise three years of EC flux measurements with contrasting hydroclimatic conditions: a “wet” year (2017), a “drought” year (2018), and a “recovery” year (2019). The specific objectives of this study are:

1. to quantify annual carbon and water exchange in an alder forest over three contrasting years;

2. to assess the influence of varying soil moisture conditions; and

3. to evaluate drought recovery and carry-over effects.

2 Methods

2.1 Study site and footprint area

The ecosystem in our study is a mature 40-year-old riparian grey alder (Alnus incana (L.) Moench) forest stand. It was established naturally following the abandonment of a wet meadow formerly used for haymaking and livestock grazing in southern Estonia. The terrain is flat, formed at the bottom of former periglacial lake systems, with an average elevation of 32 m a.s.l. and around 1 % inclination slope towards a tributary of the Kalli River. The average annual air temperature over the ten years before the study (2006–2016) was 6.6 °C with 627 mm yr⁻¹ of precipitation (Eesti Keskkonnaagentuur). The soil at the study site is classified as a Gleyic Luvisol, a hydromorphic soil type typical of seasonally waterlogged riparian zones. The humus layer thickness is 15–20 cm. The upper soil layer is moderately fertile, with a relatively high organic matter content, moderate total carbon and nitrogen concentrations, and a balanced C:N ratio (Table C1). Bulk density is relatively high, suggesting some compaction, likely due to past land use and seasonal wetting and drying cycles. Poor drainage and a fine-textured subsoil limit infiltration, making the site sensitive to both waterlogging and rapid topsoil drying during drought.

The footprint area of the tower (Fig. 1) is 1.65 ha, 85 % of which (1.41 ha) is covered by grey alder. The river, birch and spruce trees and a narrow section of the adjacent clear-cut represent the remaining area at the edges of the footprint. The average stand height is 17.5 m, the stand density is 1520 trees per ha, the mean stem diameter at breast height is 15.6 cm, and the basal area is 30.5 m² ha⁻¹ (Mander et al., 2022). The understory is dominated by herbs (Filipendula ulmaria (L.) Maxim., Aegopodium podagraria L., Cirsium oleraceum (L.) Scop., Geum rivale L., Crepis paludosa L., mosses (Climacium dendroides (Hedw.) F. Weber & D. Mohr, Plagiomnium spp. and Rhytidiadelphus triquetrus (Hedw.) Warnst. Moench), shrubs (Rubus idaeus L., Frangula alnus L., Daphne mezereum L.) and young trees (Alnus incana, Prunus padus L.).

Figure 1 Tower footprint area (10 %–80 %), Kljun model (Kljun et al., 2015); the blue line indicates the location of Kalli River. Map data: Estonian Land Board (Maa-amet), accessed via QGIS in 2023.

2.2 Instrumentation

The eddy-covariance (EC) setup consisted of a 3-D sonic anemometer Gill HS-50 (Gill Instruments Ltd., Lymington, Hampshire, UK) and enclosed CO₂ and H₂O gas-analyser LI-7200 (LI-COR Biosciences, Lincoln, NE, USA) measuring with 10 Hz frequency. The instruments were mounted on top of a 21 m scaffolding tower in spring 2017, with the first measurements starting on the 15 May 2017. Air temperature and relative humidity were measured using a Rotronic HC2A-S3 sensor (Rotronic AG, Bassersdorf, Switzerland), and shortwave radiation (Rg) was measured using a Kipp & Zonen CMP22 pyranometer (Kipp & Zonen B.V., Delft, The Netherlands). Twelve soil temperature (107, Campbell Scientific Inc., Logan, Utah, USA) and soil water content (ML3 ThetaProbe, Delta-T Devices, Burwell, Cambridge, UK) sensors were installed at 10 cm depth in the end of July 2017. Precipitation data were obtained from a nearby station, located around 2 km away from the site.

2.3 Fluxes calculation and post-processing

The fluxes of CO₂ and latent heat (LE) were calculated as a covariance between vertical wind speeds and CO₂ (or H₂O) concentrations using EddyPro software (version 6.3.0, LI-COR Biosciences, USA) and averaged over the 30 min intervals. In the absence of a storage measuring profile system, we estimated flux storage using the tower-top method, which utilised half-hourly CO₂ concentration measurements from the EC system. Net ecosystem exchange (NEE) was then calculated as the sum of eddy flux and storage. To eliminate periods of underdeveloped turbulence, we applied friction velocity filtering; the thresholds of 0.28 m s⁻¹ for 2017–2018 and 0.22 m s⁻¹ for 2019 were calculated with a moving point test (Papale et al., 2006). Fluxes during the half-hours with friction velocity values below these thresholds were removed from the analysis. To ensure adequate mixing conditions throughout the measurement period, we opted to remove not only nighttime half-hours, but also daytime NEE values associated with low friction velocity estimates.

In a previous study conducted at the same site by Mander et al. (2022), we noted that strong advection might be a feature of this forest site, with a rather dense canopy during the active vegetation period and a slight inclination towards the river tributary. To identify the periods when advection was significant, we applied the filtering method following Wharton et al. (2009) and Chi et al. (2019). Turbulence intensity parameters (I_w and I_u) were calculated for each half-hour as the ratios of vertical and horizontal wind velocity to turbulence intensity, respectively. For any half-hour, if I_w or I_u was outside of the window of mean plus one standard deviation estimated for the entire measurement period, advective conditions during this half-hour were considered non-negligible, and NEE and LE were filtered out. The remaining spikes in the dataset could be attributed to the simplification of the flux storage calculation procedure or the instrumental failure. Therefore, fluxes outside the common range (mean ±3 × standard deviation) were filtered out over a 14 d moving window (151 half-hour values). After all the filtering steps, 60 % in 2017, 66 % in 2018 and 65 % in 2019 of quality-controlled values remained.

Evapotranspiration (ET) was calculated by dividing the filtered LE by the latent heat of vaporisation (Allen et al., 1998). Energy balance closure (EBC) was 70 %, 71 % and 80 % in 2017, 2018 and 2019, respectively (Fig. A1). Given the simplified estimation of available energy in the absence of direct net radiation and ground heat flux measurements (see Appendix A for details), we chose not to apply any EBC correction to LE. This avoids introducing additional uncertainty and ensures that year-to-year comparisons of ET remain internally consistent.

In order to obtain fluxes aggregated over various time scales, we gap-filled NEE and ET using XGBoost as recommended by Vekuri et al. (2023). The hyperparameters were tuned during 5-fold cross-validation and included maximum tree depths (3, 5, 10, 15), regularisation strength with default 0, data sampling ratios (0.5, 0.75, 1), feature sampling ratios (0.4, 0.6, 0.8, 1), and minimum child weights (2, 5, 10). The hyperparameters were determined using all available data. A squared loss with a default learning rate of 0.1 was used as an objective function.

The partitioning of NEE into gross primary production (GPP) and ecosystem respiration (ER) was performed with the “nighttime” method in the ReddyProcWeb tool (Wutzler et al., 2018). Nighttime respiration was considered equal to nighttime gap-filled NEE values, while daytime ER was modelled in ReddyProcWeb using the air temperature dependence of measured nighttime values (Eq. 1)

$$\begin{array}{}\text{(1)}& \mathrm{ER}={\mathrm{ER}}_{\mathrm{ref}}{e}^{{\mathrm{E}}_{\mathrm{0}}\left(\frac{\mathrm{1}}{T_{\mathrm{ref}}-T_{\mathrm{0}}}-\frac{\mathrm{1}}{T-T_{\mathrm{0}}}\right)},\end{array}$$

where R_ref (µmol m⁻² s⁻¹) is the respiration at the reference temperature; E₀ (kJ mol⁻¹) is the activation energy; T (°C) is the measured air temperature. T_ref was set to 15 °C, and T₀ was kept constant at −46.02 °C following Lloyd and Taylor (1994).

GPP was then calculated as the difference between gap-filled NEE and modelled ER. We chose to use the nighttime flux partitioning method because, unlike the daytime method, where GPP is modelled, here GPP is derived indirectly as a residual. This approach allowed us to further calculate canopy physiological response parameters. Following the micrometeorological convention, negative flux denotes uptake, while positive flux is a release from the ecosystem into the atmosphere.

2.4 Canopy physiological response parameters

To study the physiological response of the ecosystem to varying soil moisture conditions, we calculated additional parameters: ecosystem water use efficiency (EWUE), canopy photosynthetic capacity (GPP_sat) and canopy conductance (G_c). Since these parameters characterise the vegetation activity, they were computed only for dry- (rainfall less than 1 mm d⁻¹) active-canopy (GPP < −1 g C m⁻² d⁻¹ and ET > 0.25 mm d⁻¹) days during the growing seasons. The start and end of each growing season were estimated by fitting a double-logistic curve to daily GPP sums and identifying the inflexion points, as outlined in Gonsamo et al. (2013). Canopy EWUE and G_c were calculated using only half-hours with sufficient light conditions. The threshold global radiation (Rg) value of 435 W m⁻² was computed from the bin-averaged GPP-Rg response curve in summer (JJA) dry- and active-canopy days of all three years using breakpoint analysis to identify the flattening point of the curve.

Ecosystem water use efficiency (EWUE), defined as the amount of carbon obtained by the forest per unit of water lost to the atmosphere, can serve as an indicator of a forest's adaptability to changing water availability (Huang et al., 2015; Keenan et al., 2013; Yang et al., 2016). We calculated EWUE as the ratio of the sum of absolute GPP values to the sum of ET, using two approaches. First, to compute annual and May–September EWUE, we used period sums of GPP and ET, including all data points. Second, to characterise canopy-specific EWUE, we calculated daily values focusing solely on periods of active photosynthesis and transpiration. For this, we included only half-hourly measurements taken under sufficient light conditions and restricted calculations to dry, active-canopy days within the growing seasons. Although we did not explicitly partition ET into evaporation and transpiration components, our filtering approach ensures that canopy-driven water fluxes dominate the ET.

Canopy photosynthetic capacity (GPP_sat) represents the forest's carbon uptake potential, i.e. how much carbon the ecosystem can sequester when light is not limiting (Aubinet et al., 2012; Chen et al., 2023; Fleischer et al., 2013; Musavi et al., 2017; Xu et al., 2020a). To obtain GPP_sat, we used a Michaelis-Menten photosynthetic response model (Michaelis and Menten, 1913), fitted to half-hourly daytime (Rg > 15 W m⁻²) GPP and global radiation data (Eq. 2). The fits were done in 3 d running windows, using only dry and active-canopy days, and parameters we assigned to the middle of each window.

$$\begin{array}{}\text{(2)}& \mathrm{GPP}={\displaystyle \frac{\mathit{\alpha }{\mathrm{GPP}}_{\max }{R}_{\mathrm{g}}}{\mathit{\alpha }{R}_{\mathrm{g}}+{\mathrm{GPP}}_{\max }}},\end{array}$$

where α (µmol J⁻¹) is the canopy light utilisation efficiency; GPP_max (µmol m⁻² s⁻¹) is the maximum GPP; Rg is global radiation (W m⁻²)

We chose the rectangular form of the light response curve over the more detailed non-rectangular one (Chen et al., 2023; Gilmanov et al., 2003; Musavi et al., 2017) because it demonstrated considerably better performance (a higher number of successful fits) with our dataset. However, a limitation of the simpler model is that the estimated GPP_max does not always correspond to the actual light saturation point. Therefore, we computed GPP_sat as GPP at Rg of 1000 W m⁻². Only the values from windows with significant fit parameters (p<0.05) and R²>0.5 were retained. For clarity in describing GPP_sat variability, we use its absolute values, omitting the negative sign that typically denotes flux direction.

Canopy conductance (G_c) is a representation of stomatal conductance on the ecosystem level. We computed G_c by the inversion of the Penman-Monteith equation (Eqs. 3, 4) only for the dry active canopy under sufficient light (Rg > 435 W m⁻²):

$$\begin{array}{}\text{(3)}& {\displaystyle }{G}_{\mathrm{c}}^{-\mathrm{1}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}{C}_p\left(\frac{\mathrm{VPD}}{\mathit{LE}}\right)}{\mathit{\gamma }}}+\left({\displaystyle \frac{\mathrm{\Delta }}{\mathit{\gamma }}}\mathit{\beta }-\mathrm{1}\right)g_{\mathrm{a}}^{-\mathrm{1}},\text{(4)}& {\displaystyle }{g}_{\mathrm{a}}^{-\mathrm{1}}={\displaystyle \frac{\mathit{\mu }}{{\mathit{\mu }}_{\ast }^{\mathrm{2}}}}+\mathrm{6.2}{\mathit{\mu }}_{\ast }^{-\mathrm{2}/\mathrm{3}},\end{array}$$

where ρ_a is the air density (kg m⁻³); C_p is the specific heat capacity of air (J kg K⁻¹); VPD is vapour pressure deficit (kPa); LE is latent heat flux (W m⁻²); γ is the psychrometric constant (kPa °C⁻¹), Δ is the slope of the saturation vapor pressure curve (kPa °C⁻¹); β is Bowen ratio; g_a is aerodynamic conductance (m s⁻¹); μ is wind speed (m s⁻¹), and μ_* is friction velocity (m s⁻¹).

To evaluate differences among years, we first detrended the daily data by subtracting the multi-year average seasonal cycle. The resulting anomalies were compared across years using the Kruskal-Wallis test, followed by pairwise Wilcoxon rank-sum tests with Bonferroni correction for multiple comparisons.

2.5 The impact of varying soil moisture conditions

To evaluate the influence of soil moisture conditions on ecosystem functioning, we computed the soil saturation ratio (SSR) as the ratio of measured soil water content (SWC) to its 99th percentile across the entire observation period. We then analysed the response of carbon and water fluxes to SSR variability by binning the data into SSR intervals while controlling for the main environmental drivers that could otherwise overshadow the effects of soil moisture variability. Because of the strong seasonality of fluxes and the lack of LAI data, we restricted the analysis to summer months (JJA) to ensure a fully developed canopy.

Canopy photosynthetic capacity (GPP_sat) was calculated for each 0.1 SSR bin using daytime data from dry- and active-canopy days. This was done both for data pooled across all three summers, to capture the overall response pattern, and separately for each summer to assess interannual variability. Using GPP_sat rather than GPP allowed the removal of light as the primary driving factor.

To control for temperature effects on ER, we initially attempted to obtain ER at fixed temperature by fitting temperature response curves to measured nighttime ER data (nighttime NEE) in SSR bins. However, the range of air temperatures within bins was insufficient. Instead, we used reference ER at 15 °C (ER_ref) as provided by the ReddyProcWeb partitioning procedure, which allowed the impact of SSR on ER to be analysed independently of air temperature variability.

Because ET is strongly driven by VPD, we divided ET by VPD and calculated mean values for 0.05 SSR bins. To minimise the influence of evaporation, only half-hourly data from dry- and active-canopy days and under sufficient light conditions were included. Similarly, EWUE was normalised by daytime VPD before calculating averages for 0.05 SSR bins.

To examine how soil moisture modifies the sensitivity of canopy conductance, G_c was divided into 0.1 SSR bins, using data from all three summers. For each SSR bin, reference canopy conductance (G_c,ref) was estimated by fitting the Oren et al. (1999) model (Eq. 5):

$$\begin{array}{}\text{(5)}& G_{\mathrm{c}}=-m\ln \left(\mathrm{VPD}\right)+G_{\mathrm{c},\mathrm{ref}},\end{array}$$

where G_c is canopy conductance (mm s⁻¹), m is the stomatal sensitivity, G_c,ref is reference canopy conductance at 1 kPa, VPD is vapour pressure deficit (kPa).

Although the slope ($m/G_{\mathrm{c},\mathrm{ref}}$) generally fell within the expected range of 0.5–0.7, three SSR classes with low R² values (0.05, 0.03 and 0.15) exhibited notably lower slopes (0.23, 0.23 and 0.44, respectively; Table B1, Fig. B1). To evaluate whether G_c,ref estimates were biased by poor model fits, we derived an additional set of G_c,ref with fixed $m/G_{\mathrm{c},\mathrm{ref}}=\mathrm{0.6}$ (Fig. B2). This approach improved R² values, while G_c,ref estimates remained largely unchanged (Table B1). We therefore based subsequent analyses on G_c,ref values calculated with the fixed slope, while also indicating the alternative estimates in Fig. 6b. Similar analysis was carried out for each growing season separately to assess the interannual difference in G_c sensitivity to soil moisture variability.

2.6 Drought recovery and carry-over effects

To disentangle the possible carry-over effects of the drought year from the natural interannual variability, we applied a two-step approach combining model-based analysis of ecosystem fluxes with resistance, recovery and resilience indices. To ensure temporal consistency, we restricted the analysis to a common portion of the growing season (May–September) for each year.

First, we assessed whether observed interannual variability in GPP and ER could be attributed solely to changes in their primary environmental drivers (light and temperature, respectively). To do this, we estimated Michaelis-Menten light response curve parameters (Eq. 2) within a running three-day window using half-hourly daytime GPP and Rg data for each year separately. For ER, we utilised ER_ref and E₀ parameters derived during flux partitioning in the ReddyProcWeb tool (Eq. 1). Each year's parameter set was then used to model GPP and ER across all three years using measured Rg and air temperature of each corresponding year. This cross-year modelling allowed us to test whether model parameters obtained for one year could accurately predict flux dynamics in other years. Differences between fluxes when applying parameters from a non-drought year to drought or recovery years (and vice versa) can thus highlight a possible carry-over effect, reflecting changes in ecosystem functioning that persist beyond immediate environmental conditions. We used the Kruskal-Wallis test with Bonferroni correction to assess whether differences between observed and cross-year modelled fluxes were statistically significant.

To further quantify the magnitude of the drought impact and evaluate the ecosystem's ability to recover, we calculated resistance (Rt), recovery (Rc), and resilience (Rs) indices for daily ecosystem parameters (Lloret et al., 2011; Portela et al., 2023). Performing the analysis on daily values allowed us to capture short-term fluctuations that could be partially masked by aggregating the parameters seasonally or annually. Resistance was calculated as the ratio of drought to pre-drought year, quantifying the immediate drought impact (Eq. 6). Recovery was computed as the ratio of recovery to drought year. Resilience was calculated as the ratio of recovery to the reference year, quantifying the ability of the ecosystem to return to pre-drought levels and allowing identification of potential carry-over effects.

$$\begin{array}{}\text{(6)}& \mathrm{Rt}={\displaystyle \frac{\mathrm{dry}}{\mathrm{ref}}};\mathrm{Rc}={\displaystyle \frac{\mathrm{rec}}{\mathrm{dry}}};\mathrm{Rs}={\displaystyle \frac{\mathrm{rec}}{\mathrm{ref}}},\end{array}$$

The ecosystem parameters included daily values of carbon exchange components (GPP and ER) and their main driver-normalised versions (GPP_sat and ER_ref), ET, transpiration (T), EWUE and G_c. Transpiration was estimated as the daily average of ET under sufficient light (Rg > 435 W m⁻²) for dry (rainfall less than 1 mm d⁻¹) and active-canopy (GPP < −1 g C m⁻² d⁻¹ and ET > 0.25 mm d⁻¹) days. To account for the strong day-to-day variability of these parameters, which could significantly bias average estimates, we applied non-parametric bootstrapping (n=5000), resampling daily values within each year independently. For each index, we report the bootstrapped mean and 95 % confidence intervals.

We acknowledge a limitation in the selection of a single reference year, 2017, which may not fully represent long-term baseline conditions. Consequently, both the magnitude and interpretation of the indices should be viewed in the context of this wet year reference. Additionally, we note that the interpretation of drought-induced changes in EWUE differs from that of other variables. While increases in EWUE may suggest that the ecosystem is coping under stress by maintaining carbon uptake relative to water loss, they often result from stomatal regulation and reduced transpiration, and thus may reflect a physiological stress response rather than enhanced functioning.

3 Results

3.1 Weather conditions

Meteorological conditions across the three study years exhibited a clear seasonal pattern typical for the region (Fig. 2). The mean annual air temperature in 2017 was close to the 10-year average (6.6 °C), while both 2018 and 2019 were around 1 °C warmer (Table 1). Based on differences in precipitation (P) and soil water content (SWC), we categorised the years as “wet” (2017), “drought” (2018), and “recovery” (2019), reflecting conditions before, during and after the drought. Although in-situ SWC measurements began only in mid-July 2017, observations from a nearby station (Fig. C1) and visual assessments during instrumentation setup confirmed elevated SWC levels (including standing water) also in late spring and early summer of that year.

Figure 2 Hydroclimatic conditions during the three studied years (a–d): 10 d running means of meteorological parameters; dashed lines denote the beginning and end of the corresponding growing seasons. Daily soil water content and cumulative daily precipitation (e) of May–September. x-axis ticks correspond to the beginning of each month.

Compared to the 12-year mean of 336 ± 75 mm, total May–September precipitation was 11 % higher in the “wet” year, 19 % lower in the “drought” year and similar in the “recovery” year, with all three years falling within the one standard deviation of the long-term mean (Table 1). However, the temporal distribution of the rainfall varied among the years (Fig. 2e).

In 2017, rainfall was skewed toward the second part of the season (August–September), resulting in elevated SWC (up to 0.64 m³ m⁻³) and localised flooding. In contrast, 2018 experienced extended dry spells in May and July, resulting in a pronounced soil moisture deficit (Fig. 2e). SWC declined from 0.73 to 0.28 m³ m⁻³ over 33 d in May (1 May–3 June, −0.014 m³ m⁻³ d⁻¹) and from 0.33 to 0.11 m³ m⁻³ over 32 d in July (2 July–3 August, −0.007 m³ m⁻³ d⁻¹). These periods of progressive drought coincided with elevated vapour pressure deficit (VPD), and the second part of July to early August period overlapped with high air temperatures, previously identified as a heatwave period (Krasnova et al., 2022). Total precipitation in the active season of 2018 was 27 % lower than in 2017.

In 2019, precipitation was more evenly distributed throughout the season, with no extended dry spells, resulting in intermediate SWC levels and cumulative precipitation 25 % higher than in 2018 and 10 % lower than in 2017 (Table 1). Although our SWC measurements were limited to the upper soil layer (10 cm depth), they remain representative of hydrological conditions affecting alder root activity, as roots are predominantly confined to shallow depths due to adaptation to waterlogged, compacted soils.

3.2 Interannual differences in accumulated fluxes

In each of the three study years, the alder forest acted as a net carbon sink (cumulative annual NEE < 0; Table 1). Cumulative NEE in the active vegetation season (May–September) accounted for 96 % of the total annual flux (97 % in 2018 and 95 % in 2019). Based on this seasonal share, we estimated the total annual NEE for 2017 to be −600 g C m⁻² yr⁻¹, reflecting a smaller net CO₂ uptake than in 2018, but higher than in 2019. Following the same approach, we obtained estimates for annual GPP in 2017 (May–September GPP also accounting for an average of 96 % of the total), and subsequently calculated ER as the difference between GPP and NEE. Similarly, as the majority of evapotranspiration (ET) occurred during the active season (94 % in 2018 and 91 % in 2019), total annual ET and ecosystem water use efficiency (EWUE; GPP/ET) were estimated for 2017. Across the study period, the three-year means (±SD) for annual NEE, GPP, ER, and ET were −586 ± 85, −1330 ± 82, 743 ± 166 g C m⁻² yr⁻¹ and 264 ± 74 mm yr⁻¹, respectively.

Table 1 Annual and May–September average air temperature (Ta), soil temperature at 10 cm depth (Ts), global radiation (Rg), and soil water content at 10 cm depth (SWC); cumulative precipitation (P), net ecosystem exchange (NEE), gross primary production (GPP), ecosystem respiration (ER) and evapotranspiration (ET), ecosystem water use efficiency (EWUE) in the wet (2017), drought (2018) and recovery (2019) years. May–September values are shown in parentheses. Annual fluxes for 2017 (in italics) were estimated from the ratio of May–September values to the full year of measurements (see text for details). EWUE was calculated from the sums of GPP and ET for the corresponding periods.

^a data from the Estonian National Weather Service; ^b starting from 24 July 2017.

The net carbon uptake of the drought year was highest among the three years, with May–September NEE exhibiting a 35 % increase relative to the wet year (p<0.001). This enhanced net sink was a result of a significant 36 % reduction in ER (p<0.001), while the 3.4 % decline in cumulative GPP was not significant. In the recovery year, NEE during the active season was 18 % lower than in the drought year (p<0.001) yet remained similar to the wet year (+10 %, p=0.02). ER featured the largest interannual change, significantly increasing by 62.5 % (p<0.0001) in 2019 relative to 2018, but only 4.7 % above 2017, with no significant difference in daily fluxes during the active season (p>0.1). GPP increased moderately in the recovery year (+11 % compared to 2018, p=0.0004), while remaining statistically similar to the wet year (+7.4 %, p>0.1).

During the drought year, total evapotranspiration (ET) for the active season decreased significantly (p<0.001) by 23 % relative to the wet year. ET in the recovery year was the highest of the three years, exceeding the drought and wet years by 71 % and 32 %, respectively. Interannual differences in daily ET were significant across all years (p<0.001). The active season $\mathrm{ET}/P$ ratio was similar in wet and drought years (0.69 and 0.67, respectively) but increased in the recovery year (0.92). EWUE peaked in the drought year, reflecting stable carbon uptake under reduced water loss, and was lowest in the recovery year, both annually and seasonally (Table 1).

3.3 Seasonal dynamics of carbon and water fluxes, and canopy physiological response parameters

Carbon exchange components (NEE, GPP, ER) and ET exhibited clear seasonal patterns in all years, with a sharp increase in spring, peak rates around mid-summer, and a decline toward autumn, reflecting the phenological cycle of a deciduous forest in the hemiboreal zone (Fig. 2). However, the timing, magnitude and duration of flux peaks varied among the study years.

Figure 3 Seasonal dynamics of net ecosystem exchange (NEE, a), ecosystem respiration (ER, b), gross primary production (GPP, c), evapotranspiration (ET, d) represented by 10 d running means. Vertical dashed lines are the borders of growing seasons. x-axis ticks correspond to the beginning of each month.

In 2017, lower spring and summer temperatures (Fig. 2) delayed the onset of ER, GPP and ET (Fig. 3b–d), all of which increased at a slower rate compared to 2018 and 2019. High SWC supported higher than in 2018 ET throughout the season, but lower VPD limited evaporative demand, resulting in ET rates below those in the recovery year (2019). All three fluxes peaked in August, following a period of warmer temperatures and clear-sky conditions, with ET exceeding the recovery year levels. Sustained ER during September-October contributed to an earlier decline in NEE compared to 2018 and 2019 (Fig. 3a).

In contrast, the 2018 drought year was characterised by higher spring and summer air and soil temperatures, elevated VPD, and a progressive decline in SWC, particularly in May and again from July onward. These conditions contributed to an earlier rise in ER, GPP and ET, followed by a sharp suppression of all fluxes once SWC became limiting. Despite high atmospheric demand, ET noticeably declined in July and August, consistent with water limitation. ER declined more than GPP in late summer, resulting in a more negative NEE, reflecting increased net carbon uptake toward the end of the growing season.

In the recovery year 2019, spring conditions resembled those of 2018, but lower air temperatures and a more uniform precipitation distribution prevented soil moisture depletion. Warm summer with peaking air temperatures in June supported earlier peaks in ER and GPP, with cooler and wetter conditions in July and August co-occurring with a moderate decline in both fluxes. ET exhibited two distinct peaks, in June and again in August–September, and remained almost consistently higher than in the other two years, supported by moderate VPD, ample SWC, and frequent rainfall events.

The seasonal dynamics of ER and GPP are primarily governed by their main environmental drivers, temperature and light, respectively. Moreover, daytime ER is modelled based on temperature, which can bias direct comparisons across years. Therefore, ER_ref (ER at reference temperature) and GPP_sat (canopy photosynthetic capacity, i.e. GPP at saturating light) are more objective measures for evaluating interannual differences in seasonal variability (Fig. 4).

Figure 4 Seasonal dynamics of reference ecosystem respiration (ER_ref, a) and canopy physiological response parameters: canopy photosynthetic capacity (GPP_sat, b), water use efficiency (EWUE, c), and canopy conductance (G_c, c). Lines denote 10 d running means. x-axis ticks correspond to the beginning of each month. EWUE markers are daily values that were calculated from sums of GPP and ET filtered for the active photosynthesis and transpiration under sufficient light (Rg > 435 W m⁻²). GPP_sat makers are the centres of 3 d running windows. G_c markers are the daily average values from the half-hourly estimates under sufficient light. All canopy parameters were obtained only for the days with <1 mm of rain, GPP > 1 g C m⁻² d⁻¹ and ET > 0.25 mm d⁻¹ for three growing seasons.

ER_ref exhibited clear seasonal patterns across three years, with seasonal maxima in August 2017, May 2018, and June 2019 (Fig. 4a). The severe suppression of ER_ref from June to October 2018 confirms that factors beyond temperature strongly influenced ecosystem respiration during the drought year. In contrast, ER_ref remained relatively stable throughout active season in 2017 and 2019. Average ER_ref were 4.7 ± 1.3, 2.5 ± 1.2 and 4.3 ± 1.3 µmol m⁻² s⁻¹ in 2017–2019, respectively, with no significant difference detected between the wet and recovery years (p>0.1).

GPP_sat peaked in June in all three years (Fig. 4b), followed by a sharp mid-summer decline coinciding with VPD peaks (Fig. 2c). In May and early June 2018, it was higher compared to the other years, but the late-summer rebound observed in 2017 and 2019 was absent. Despite these seasonal differences, the average values were similar 24.9 ± 9.1, 22.4 ± 7.5 and 23.5 ± 8.2 µmol m⁻² s⁻¹ in 2017, 2018, and 2019 (p>0.1). Canopy EWUE varied between the years (Fig. 4c), with the drought year demonstrating higher and more variable values (5.2 ± 2.9 g C kg H₂O⁻¹) than both 2017 (4.1 ± 1.0 g C kg H₂O⁻¹) and 2019 (3.2 ± 0.7 g C kg H₂O⁻¹) and peaking in July and September 2018. The difference between all three years was significant (p<0.001).

Canopy conductance (G_c) followed a similar seasonal pattern in the wet and recovery years, with higher values during summer relative to early and late in the growing season (Fig. 4c). Nevertheless, the difference between 2017 and 2019 was significant (p<0.001), with the recovery year exhibiting a lower average G_c (9.5 ± 4.4 mm s⁻¹) compared to the wet 2017 year (11.1 ± 4.9 mm s⁻¹). In contrast, G_c remained consistently low throughout the active season of the drought year, averaging 4.6 ± 1.4 mm s⁻¹, significantly lower than both other years (p<0.001).

3.4 The impact of soil moisture variability

Over the three contrasting summers, only the data from the recovery year covered the majority of soil saturation ratio classes (SSR, SWC normalised to the 99th percentile), with the lowest SSR class occurring only in the drought year (Fig. 5). GPP_sat demonstrated an optimum at 40 %–50 % saturation, declining at both lower and higher SSR. ER_ref increased from low to medium SSR and remained relatively stable with further increase in SSR. However, in 2017, ER_ref remained higher at high SSR than in the recovery year.

Figure 5 The variation of canopy photosynthetic capacity (GPPsat, a), reference ecosystem respiration (ERref, b), normalised canopy evapotranspiration ($\mathrm{ET}/\mathrm{VPD}$, c) and normalised water use efficiency (EWUE × VPD, d) under different soil saturation ratios (SSR) in summer months (JJA) of three study years. GPPsat (a) is shown as the value estimated for each SSR class and year; shaded areas represent the 95 % confidence intervals. Other parameters (panels b–d) are shown as averages within 0.05 SSR bins; shaded areas denote ± 1 standard deviation. Dashed lines in all panels indicate results for all data pooled across years. For details on data filtering and calculations, see the corresponding section of the Methods.

Because both ET and EWUE are strongly influenced by VPD, we normalised them to isolate soil moisture effects (Fig. 5c–d). ET increased with SSR, reaching maximum values near 50 % before stabilising, similar to the pattern observed in ER_ref. In contrast, VPD-normalised EWUE declined with increasing SSR, but with distinct interannual differences. In 2017, EWUE remained low and relatively constant at moderate saturation (∼ 40 %–60 %), decreasing further only from 70 % saturation. In 2018, EWUE was higher across most SSR classes compared to other years, with only a minor decline from 10 % to 80 % and a drop near full saturation. The recovery year (2019) showed a unimodal response, with low EWUE at 20 %–50 %, a moderate rise at 60 %, and a decline at higher SSR.

To examine stomatal regulation under varying soil moisture conditions, we assessed the response of G_c to VPD across the SSR classes, both for all three growing seasons combined (Fig. 6a–b) and for each year separately (Fig. 6c). Reference G_c (G_c,ref, canopy conductance at 1 kPa VPD) was the highest at moderate saturation (45 %) and the lowest at the driest (15 %) and wettest (85 %) conditions when data from all years were pooled. This unimodal response was also evident in the wet and recovery years analysed separately. In contrast, during the drought year, G_c,ref remained low across all SSR classes, with only a small non-significant (overlapping CIs) increase at 35 %–55 %. Overall, G_c,ref in 2018 was about half the magnitude observed in the wet and recovery years at low to medium SSR, approaching recovery-year values only at SSR ≥65 %, when the G_c,ref of the other two years was suppressed.

Figure 6 (a) Sensitivity of canopy conductance (G_c) to vapour pressure deficit (VPD) across different soil saturation ratio (SSR) classes, indicated by shades of blue, with darker colours representing higher SSR. Curves represent the Oren et al. (1999) model (Eq. 5) with a fixed slope of 0.6. (b) Reference canopy conductance (G_c,ref) for each SSR class, where SSR values denote the midpoint of each 0.1 interval. Error bars represent 95 % confidence intervals. Squares indicate G_c,ref estimates derived from the variable-slope model (see Appendix B for more details). (c) Estimated G_c,ref per SSR class for each study year. Bars with lighter shading indicate lower quality model fits (R²<0.2). Error bars represent 95 % confidence intervals.

Figure 7 The change in canopy physiological response parameters (photosynthetic capacity, GPP_sat; canopy conductance, G_c; water use efficiency, EWUE) with the progressive soil drying in 2018. Shaded bars are the daily soil saturation ratio (SSR, blue) and the daily air temperature (red). All other values are normalised 3 d running means.

In May and July 2018, the forest under study experienced a progressive drought, with soil moisture declining in the absence of rainfall over multiple consecutive days (Fig. 7). SSR decreased from full saturation in early May to 37 % by the end of the month, remained between 30 %–50 % in June, and dropped further to 15 % over July, persisting at low levels through mid-August (Fig. 7). While the soil was drying in May, all canopy physiological response variables increased, with maximum values reached under moderate SSRs in June. During the progressive drought in July, most variables declined, with the exception of VPD-normalised EWUE, which reached its highest values under the driest conditions. In early August, when SSR remained at its minimum, GPP_sat stabilised at approximately half of its June peak, while VPD-normalised EWUE fell from 75 % to 45 % of its maximum.

3.5 Drought recovery

To assess how the ecosystem recovered from drought, we compared daily GPP and ER of each year with cross-year modelled values from the other two years. The difference between observed and cross-year modelled daily GPP (Fig. 8a–c) was significant across all years (p<0.001), although the values of 2017 and 2019 GPP closely matched in August–September. In 2018, the suppression of GPP starting mid-June was also present in cross-year modelled values (Fig. 8a, c; dashed orange lines). Observed and cross-year modelled ER in 2017 and 2019 (Fig. 8d–f) did not differ significantly (p>0.1), while ER in 2018 was lower than cross-year modelled values (p<0.001).

Figure 8 Gross primary production (GPP; panels a–c) and ecosystem respiration (ER; panels d–f) during the growing seasons of 2017 (a, d), 2018 (b, e), and 2019 (c, f). Solid lines are observed GPP (upper panels) and ER (lower panels) for each year, while dashed lines are cross-year modelled fluxes calculated using the measured global radiation (for GPP) or air temperature (for ER) of each year with parameter sets derived from other study years. All lines represent 10 d running means.

To further quantify the drought impact and recovery of the daily ecosystem parameters (GPP, GPP_sat, ER, ER_ref, ET, T, EWUE, G_c), we estimated resistance (Rt), recovery (Rt) and resilience (Rs) indices (Fig. 9). According to overlapping 95 % confidence intervals, no significant differences were detected between carbon fluxes and their driver-normalised versions (GPP vs. GPP_sat; ER vs. ER_ref) for any of the indices.

The ecosystem's resistance (Rt) to drought exhibited considerable variation among the studied parameters. Carbon uptake (GPP and GPP_sat) maintained moderate resistance (0.84 and 0.88, respectively), indicating partial suppression of daily photosynthetic activity in 2018. In contrast, respiratory fluxes (ER and ER_ref) exhibited much lower resistance (0.57 and 0.53, respectively). Water fluxes were likewise reduced, with Rt = 0.63 for evapotranspiration (ETtot) and Rt = 0.71 for transpiration (T). EWUE exhibited the highest resistance (∼ 1.38), whereas G_c was most affected, with Rt = 0.39.

All ecosystem parameters demonstrated full recovery (Rc > 1), with EWUE declining to values below the drought level (Rc = 0.58). Resilience (Rs) indices indicated a generally strong rebound in the year following drought (2019), with carbon and water fluxes returning to or exceeding pre-drought levels observed in 2017, and EWUE declining to even lower levels (Rs = 0.79), indicating a return to less constrained conditions. However, G_c exhibited lower resilience (Rs = 0.88), which may modulate ecosystem responses to future stress events.

Figure 9 Resistance (Rt), recovery (Rc) and resilience (Rs) indices of the gross primary production (GPP), canopy photosynthetic capacity (GPPsat), ecosystem respiration (ER), reference ecosystem respiration (ERref), evapotranspiration (ETtot), transpiration (T, calculated as filtered ET with maximum share of T), canopy water use efficiency (EWUE) and canopy conductance (G_c). Error bars are 95 % confidence intervals.

4 Discussion

4.1 Alder forest as a strong net carbon sink

Alder is a widely distributed tree species across hemiboreal and temperate zones, commonly found in riparian buffers, yet data on alder forest ecosystem carbon exchange remain surprisingly limited. In a chronosequence of alder forest stands in Estonia studied by Uri et al. (2017), the two oldest sites (“Agali” and “Kolleste 2”) were of similar age to the current one (∼ 35 and ∼ 40 years old, respectively) at the time of measurements. Uri et al. (2017) applied a biometric methodology, which integrates stand biomass, production, litterfall, and monthly soil respiration measurements to estimate carbon exchange components. Although this approach differs from the eddy-covariance method used here and the results are not directly comparable, it is possible to approximate NEE, GPP, and ER from their data under a set of assumptions (see Appendix D for details). The calculated values of NEE, GPP, ER, and total soil respiration (Rs) from Uri et al. (2017), using two-year averages for Agali and one-year data for Kolleste 2, are summarised in Table 2 together with the sites' characteristics.

Table 2 Site characteristics and carbon exchange parameters from the current study and from two mature alder forests reported in the chronosequence study by Uri et al. (2017). NEE, GPP, and ER for “Agali” and “Kolleste 2” were calculated using Rs data from Table 8 and Rh, NEP, and NPP data from Table 9 in that study, following the method described in the Appendix D.

While our site acted as a strong net carbon sink, the two sites from Uri et al. (2017) were a much weaker sink (Agali) and a net carbon source (Kolleste 2). The GPP of our site was intermediate between the two, yet ecosystem respiration was 1.5 to 1.7 times lower. These differences can partly stem from contrasting soil properties: the current site's Gleyic Luvisol exhibited much higher bulk density (1.7 g cm⁻³), moderate acidity (pH 5.3), and lower soil organic matter content (6.5 %) compared to the other sites, which had lighter soils, higher SOM, and, in the case of Kolleste 2, strongly acidic conditions (pH 3.7). Higher bulk density likely limits soil aeration and microbial activity, reducing respiration rates and favouring net carbon uptake, while more acidic soils and higher organic matter at Kolleste 2 may promote microbial respiration, resulting in enhanced soil respiration and, consequently, net carbon release. However, it should be noted that the study periods differ; Uri et al. (2017) conducted measurements between 2011 and 2014, whereas the current study covers 2017 to 2019, and interannual weather variability during these periods may have contributed to observed differences in carbon fluxes.

In the same years of measurement, a mature upland pine forest growing on sandy soil in Estonia exhibited similar GPP but higher ER, resulting in a weaker net carbon sink (Table 3). Compared to previously reported values for various broadleaved forests in boreal and hemiboreal zones, the NEE at our site exceeds most estimates but aligns with fluxes observed in more southern broadleaved and coniferous forests. While GPP at our site was comparable to that of boreal and hemiboreal forests, ER was, again, notably lower. In contrast, forests with a similar NEE range exhibited higher GPP but also greater ER, likely driven by their warmer climate with a longer active vegetation season.

Very low ER in our study likely reflects oxygen limitation in a compact, frequently wet mineral soil, rather than nutrient shortage. The slightly acidic soil with high bulk density forms conditions that reduce gas diffusion and favour anoxic microsites, suppressing microbial decomposition despite moderate total C. Together with restricted fine-root activity under dense, saturated conditions, these factors could lead to low ecosystem respiration. We note, however, that our observations are limited to three years with exceptional weather conditions, which may not fully capture the “typical” respiration rates of this forest.

Table 3 Net ecosystem exchange (NEE), gross primary production (GPP), and ecosystem respiration (ER) reported for broadleaf forests and sites with values comparable to this study. “(s. yrs.)” indicates that averages were calculated using only the same years as in this study (2017–2019). All values are in g C m⁻² yr⁻¹, mean ±SD where available.

Anoxic conditions, combined with fluctuating soil moisture levels, are very favourable for methane (CH₄) production (Dlugokencky et al., 2011; Feng et al., 2020; Flanagan et al., 2021). In addition, the high nitrogen content typical of alder forests could promote nitrous oxide (N₂O) emissions (Butterbach-Bahl et al., 2013; Davidson et al., 2000). In principle, these non-CO₂ fluxes could offset the strong net carbon sink observed in our study. However, chamber-based and eddy-covariance measurements at the site (Mander et al., 2021, 2022) indicate otherwise. On an annual scale, the alder stand functioned as a very minor CH4 sink, contributing merely 0.1 % to the total GHG global warming potential (Table E1). While annual N₂O fluxes were positive, they represented just 1.1 % of total NEE in CO₂-equivalent units (or 7.9 % based on chamber-derived estimates), too small to negate the forest's role as a substantial net CO₂ sink.

4.2 Lower than expected evapotranspiration

Evapotranspiration (ET) integrates physical evaporation from soil and wet surfaces with plant transpiration, making it inherently complex (Brümmer et al., 2012; Jarvis, 1986; Jassal et al., 2009; Massmann et al., 2019) and difficult to quantify accurately (Fisher et al., 2017). Eddy covariance estimates are further affected by incomplete energy balance closure (Appendix A) (Amiro, 2009; Foken, 2008; Mauder et al., 2018, 2020). We chose to report the ET based on directly measured LE, as net radiation and ground heat flux were not available for accurate adjustment (Mauder et al., 2018).

Annual ET in the wet year 2017 (255 mm yr⁻¹) and especially in the drought year 2018 (194 mm yr⁻¹) were lower than expected for the boreal and hemiboreal region (Launiainen et al., 2022; Lindroth et al., 2020; Wang et al., 2021) and much lower than in the various riparian forests (Kochendorfer et al., 2011; Kowalska et al., 2020; Scott et al., 2004; Xu et al., 2020b). However, ET in the recovery year 2019 (342 mm yr⁻¹) was closer to previously reported annual ET values in boreal and hemiboreal forests in 2015–2018 (Lindroth et al., 2020) and multi-year averages reported by Wang et al. (2021) (384 ± 12 mm yr⁻¹) and Launiainen et al. (2022) (348 ± 26 mm yr⁻¹). The average precipitation during the growing season in the latter study (383 ± 83 mm) was comparable to our wet (372 mm) and recovery (338 mm) years but exceeded the drought year (271 mm), indicating that lower atmospheric water supply likely contributed to the reduced ET in 2018. Lower ET under higher precipitation in the wet year may be explained by cooler summer temperatures reducing evaporative demand. On the other hand, lower energy balance closure levels in 2017 and 2018 (70 % and 71 %, respectively) might contribute to the ET underestimation (Fig. A1). Furthermore, the low ET with moderate GPP results in notably higher annual EWUE (5.3 ± 1.2 g C kg H₂O⁻¹), which exceeded values (0.9–4.1 g C kg H₂O⁻¹), previously reported for various forests (Jin et al., 2023; Niu and Liu, 2021; Xie et al., 2016; Zhou et al., 2014).

The evaporative index ($\mathrm{ET}/P$) in 2017 and 2018 (0.69 and 0.67, respectively) was consistent with values reported for other forested ecosystems (Williams et al., 2012), whereas a higher ratio of 0.92 in 2019 likely reflects the combined effects of a warm growing season and well-distributed precipitation, which stimulated both photosynthesis (Table 1) and transpiration (Fig. 9). As previously demonstrated by Eschenbach and Kappen (1999), alder's high leaf stomatal conductance supports enhanced transpiration under sufficient water supply, implying that conditions in 2019 may have been near optimal for maximising water and carbon exchange. We have to note that the evaporative index remained below one in all years of our study, which is surprising for a riparian forest that typically has access to additional water through lateral inputs, and thus ET would exceed P (Kochendorfer et al., 2011; Xu et al., 2020b). The shallow and narrow forest river near our study site likely provided only limited supplemental water. While this forest is experiencing seasonal flooding, it might be a result of historically formed high-density soils, rather than the river's impact in the three studied years.

4.3 Moderate soil water saturation enhances ecosystem fluxes

Soil moisture variability plays an important role in modulating ecosystem carbon exchange, although its effect is usually more pronounced in water-limited, rather than radiation-limited regions (Green et al., 2019; Kannenberg et al., 2024). Nevertheless, variable soil moisture conditions are intrinsic to riparian forests, where seasonal flooding and fluctuating groundwater table create a dynamic hydrological regime (Kowalska et al., 2020; Portela et al., 2023; Singer et al., 2014). On our site, the variation in soil saturation (SSR) levels caused a non-linear response of canopy gas exchange, with both very dry and very wet conditions constraining reference stomatal conductance (G_c,ref, G_c at 1 kPa of VPD), representing the maximum stomata aperture (Fig. 6). Soil moisture extremes impose both hydraulic and metabolic constraints: low saturation limits water supply to leaves, while high saturation can cause oxygen limitation in the rhizosphere, impairing root function and nutrient uptake (Kochendorfer et al., 2011; Kozlowski, 1997). The physiological optimum for canopy conductance at a moderate soil saturation (∼ 40 %–50 %) was reflected in photosynthetic capacity and transpiration, and enhanced carbon and water fluxes in the corresponding growing season periods of different years.

Relatively high values of photosynthetic capacity and ET, sustained into wetter ranges during summer months, especially during the colder year, may reflect adaptations to periodic flooding typical of riparian forests. Canopy EWUE declined with increasing soil water content in the wet year, consistent with a “relaxed” physiological state when water is not limiting. Reference respiration (ER_ref) also followed a similar saturation curve, with moderate soil moisture promoting optimal metabolic activity, while substrate constraints under lower saturation ratios reduced respiration rates.

4.4 Alder forest in the 2018 drought year

In 2018, low precipitation caused widespread soil moisture deficits across Europe, while extremely high air temperatures further intensified drought conditions through elevated VPD (Fu et al., 2020; Lindroth et al., 2020; Smith et al., 2020). Estonia also experienced two extended dry spells (in May and July), accompanied by an exceptional heat anomaly from mid-July to early August (Krasnova et al., 2022). These conditions produced a progressive summer drought at our site (Fig. 7).

In May, while the soil drying was faster than in July, it reached only ∼ 50 % soil saturation, which persisted through June. We found this moisture level to be optimum for ecosystem fluxes (Fig. 5) and plant stomatal activity (Fig. 6), although still constrained by the early stages of the growing season. The relatively low canopy conductance at that time likely reflected ongoing alder leaf development. Under these favourable early-season water conditions and warmer-than-average May temperatures, all fluxes and canopy physiological parameters gradually increased (Fig. 7).

Higher spring temperatures can enhance annual net carbon uptake by extending the growing season (Keenan et al., 2014; Wolf et al., 2013) and offsetting the influence of the forthcoming summer drought on the annual carbon balance (Angert et al., 2005; Kljun et al., 2006; Smith et al., 2020; Wolf et al., 2016). A similar pattern was reported for a riparian mixed broadleaf forest in the Czech Republic, where an anomalously warm spring in 2018 led to an increase in both GPP and ET, counteracting the negative effects of the summer drought (Kowalska et al., 2020). In boreal and hemiboreal regions, moderate spring warming typically coincides with ample soil moisture from snowmelt, ensuring sufficient water supply for early-season assimilation. However, enhanced spring productivity and transpiration can also accelerate soil water depletion, increasing susceptibility to summer drought (Bastos et al., 2020).

In July 2018, the soil moisture decreased further, reaching its minimum by August. Combined with extremely high temperatures and VPD peaks, the progressive drought suppressed all gas fluxes and vegetation activity (Fig. 7). The drought resistance indices, calculated for May–September, indicated average reductions of 15 %, 37 %, and 43 % in daily GPP, ET, and ER, respectively, relative to 2017 (Fig. 9). Although the resistance-based GPP decline in 2018 (Rt = 0.85) appears to contradict the negligible difference between cumulative May–September values in 2017 (−1258 g C m⁻² per period) and 2018 (−1215 g C m⁻² per period), this discrepancy results from the difference in temporal scale. Seasonal sums smooth out short-term variability: in 2018, higher GPP during the early growing season partly compensated for the pronounced reductions in mid- to late summer (Fig. 3), resulting in comparable seasonal totals. In contrast, the resistance index, derived from daily values, captures these episodic declines more accurately, reflecting the stronger suppression of photosynthesis during the drought period.

The reduction of GPP in summer 2018 is in line with observations from multiple sites across Europe (Fu et al., 2020; Lindroth et al., 2020) and can be attributed to stomatal regulation under the lack of soil water availability. Indeed, we estimated a 61 % decline in daily G_c over May–September 2018. After the initial increase in May-June, it continuously declined through July and remained low until the end of the growing season (Figs. 4d, 7). This explains the suppressed canopy photosynthetic capacity in August 2018, especially when compared to the peaking values in the reference year (Fig. 4b).

The G_c suppression was likely driven by high atmospheric demand (i.e. increasing VPD) rather than soil moisture depletion, as indicated by uniformly low values of G_c,ref across all SSR classes in the active season of 2018. High VPD can override soil moisture gradients, forcing sustained stomatal downregulation regardless of soil moisture variation (Novick et al., 2016). The reduction of stomatal conductance to prevent water loss has been previously documented across multiple plant species and forest types (Farquhar and Sharkey, 1982; Fu et al., 2020; Lindroth et al., 2020; van der Molen et al., 2011; Novick et al., 2016; Reichstein et al., 2002), reflecting a conservative hydraulic strategy of our study site.

Since both GPP and ET are mechanistically connected with stomatal regulation, the low ET resistance is not surprising; however, its sensitivity was much higher than that of GPP (37 % decline compared to only −15 %). Total ET includes both transpiration (T) and evaporation, though the latter is limited under drought. However, even when filtering only for the periods with maximum T contribution, the drop in daily values (−29 %) still exceeded that of GPP (Fig. 9). This additionally contradicts the findings of Lindroth et al. (2020), where the majority of sites demonstrated an increase in ET in the drought year. Boese et al. (2019) found that the sites with high seasonal dryness variability experienced a lower ET decrease rate during the progressive drought due to plant adaptations such as deeper root systems to access the water. However, at our site, high soil moisture variability is skewed towards flooding rather than drying, which is consistent with the sharp drop in ET over the course of the July progressive drought (Figs. 2, 3).

The greater ET sensitivity compared to GPP resulted in enhanced daily EWUE (+38 %, Fig. 9) over May–September. An elevated EWUE during drought has been previously reported, for example, for a boreal aspen stand in Canada (Krishnan et al., 2006) and a mixed deciduous forest in Switzerland (Wolf et al., 2013). However, responses appear to be species- and site-dependent; for example, no change in EWUE was observed in a Finnish forest under low rainfall conditions (Ge et al., 2014), while a decline in EWUE was reported for a pine forest in Finland under severe drought stress (Gao et al., 2017).

Although daily EWUE declined during progressive drought in July, VPD-normalised EWUE remained elevated throughout the drought, indicating that high atmospheric demand combined with soil moisture limitation drove the observed water use efficiency dynamics. Similar increases in VPD-normalised EWUE under moderate drought have been reported across forest and grassland ecosystems (Beer et al., 2009).

Interestingly, ER resistance was nearly twice as low as that of GPP (Fig. 9), pointing to a strong drought impact on soil microbial and root respiration. This reduction in ER effectively lengthened the period of ecosystem net carbon uptake and, combined with the favourable early-season conditions in May, resulted in the highest annual net carbon uptake of the three study years.

4.5 Drought recovery and carry-over effects

Drought can influence ecosystem functioning well beyond the event itself, with lagged effects persisting for years after water stress has ended (Kannenberg et al., 2020). In forests, such “drought legacy effects” are common and often span three to four years (Anderegg et al., 2015). They could be caused by the carbon depletion due to reduced uptake during the drought (Bréda et al., 2006; McDowell et al., 2008), the cost of repairing hydraulic damage (Anderegg et al., 2015; Kannenberg et al., 2019), changes in the nutrient cycle (Houle et al., 2016; Schlesinger et al., 2016), or from shifts in carbon allocation towards root development or canopy restoration (Arain et al., 2022; Doughty et al., 2014; Hikino et al., 2022), all of which can constrain subsequent tree growth and ecosystem functioning.

Although our study period was too short to assess long-term drought legacies, it allowed us to evaluate recovery and possible carry-over effects in the year following the 2018 drought. In 2019, GPP, ER, and ET reached their highest values of the three study years, both annually and over the active season (Table 1). Recovery and resilience indices indicated full recovery of all fluxes, with daily active-season ET and T even exceeding that of the reference year (Fig. 9).

While soil moisture was declining over the recovery year summer, the evenly distributed precipitation kept the favourable soil saturation rates over all months, contributing to the ecosystem recovery (Fig. 2). Soil water depletion in riparian systems can vary considerably depending on groundwater connectivity, precipitation patterns, and vegetation water use (Capon et al., 2013). In systems with strong hydrological connectivity to groundwater, depletion may be minor; however, under drought conditions or in systems with limited lateral or vertical recharge, significant drawdown can occur (Rohde et al., 2021; Rood et al., 2008). At our site, in the absence of runoff or drainage measurements, we cannot fully quantify the water balance, and our interpretation of groundwater connectivity remains speculative.

The nearly 30 % increase in annual ER in 2019 relative to the pre-drought year, and 55 % relative to 2018, likely reflected a combination of higher spring temperatures in the recovery year and a pronounced June–July peak in ER_ref, which cannot be explained by temperature alone (Fig. 4a). This interpretation is supported by the absence of differences in active-season daily ER between 2017 and 2019 when modelled using each other's temperature response parameters (Fig. 8d–e). Similarly, the resilience index of ER and ER_ref indicated a full rebound (Fig. 9). The suppressed ER observed during the 2018 drought may have led to the temporary accumulation of labile carbon within the ecosystem, providing a larger pool of carbon that could be available for subsequent physiological activity or decomposition. An increase in ecosystem respiration in the year following the drought was also observed in a beech forest in Denmark (Pilegaard and Ibrom, 2020), attributed to accumulated soil organic matter following suppressed heterotrophic respiration during the drought and autumn months, and for a ponderosa pine forest in the USA (Thomas et al., 2009), where the effect was linked to enhanced litter decomposition.

Likewise, the increase in GPP, though less pronounced than in ER, was caused by a combination of vegetation activity and meteorological conditions. In 2019, elevated temperatures and higher radiation compared to 2017 advanced the onset of the growing season. Enhanced GPP_sat suggests that optimal temperature and VPD, rather than light, were the main drivers at that stage. Later in the season, GPP_sat was lower than in 2017, yet cross-year modelled values converged, particularly in August (Fig. 8a), pointing to lower radiation as the primary carbon uptake constraint. As with ER, both GPP_sat and GPP demonstrated full recovery (Fig. 9).

While transpiration and photosynthesis recovered, stomatal conductance did not fully return to the pre-drought level (Rs = 0.88), indicating that subtle physiological constraints persisted despite overall functional recovery, potentially limiting tolerance to future droughts. However, as these indices are based on a single pre-drought reference year, interannual variability in meteorological conditions may bias interpretation. Favourable conditions in 2019, including evenly distributed precipitation and a warm growing season, likely facilitated the rapid recovery, consistent with observations across diverse ecosystems (Schwalm et al., 2017).

In contrast, strong legacy effects on the carbon cycle have been observed following the 2018 drought in other European forests. For example, in a mixed deciduous forest in central Germany, NEP declined by 150 g C m⁻² yr⁻¹ in 2019, with reductions in both GPP (−281 g C m⁻² yr⁻¹) and ER (−132 g C m⁻² yr⁻¹) compared to the previous year (Pohl et al., 2023). European beech forests have exhibited particularly high sensitivity to drought, with observed tree mortality linked to hydraulic failure (Schuldt et al., 2020). More broadly, drought-induced tree mortality can have long-lasting consequences, with post-drought effects often persisting for months or years (Brodribb et al., 2020; Schwalm et al., 2017). Drought-related growth decline and canopy dieback have also been documented in various riparian trees (Kibler et al., 2021; Schnabel et al., 2022; Singer et al., 2013; Stella et al., 2013; Valor et al., 2020). Our site provided no visual or numerical evidence of increased tree mortality in the year following the drought. However, given that drought-induced mortality can manifest with a delay, it remains possible that long-term effects could emerge beyond the period of our study. Future monitoring would be critical to assessing whether the observed recovery is sustained or whether cumulative drought stress could compromise forest resilience over time.

5 Conclusions

The mature riparian grey alder forest under study remained a strong and consistent net carbon sink over three years with contrasting soil moisture conditions. While GPP was comparable to that of similar ecosystems, ER was generally lower, likely due to dense, poorly aerated soils and periodic flooding.

Moderate soil water saturation ratio (40 %–50 %) enhanced ecosystem fluxes, with flux rates generally persisting even at higher saturation levels. In contrast, the 2018 progressive drought mildly reduced GPP and, to a much greater extent, ER, while also suppressing ET. High EWUE and reduced G_c during the drought indicate stomatal regulation that minimised water loss while maintaining carbon uptake. The co-occurrence of elevated temperatures (driving high VPD) and prolonged dry spells (causing progressive drought) in 2018 was particularly detrimental, with the canopy conductance suppression primarily driven by elevated VPD, while soil moisture variation played a larger role in the other two years.

In the year following the drought, the forest exhibited an overall recovery, supported by high, but not extreme, temperatures and evenly distributed precipitation. The intermediate cumulative NEE was a result of elevated ER in spring to early summer, likely due to decomposition of residual organic matter from the preceding year. GPP also increased, although to a lesser extent, with canopy conductance remaining partially suppressed, suggesting a potential vulnerability to multi-year drought events.

Overall, our results demonstrate that this mature alder forest maintained both productivity and resilience following the 2018 drought, with only partial suppression of canopy conductance, highlighting the stability of this ecosystem under variable hydroclimatic conditions.

Appendix A: Energy balance closure

To assess the performance and consistency of turbulent energy flux measurements, we evaluated the energy balance closure (EBC) on a daily timescale for June-August of each study year (2017–2019). The turbulent fluxes were defined as the sum of latent (LE) and sensible (H) heat fluxes measured using the eddy covariance system. Since the components of available energy, net radiation (Rn) and ground heat flux (G), were not measured at our site, we used the following approach.

Rn was approximated using measured incoming shortwave radiation (Rg) and daily albedo values:

$$\mathrm{Rn}=(\mathrm{1}-\mathit{\alpha })\times \mathrm{Rg},$$

where α is surface albedo and Rg is daily incoming shortwave radiation in MJ m⁻² d⁻¹.

Rg was measured at the site, and daily albedo values were derived from MODIS (MCD43A3 v061), (Schaaf and Wang, 2021) using Google Earth Engine. Extracted albedo values were averaged black-sky and white-sky shortwave albedo components to approximate actual albedo under mixed sky conditions. The data were quality-controlled using MODIS-provided QA flags and seasonally averaged, resulting in mean albedo values of 0.161 ± 0.009, 0.154 ± 0.007 and 0.151 ± 0.007 in 2017, 2018 and 2019, respectively.

G was estimated as 5 % of Rn, following common practice for forest soils. We acknowledge that this approximation may slightly overestimate G under certain soil and moisture conditions but provides a reasonable estimate for comparative purposes. Only daytime (Rg > 15 W m⁻²) half-hourly records were included in the daily energy sums to ensure that energy components reflected active turbulent exchange.

While this approach is limited by the absence of direct measurements of Rn and G, it still provides a consistent method for comparing EBC across years. Since longwave radiation was not accounted for, the available energy may be slightly overestimated. Therefore, the true energy balance closure at the site is expected to be slightly higher than reported. Despite these limitations, and although the EBC values are on the lower end, they fall within the expected range for forested ecosystems. Moreover, the similar closure observed in 2017 and 2018 supports the interpretation that reduced evapotranspiration (ET) in the drought year (2018) was not driven by lower EBC but likely reflects actual physiological or environmental responses.

Figure A1 Energy balance closure in June-August of the three study years. Available energy is represented by net radiation (Rn) minus soil heat flux (G), turbulent energy is a sum of sensible (H) and latent heat (LE) fluxes, all aggregated over daytime periods (Rg > 15 W m⁻²). The energy balance closure is expressed as the slope of the least squares regression, shown in red colour in each panel. The dashed lines are 1:1.

Appendix B: Sensitivity of canopy conductance to VPD in different soil saturation classes

Table B1 Parameters of Oren et al. (1999) canopy conductance sensitivity model.

^* Soil saturation ratio (SSR) classes, with a low slope corresponding to the low quality of the fit.

Figure B1 Oren fits with variable m to different soil saturation ratio (SSR) classes.

Figure B2 Oren fits with a fixed $m/G_{\mathrm{c},\mathrm{ref}}=\mathrm{0.6}$ to different soil saturation ratio (SSR) classes.

Appendix C: Site soil properties

Figure C1 Soil water content at the study site (blue) and SMEAR Estonia station (red).

Table C1 Soil physical-chemical properties of the study site. Mean and standard error values (in parentheses). DM – dry matter, SOM – soil organic matter, TC – total carbon, TN – total nitrogen.

Appendix D: Estimation of NEE, GPP and ER from carbon budget data

For comparison with eddy-covariance estimates, net ecosystem exchange (NEE) from the biometric, inventory-based carbon budget was obtained as the negative of net ecosystem production (NEP):

$$\mathrm{NEE}=-\mathrm{NEP}$$

Gross primary production (GPP) and total ecosystem respiration (ER) were derived from annual net primary production (NPP) and heterotrophic respiration (R_h) estimates using a fixed autotrophic respiration fraction. Autotrophic respiration (R_a) was assumed to represent 57 % of GPP, following the synthesis of boreal and temperate forest carbon budgets by Lindroth et al. (2020). Under this assumption, NPP can be expressed as:

$$\mathrm{NPP}=\mathrm{GPP}-R_{\mathrm{a}}=\left(\mathrm{1}-a\right)\times \mathrm{GPP}$$

where a is the fractional contribution of autotrophic respiration to GPP (here a=0.57). GPP was therefore calculated as:

$$\mathrm{GPP}={\displaystyle \frac{\mathrm{NPP}}{\mathrm{1}-a}}$$

Autotrophic respiration was then obtained as:

$$R_{\mathrm{a}}=a\times \mathrm{GPP}$$

Total ecosystem respiration was computed as:

$$\mathrm{ER}=R_{\mathrm{h}}+R_{\mathrm{a}}$$

NEP from the carbon budget was then used as a consistency check:

$$\mathrm{NEP}=\mathrm{GPP}-\mathrm{ER}$$

Appendix E: Contribution of N₂O and CH₄ to the CO₂-based global warming potential of the mature alder forest

Table E1 The annual balance of the alder forest in 2018 and 2019 expressed in CO₂ eq. (GWP100 = 1 for CO₂, GWP100 = 27.9 for CH₄ and GWP100 = 273 for N₂O; IPCC, 2021).