Interannual Variability of Water and Heat Fluxes in a Woodland Savanna (Cerrado) in Southeastern Brazil: Effects of Severe Drought and Soil Moisture

Abstract

The Brazilian Cerrado biome is known for its high biodiversity, and the role of groundwater recharge and climate regulation. Anthropogenic influence has harmed the biome, emphasizing the need for science to understand its response to climate and reconcile economic exploration with preservation. Our work aimed to evaluate the seasonal and interannual variability of the surface energy balance in a woodland savanna (Cerrado) ecosystem in southeastern Brazil over a period of 19 years, from 2001 to 2019. Using field micrometeorological measurements, we examined the variation in soil moisture and studied its impact on the temporal pattern of energy fluxes to distinguish the effects during rainy years compared to a severe drought spell. The soil moisture measures used two independent instruments, cosmic ray neutron sensor CRNS, and FDR at different depths. The measures were taken at the Pé de Gigante (PEG) site, in a region of well-defined seasonality with the dry season in winter and a hot/humid season in summer. We gap-filled the energy flux measurements with a calibrated biophysical model (SiB2). The long-term averages for air temperature and precipitation were 22.5 °C and 1309 mm/year, respectively. The net radiation (Rn) was 142 W/m², the evapotranspiration (ET) and sensible heat flux (H) were 3.4 mm/d and 52 W/m², respectively. Soil moisture was marked by a pronounced negative anomaly in the 2014 year, which caused an increase in the Bowen ratio and a decrease in Evaporative fraction, that lasted until the following year 2015 during the dry season, despite the severe meteorological drought of 2013/2014 already ending, which was corroborated by the two independent measurements. The results showed the remarkable influence of precipitation and soil moisture on the interannual variability of the energy balance in this Cerrado ecosystem, aiding in understanding how it responds to strong climate disturbances.

1. Introduction

The Brazilian Cerrado is one of the largest biomes in South America, covering 1.5 million km² or 24% of the Brazillian land, composed of grasses and savanna forests and playing a fundamental hydrological role by containing headwaters of large rivers that recharge large deep aquifers, with a great spatial climatic variation with average temperatures ranging from 20 °C to 26 °C and precipitation from 800 to 1800 mm/year [1]. The vegetation of the Cerrado is characterized by strong seasonality, with a humid and hot summer and a dry and mild winter, which markedly controls evapotranspiration and primary productivity [2,3]. This seasonal control is generally noted in the surface energy balance pattern, especially in the partitioning of turbulent atmospheric fluxes of sensible heat H and latent heat LE, which is very different from the Amazonian forests [4]. The partitioning of available energy, for example, expressed in the Bowen ratio (β = H/LE), shows high values in the dry season and small values in the rainy season, respectively [4], associated with the high soil water deficit during the dry season [5]. The most common vegetation physiognomies in the Cerrado biome are the campo Cerrado and Cerrado sensu stricto, structurally composed of herbaceous and arboreal strata that coexist simultaneously [6]. During the dry season, many Cerrado trees show adaptation mechanisms in their senescence process to try to use moisture from deeper soil levels [7], while grasses undergo dormancy processes. Fire events in the Cerrado, which are part of its adaptive mechanisms in its evolutionary history, also influence the energy balance and soil water status temporarily [5,8].

Furthermore, on a regional scale and especially at the border of the Amazon and Cerrado biomes, deforested areas replaced by cultivated pastures have possibly promoted the reduction of evapotranspiration (ET), which has possibly influenced the reduction of rainfall at the beginning of the rainy season in recent decades [9]. In general, in terms of natural precipitation variability, Brazil has been subject to large regional-scale meteorological droughts in the last two decades, such as in the Pantanal biome in 2020 [10], in the Cerrado biome of the central-west region between 2016 and 2017 [11] and particularly in the southeastern region from 2013 to 2014 [12] with a rainfall reduction of 44%, considered an exceptional index that disrupted the socio-economic condition of 80 million people due to water consumption restrictions [13].

Soil moisture is a condition of the Cerrado state that is inherently linked to vegetation seasonality, energy partitioning at the surface and possibly interannual climate variability, whose variation depends on the precipitation, evapotranspiration and water storage in the soil [3,4,5,14]. The measurement of soil moisture can use different methods, such as direct current resistivity [15], electromagnetic induction [16] and time-domain reflectometry (FDR) [17], which, however, sample a small volume of soil at a point scale (ANNISS, 2021). The cosmic ray neutron soil moisture probe (CRNS) measures neutrons that have been thermalized by hydrogen atoms present in the soil [18], potentially within a footprint of 130–240 m radius around the probe and extending to depths from 0.15 m up to 0.83 m [19]. Joint measurements of soil moisture and energy fluxes in the Cerrado are scarce in the literature, and knowledge of the interactions between soil moisture and vegetation remains an open question, with several implications: How sensitive are the tree and grass layers to the seasonal water deficit? How does soil water deficit control the vegetation response during season transitions (from wet to dry and from dry to wet)? What is the impact of major meteorological droughts on soil moisture and energy partitioning?

Our work aimed to evaluate the seasonal and interannual variability of the surface energy fluxes in an area with woodland savanna ecosystem (Cerrado) in southeastern Brazil, using micrometeorological measurements, in a range during 19 years from 2001 to 2019, and discusses the soil moisture variability in distinguishing the effects in rainy years compared to severe drought years.

2. Materials and Methods

2.1. Study Area and Micrometeorological Measurements

Our study area is located in the city of Santa Rita do Passa Quatro, São Paulo state, Brazil, in the Pé de Gigante Forest Reserve (PEG), with an area of 1060 hectares, covered with woodland savanna vegetation (Figure 1, with dense Cerrado predominantly of Cerradão and Cerrado sensu stricto physiognomies [20,21] and other small areas of campo and deciduous forest, with a surrounding eucalyptus plantation to the west and north, and crops to the east and south.

The site is located ≃ 300 km from the coast and 700 m above sea level, with climatological normals (interval 1961–1990) at the Santa Rita do Passa Quatro meteorological station (21°43′09″ S; 47°28′22″ W, 715 m) reporting Cwa climate (Koeppen classification), average temperature ranging from 17.6 °C (July) to 23.5 °C (February), a dry season centered in winter from May to August and a rainy season centered in summer, from September to April with average precipitation of 1478 mm/year. The canopy height of vegetation in PEG varies from approximately 10 to 12 m. The soil in the region is predominantly Neossolo Quartzarênico (94%), sandy with good rooting capacity and low natural fertility, with reduced levels of nutrients. Its superficial layer horizon A is characterized by the decomposing organic layer, and horizon B presents a higher concentration of clay [23].

Micrometeorological measurements (

Table 1. Field measurements at PEG site.

Variable (Unit)	Instrumentation
Air temperature and humidity (°C,%)	CSI HMP45C
Precipitation (mm)	Hydrological Services TB4
Horizontal wind speed (m/s)	RM Young
Net radiation (W/m2)	REBS Net-Lite
Latent and sensible heat flux (W/m2)	Eddy correlation using sonic anemometer CSAT3 CSI and Open path gas analyzer LI-7500
Soil moisture (m3/m3)	Reflectometer FDR CS615 (2001–2009) and CS616 (2010–2018); CRNS (Hydroinnova LLC)

) were performed automatically between 2001 and 2019 on two micrometeorological towers, the first after 2001 (21.619238° S, 47.632327° W) at 21 m height, and the second (21.618609° S and 47.633810° W) from 2010 on and about 100 m southeast of the first, at 13 m height. The positioning was designed so that the tower fetch took advantage of the predominant wind direction, usually from the east and northeast during the day and from the southeast at night, and thus representative of the Cerrado vegetation area.

With the installation of a 2.5 m horizontal mast at 18 m high, we measured net radiation (Rn, LITE, Kipp and Zonen, Delft, The Netherlands), global solar incoming radiation (Rg, CM3, Kipp and Zonen) and the incident and reflected fluxes of photosynthetically active radiation (PAR, LITE, Kipp and Zonen). At a height of 21 m, we measured precipitation (TB4, Hydrological Services Pty. Ltd., Sydney, NSW, Australia), wind speed (RM Young) and air temperature and humidity (HMP45 Vaisala, Helsinki, Finland). The soil heat flux (G) was measured using four heat plates (REBS, Seattle, WA, USA). The turbulent heat fluxes (LE) and (H) were estimated using the Eddy Covariance (EC) method, using a sonic anemometer (CSAT3 Campbell Sci) and an open-path gas analyzer (LI-7500 Licor Inc., Lincoln, NE, USA), with acquisition at 10 Hz with a CR1000 datalogger. The data processing was performed using EddyPro 7 software (Li-Cor, Lincoln, NE, USA). The data were converted from binary format to ASCII using Card-Convert software (Campbell Inc., Logan, UT, USA) and processed in 30 min blocks. To calculate the fluctuations, block averaging, linear trend elimination and moving average were applied. Time lag compensation was performed by maximizing the covariance, axis rotation for the correction of inclination through double rotation calculation, fluctuation and density compensation [24], sensor heating correction and high-pass and low-pass filters https://www.licor.com/env/support/EddyPro (accessed on 10 January 2001).

To address uncertainties in the estimation of energy fluxes due to energy balance closure, we used the correction of [25], assuming that the available energy measurements represented the coverage area of the EC and using criteria based on the observed Bowen ratio:

β = Ho/LEo

(1)

where Ho = observed sensible heat flux; LEo = observed latent heat flux; and the corrected fluxes Hc and LEc were estimated as:

LEc = (Rn − G − Sa − Sb)/(1 + β)

(2)

Hc = Rn − G − LEc − Sa

(3)

where Rn = net radiation; G = soil heat flux; Sa = sensible and latent heat storage change in the air column.

2.2. Soil Moisture

The soil moisture was measured by two instruments, one using vertical profile scanning with frequency domain reflectometry (FDR) (CS615 Campbell Sci) at depths of 10, 20, 50, 80, 100, 150, 200 and 250 cm, from October/2000 to 2009 and CS616 (Campbell Sci) from 2010 to 2019, and the second using horizontal/vertical scanning with a cosmic ray neutron sensor (CRNS) (Hydroinnova LLC, Albuquerque, NE, USA), installed in July/2011 and available until August/2018. The sensors were calibrated with local soil gravimetric samples (BRUNO 2004). The CRNS sensor represents a technological innovation that reaches a much wider horizontal extent of surface soil moisture and prioritizes estimates at shallower soil depths. We are not aware of previous reports of publications using CRNS measurements in savanna areas in Brazil. The estimation of soil moisture by CRNS went through calibration steps, due to interactions with water molecules present in the atmosphere, biomass and minerals. The raw data from the CRNS sensor were obtained from the website http://cosmos.hwr.arizona.edu/Probes/StationDat/044/index.php (accessed on 2 February 2012), with the following variables used: neutron count (MOD), correction factor that considers the size of the probe (PROBE), correction factor due to atmospheric pressure (PRESS) and correction factor that considers changes in the intensity of cosmic rays as a function of time (INTEN). Air temperature, relative humidity and precipitation data from the EC tower located at PEG, and the calculation of absolute air humidity (ρ_v) were also used in the correction process. In the first correction step, based on the neutron count according to vapor pressure (fh), it was assumed that:

Fh = 1 + 0.0054 × ρ_v;

(4)

The second step, biomass, was based on the process described by [26]. Literature values for the Brazilian Cerrado region [27] were also used to establish the average soil organic matter content. As a result of the processes described above, we obtained the corrected neutron count (Npihv). Thus, soil moisture was obtained by the equation:

Θ = (α0/(Npihv/N_0 − α1) − α2 − lw − SOC) × BD;

(5)

α0 = 0.0808 (cm³/g), α1 = 0.372 and α2 = 0.115 (cm³/g) are fixed, with lw = 1.7% being the water content in soil minerals, SOC is the water content in soil organic matter, BD is the soil bulk density and N_0 is obtained from the available calibration data at (http://cosmos.hwr.arizona.edu/Probes/StationDat/044/calib.php, accessed on 25 May 2024).

After correction, the calibration process was based on soil moisture data obtained by the gravimetric method.

The daily averaged volumetric water content data (θ, in m³/m³) were used to evaluate the Soil Wetness Index (SWI) and soil wetness (W), respectively, given by:

SWI = ((θi − θmin))/((θmax − θmin))

(6)

W = θi/θs

(7)

θi = soil moisture measured on day i; θmin = absolute minimum moisture content of the daily average series; θmax = absolute maximum moisture content of the daily average series; and θs = soil porosity, described here as θmax.

We noticed that the comparison of the CRNS sensor showed a strong correlation with the FDR sensor measurements at the shallower soil depths of 10 cm and 20 cm, respectively, as shown in Figure 2, while at deeper FDR depths (50 to 250 cm), the correlation was weak (not shown). Knowing in advance that the vertical sampling of the CRNS is relatively shallow, it was useful to note that the FDR sensors, which have a weakness in horizontal sampling, compared well with the CRNS at the shallower depths. This suggests that there was adequate representation of surface moisture in the study area with both measurement approaches, at least in the shallow soil.

2.3. Gap-Filling the Observational Series

To fill in the gaps in the historical series of energy fluxes at hourly resolution, we used the SiB2 model [28] as a predictor of net radiation and atmospheric fluxes H and LE. The SiB2 model was run in offline mode, forced with field meteorological data (incoming solar radiation, air temperature and humidity, horizontal wind speed, precipitation and CO₂ concentration) at a time step of 1 h, and it calculated the fluxes of net radiation, sensible heat flux, latent heat flux and its components (soil evaporation, transpiration, and rainfall interception evaporation), soil heat flux, soil temperature and soil moisture and gross and net CO₂ assimilation. The meteorological data used to force the model were checked for consistency to eliminate spurious values due to sensor malfunctioning. The inconsistencies and gaps were filled using data from the ERA5-Land reanalysis [29], downscaled to the micrometeorological tower scale, with a linear regression method of hourly percentiles of measurements with the reanalysis. From 2001 to 2019, the filling totaled 5.1% for incoming solar radiation, 4.4% for air temperature, 4.6% for water vapor pressure and 7.6% for horizontal wind speed.

2.4. Estimating Model Parameters with Remote Sensing

The approach involved prescribing some model parameters from remote sensing of leaf area index (LAI) and the fraction of absorbed photosynthetically active radiation (FPAR) [29]. As the first step in parameter estimation during calibration, the time-variant parameters (monthly based) of LAI, FPAR, canopy greenness fraction (N) and land cover fraction (V) were estimated using the adapted MAPPER algorithm [30], as:

LAI = LAImax * NDVI

(8)

FPAR = Vc N [1 − exp(−kLAI/V)];

(9)

N = LAI/(LAI + Lst)

(10)

V = LAI/LAI_max

(11)

where k = mean extinction coefficient of photosynthetically absorbed radiation in the vegetation canopy; Lst = stem and trunk area index; LAI_max = maximum leaf area index; NDVI = normalized difference vegetation index.

The NDVI data were obtained of Landsat-8 time series from 2013 to 2020, with a 30 m spatial resolution. Due to the relevance that vegetation cover has on the energy partitioning of ecosystems [31], we used maps with detailed vegetation type [32].

The parameters described in

Table 2. Initial and optimized parameters for optimization of net radiation and for the LE flux (physiology).

Parameter—Radiation	Default		Optimized
Leaf transmittance (PAR green leaf)	0.050		0.0100
Leaf transmittance (NIR green leaf)	0.250		0.3050
Leaf transmittance (PAR dry leaf)	0.001		0.0050
Leaf transmittance (NIR dry leaf)	0.001		0.0050
Leaf reflectance (PAR green leaf)	0.100		0.0100
Leaf reflectance (NIR green leaf)	0.450		0.1000
Leaf reflectance (PAR dry leaf)	0.130		0.2050
Leaf reflectance (NIR dry leaf)	0.390		0.3500
Soil reflectance (PAR)	0.080		0.1550
Soil reflectance (NIR)	0.200		0.2050
Leaf angle distribution factor	0.250		0.0500
Parameter—physiology	Default JFM/JJA	Optimized JFM	Optimized JJA
Photosynthetic-conductance slope factor (m)	8.00	8.29	12.95
Average leaf turnover time (gmudmu)	0.90	0.51	0.30
Canopy greenness fraction (N)	0.90/0.60	0.99	0.43
Maximum Rubisco assimilation (Vmax, μmolm−2s−1)	100.0	120.0	119.9

influence the time-mean PAR extinction coefficient $\overline{k}$ (12), and consequently, the fraction of canopy PAR absorbed (FPAR) (13), the CO₂ assimilation (14–16) and the water vapor canopy conductance (17), as follows:

$$\overline{k}=\overline{\left[\frac{G\left(\mu \right)}{\mu }\right]}{\left[1-w_{\pi }\right]}^{1/2}$$

(12)

where $\overline{\left[\frac{G\left(\mu \right)}{\mu }\right]}$ = time-mean leaf projection; $w_{\pi }$ = leaf dispersion coefficient of photosynthetically active radiation absorbed by the canopy.

$$\Pi \approx \frac{FPAR}{\overline{k}}$$

(13)

where $FPAR$= fraction of photosynthetically active radiation absorbed by the canopy; k = term belonging to the extinction coefficient (($\overline{k}$)).

$$A=A_{n0}{\int }_0^{LAI/V}VN{e}^{\overline{-kLAI}}dL=A_{n0} \Pi$$

(14)

A_n0 = f(Wc,We,Ws) − Rp

(15)

where $A_n$ = net leaf CO₂ assimilation; Wc = photosynthetic enzyme RuBisCO efficiency; We = amount of PAR captured by chlorophyll; Ws = leaf capacity to export products of photosynthesis; and Rp = plant respiration that depends on the Rubisco-limited rate of assimilation $v_m$ (mol/m²s),

$$v_m=v_{max}f\left(Tc\right)f\left(W\right)$$

(16)

where $v_{max}$ = maximum catalytic capacity of Rubisco (mol/m²s); $f\left(Tc\right)$ = canopy temperature stress factor; $f\left(W\right)$ = soil moisture stress factor, and

$$gc=m \frac{A_n}{Cs}{h}_{s}p+b$$

(17)

where $gc$ = leaf stomatal conductance (m/s); $m$ = slope coefficient proportional to leaf water use efficiency; $Cs$ = partial CO₂ pressure on the leaf surface (Pa); $p$ = atmospheric pressure (Pa); and b = intercept coefficient = 0.01 for C3 plants and 0.04 for C4 plants.

The temporal and spatial average NDVI field in summer and winter (Figure 3) showed a large deviation between them within the polygon boundaries of PEG site, predominantly with high values in January around 0.85 and lower values in August around 0.70, as expected due to the strong seasonal effect on Cerrado vegetation. In both months, a small area with a field physiognomy, close to the center of the polygon, showed a much-reduced NDVI compared to the total area, with the absence of trees. Except for this localized minimum, there was a certain spatial heterogeneity of NDVI, for both months, but that did not correspond to a very clear association with the patterns of other dominant physiognomies in the area (Figure 1D). The establishment of correlation between the Cerrado physiognomies proves to be complex, not only in relation to NDVI as previously mentioned, but also floristically and by the analyzed soil characteristics [33]. Reference [23] also reinforced that for PEG, the detailed taxonomic classes of soils were also insufficient for a complete understanding of the spatial variation of vegetation.

Similarly, the averaged LAI (Figure 4) and FPAR (Figure 5) in summer and winter also showed modest but significant seasonal variation. In the total area covered with Cerrado, the LAI varied from approximately 4.60 m²/m² in January to 4.25 m²/m² in winter, as shown in Figure 4b. The FPAR parameter also varied in the same direction, from approximately 0.975 in January to 0.955 in August. The seasonality of vegetation in the site was already noticed in situ, with increasing litterfall during the dry season [34]. The leaf fraction, the most important component of litter, is extremely sensitive to meteorological variables such as precipitation, humidity and air temperature, that are possibly associated to structural changes of the canopy [20].

2.5. Model Calibration

The model calibration proceeded with the prescription of time-variation parameters (Section 2.4) and soil physical properties (Table 2). The next steps were the optimization of parameters with the calculation of the net radiation (Rn) and latent heat flux (LE). The optimization was based on the method of Non-Linear Least-Square Minimization and Curve-Fitting for Python (LMFIT) [35], which minimizes the squared deviation of the simulated variable from the measurement. For a more accurate and computationally efficient search for the minimum error space across the range of each parameter to be optimized, the fit was performed in two steps. In the first step, we calculated the deviations across a thresholded range of the parameter (min and max) for discrete values of the parameter equally spaced in the range, called the brute computational estimate. In the second step, from the results of the first step, a sub-range was selected for refined search of the minimum deviation, using the Levenberg–Marquardt method (leastqr routine), under a nonlinear procedure adjustment. In summary, the calibration followed serial steps: (i) estimation of time-variant parameters LAI, N and FPAR with remote sensing; (ii) prescription of time-invariant parameters; (iii) optimization of optical parameters for net radiation; (iv) optimization of physiological parameters for latent heat flux.

The calibration of optical parameters with optimization of the net radiation (Rn) used data in the time interval from 2010 to 2017 and was performed simultaneously for 11 parameters (Table 2). For the LE flux, the optimization was performed in the summer months from January to March, and in the winter months from June to August, to accommodate the seasonality of the optimized parameters and the soil moisture initialization in the year 2011 with selected days without precipitation. The optimized parameters of slope of photosynthetic conductance, average leaf time projection and green fraction of the canopy were interpolated for the other months of the year. The improvement in performance with calibration was quite adequate for Rn (Figure 6) and LE (Figure 7).

3. Results

3.1. Climatic Characterization

We discussed the climate variability of monthly average field measurements of incoming solar radiation, air temperature and humidity, wind speed and precipitation between 2001 and 2019. In general, all variables showed a well-marked seasonality (Figure 8), with pronounced differences between dry and rainy seasons (

Table 3. Long-term average of meteorological and flux variables for the range 2001–2019.

Season	Air Temperature (°C)	Precipitation (mm/year)	Rn (W/m2)	LE (W/m2)	H (W/m2)	G (W/m2)	ET (mm/d)	Et (mm/d)	Es (mm/d)	Ei (mm/d)
Annual	22.1	1309	127.9	100.6	48.8	−1.8	3.5	2.3	0.7	0.4
Rainy season	23.3	1108	147.3	118.7	51.8	−1.6	4.0	2.6	0.9	0.5
Dry season	20.0	128	88.6	65.1	42.9	−2.0	2.2	1.7	0.3	0.2

), starting with the solar radiation, which ranged from ≃ 150 to 300 (W/m²) from winter to summer, respectively, representing a two-fold factor. The seasonality of rainfall, humidity and wind was also quite pronounced, with vapor pressure maxima/minima occurring in January/August, ranging from ≃10 to 21 hPa (Figure 8c), and wind speed maxima/minima occurring in August/February, ranging from ≃ 2 to 3.5 m/s on average (Figure 8d).

The average rainfall was 1309 mm/year (Figure 8e, Table 3), slightly below the climatological normal (1478 mm/year) (Figure 9), that was most recurrently caused from the year 2010 on, when the rainfall was above normal only in 2016. The year 2014 was notable for the meteorological drought and the largest negative deviation of 488 mm/year, followed by several years of below-average rainfall that persisted until 2019.

Throughout the period, the average temperature was 22.1 °C, with extremes ranging from ≃17 °C to 26 °C on a monthly basis. During 2014, we noted the highest annual average temperature of 23.2 °C (Figure 9).

We also noted a pronounced variability in the mean diurnal cycle and concurrent with seasonality, based on the comparison of the hot/wet months January–February–March (JFM) and the cold/dry months July–August–September (JJA) (Figure 10). In JFM, the incoming solar radiation (Figure 10a) was higher, caused by higher elevation angle compared to JJA and by the greater hourly variability due to the presence of cloudiness.

The hourly median temperature (Figure 10c) varied between 22 °C and 28 °C in JFM, and from 15 °C to 25 °C in JJA, respectively, showing a larger amplitude in winter due to the cooler minima. The water vapor pressure median was generally lower in winter, ranging from about 22 hPa in JFM to 12 hPa in JJA, with little daily variation in both seasons (Figure 10b). Likewise, the minimum daily median relative humidity varied between about 37% in JJM, with many events reaching as low as 20%, and 60% in JJA, when it never fell below 40% (Figure 10e). The average wind speed was 3.5 m/s in JJA, higher than the 2.5 m/s observed in JFM (Figure 10d), reflecting the seasonality of large-scale wind patterns. During summer, wind speeds tended to increase during the daytime due to stronger radiative forcing.

3.2. Soil Moisture

We compared the mean monthly soil moisture measured by CRNS with FDR measurements estimated as an average between depths of 10 and 100 cm (referred to as FDR10–100 cm) and an average between depths of 150 and 200 cm (referred to as FDR150–200cm) in Figure 11. The seasonality of SWI (Soil Water Index) demonstrated a similar variability pattern between the CRNS, FDR10–100 cm and FDR150–200 cm. The minimum Soil Wetness Index (SWI) typically occurred around August, which usually marked the end of the dry season. With the onset of the rainy season after September, soil recharge began, usually reaching its maximum in January, which coincided with the rainiest month. From January onwards, there was a gradual depletion of soil moisture in both the CRNS and FDR above 200 cm, concurrent with a reduction in rainfall.

For soil moisture at greater depths down to 250 cm, the pattern indicated that soil recharge, after the onset of the rainy season, was slightly slower than at shallower depths. This delay is typical in tropical regions with highly permeable soils, as also observed in Amazonia [36].

Given the distinct seasonal patterns above and below the 200 cm depth, we restricted the calculation of the average annual SWI for the aforementioned layers to the months of January–February–March (JFM), July–August–September (JAS) and October–November (ON). These ranges correspond to stages of maximum and minimum seasonal moisture, and soil recharge, respectively. Overall, there was a pronounced decline in soil moisture in 2014 due to the meteorological drought, particularly in the summer (JFM), as estimated by the FDR (Figure 12a—black line) across all soil layers and by the CRNS (Figure 12b—black lines). During the winter (JAS) and seasonal transition (ON), 2014 also showed low moisture levels; however, the depletion of moisture apparently extended into the following years (2016 and 2017), in contrast to the summer when there was some recovery.

3.3. Turbulent Fluxes of Heat and Water Vapor

On average, over the interval 2001–2019, the evapotranspiration (ET) calculated by the SiB2 model was partitioned into the components, ET = Et + Es + Ei, that showed transpiration (Et) with the greatest contribution of 68%, while soil evaporation (Es) and interception losses (Ei) contributed with 20% and 12%, respectively (Figure 13).

In terms of seasonality, the component partition was close to the annual average during the rainy season; however, in the dry season the percentages change a little, with a drop in Es (11%) and Ei (7%), while Et increased (81%) (Figure 13a) as a result from less rainfall and surface soil moisture.

The mean calculated ET was 3.46 mm/d, which is close to another estimate (3.36 mm/d) at the same experimental site, despite of the different temporal interval used (2010–2012) [14].

The variation of monthly net radiation, LE and H, showed a clear seasonal pattern and significant difference in magnitude between LE and H fluxes (Figure 14). The variation in monthly net radiation, LE and H showed a clear seasonal pattern and a significant difference in magnitude between LE and H (Figure 14).

During the hottest and most humid months, net radiation was higher due to increased global solar radiation. Combined with greater soil moisture, this resulted in a higher latent heat flux. A variation in the median LE was observed around 125 W/m², with a wide range of intra-month variance. Additionally, during this period, the H flux tended to be lower throughout the year (Figure 14).

In winter and during the late dry season (August–September), the median H flux varied between approximately 40 and 80 W/m², with intra-monthly values reaching up to 140 W/m². Meanwhile, LE showed minimum values throughout the year, close to 60 W/m², and exhibited low intramonthly amplitude (Figure 14).

A break in the temporal pattern of the LE and H fluxes occurred particularly in 2014, during the drought, which extended into 2015. In the dry seasons of 2014 and 2015, the H term was significantly higher than in other years, while the LE term was lower, as noted in the monthly variation (Figure 15). A common factor in these two years was the marked reduction in soil moisture (Figure 15), indicating the lasting effects of the 2014 drought that also manifested in the following year’s dry season.

This disturbance in turbulent energy flux due to drought was further analyzed using two other indices: the evaporative fraction (EF = LE/Rn) and the Bowen ratio (β = H/LE). Before 2014 and again from 2016 onwards, the EF term typically varied between 0.5 in the dry season and 0.75 in the rainy season, indicating greater efficiency in using available radiative energy for evapotranspiration (ET) during the summer, and the opposite in winter. Similarly, the β term varied from approximately 0.25 in the rainy season, exhibiting quite stable patterns, and increased to values around 1 at the end of the dry season, usually in the form of peaks (Figure 15). During anomalous events associated with meteorological drought, it was noted that the β value exceeded 2 during the dry season of 2014 and went above 4 in 2015, indicating a pronounced inter-annual oscillation. Between 2009 and 2012, no estimates of β above 2 were reported [14]; these estimates only manifested later with the event of a great drought.

4. Discussion

Over the 19-year sampling period, we observed a decrease in annual precipitation in our investigation site since 2010, with more emphasis from 2014 on, when recurrent negative rainfall anomalies were manifested. These were not just local events but were likely serial droughts on a regional scale in southeastern Brazil, that is, of interannual natural variability and manifested as hydrological droughts at larger scales [37]. Our results have shown a significant influence of soil moisture on the response of ET in the Cerrado denso ecosystem. This effect was observed not only on a seasonal scale, as is well documented [2,4,5,14], but also on an inter-annual scale, where clear temporal propagations were evident.

In particular, during the dry seasons of 2014 and 2015, we observed minimum β owen ratio (β) values above 2. It is acknowledged that ET in the Cerrado ecosystem is significantly influenced by the available surface energy [4,14]. However, soil moisture is also closely correlated to the net radiation, resulting in complex and interdependent control mechanisms where variables do not act independently [3]. This understanding is partially attributed to the vegetation’s adaptive characteristics, particularly to water stress. The seasonal phenological changes in the vegetation, including foliar senescence in most tree species and dormancy in the herbaceous layer, typically occur during the dry season. These changes lead to a reduction in the vegetation’s photosynthetic capacity and transpiration rates [2,34].

The temporal variability in green leaf cover composition could account for some of the observed phenomena where heightened soil moisture levels did not lead to proportionate increases in ET. For instance, we documented several isolated rainfall events during the dry season, resulting in significant soil moisture pulses lasting several days, such as those in June 2013, July 2014 and June 2017. Despite these disruptions in soil moisture levels, the declining trend of ET throughout the dry season appeared unaffected, possibly due to a pronounced restriction in transpirational capacity dependent on the state of the vegetation’s foliage.

Various forms of disturbance are widespread throughout the Cerrado ecosystem on a continental scale, including occurrences such as fires followed by subsequent regeneration [38,39,40]. Recurrent fires have the potential to alter the physiognomy of the Cerrado, often leading to the transformation of its landscape into open fields [39], thereby impacting the water balance [40]. At the PEG site, we have no records of burning for approximately 50 years, and thus we precluded this type of influence on our measurements.

5. Conclusions

Our study aimed to offer an original perspective on heat and water fluxes within a dense Cerrado ecosystem. We examined surface energy partitioning and the influence of soil moisture on its temporal variability using a contemporary array of field instruments over an extensive time span. Soil moisture measurements were conducted using two independent systems, with CRNS measurements demonstrating strong agreement with FDR sensors at shallow depths, thereby bolstering their suitability in representing local conditions. Furthermore, CRNS data exhibited temporal synchrony with FDR measurements down to depths of 100 cm, indicating its capacity to effectively capture water conditions across a substantial soil profile of up to 1 m in depth. This enabled us to explore the impact of significant climate events, such as a major meteorological drought.

In general, the estimates showed the long-term average Rn of 142 W/m², ET of 3.5 mm/d and the H flux of 52 W/m². ET was heavily seasonal, influenced by surface energy availability, but soil moisture was also tightly correlated to net radiation, creating complex control mechanisms. The vegetation’s adaptive traits, also in response to water stress, contributed to this understanding.

We emphasize the soil moisture’s sensitivity to significant rainfall anomalies, such as those observed in 2014. Interestingly, the soil moisture anomaly persisted into the subsequent year despite the alleviation of the meteorological drought. This anomalous soil moisture likely contributed to the sustained elevation of the Bowen ratio, which persisted into the following dry season.

In light of the scientific insights gained, we conclude that a thorough assessment of soil moisture and surface energy fluxes can enhance our understanding of natural tropical ecosystem dynamics and their responses to climate disturbances.