Interannual and Intraseasonal Variability of Drought and Heatwaves and Their Effects on Expanding Soybean Production Regions in Brazil

1. Introduction

One of the most significant challenges of the 21st century is ensuring food security for a global population projected to exceed 10 billion people by 2050 [1]. As one of the most important crops worldwide, soybean (Glycine max) plays a critical role in addressing this challenge by serving as a major source of protein and oil [2,3]. However, soybean production is increasingly threatened by climate variability and the rising frequency, intensity, and duration of extreme hydrometeorological and climate events (EHCEs) [4]. Among these, drought and heatwave events (DHE) are significant constraints, particularly for tropical rainfed agricultural systems that are highly sensitive to climatic anomalies [5,6].

This challenge is especially relevant to Brazil, which has emerged as the world’s largest soybean producer, harvesting 169 million tonnes (~55% of global exports) in the 2024/25 season [7,8]. Although cultivated across a wide latitudinal range, from subtropical southern states to equatorial regions, Brazilian soybean yields remain highly vulnerable to climate extremes, particularly droughts [9,10,11]. Such conditions cause substantial interannual yield variability [7]. According to [12], in Southern Brazil, recurrent climatic anomalies have resulted in exceptional losses in grain production. Consecutive cropping seasons with diminished yields, particularly in 2020, 2022, and 2023, led to unharvested grains valued at approximately 8.2 billion USD (about 45 billion BRL).

The intensity and duration of droughts are rising amid increasing climate variability, underscoring the need for refined tools to monitor spatiotemporal patterns [13,14]. Process-based crop models have been widely used to simulate yield responses [15,16,17]. Nevertheless, their reliance on site-specific inputs (soil properties, cultivar traits, and management practices) often introduces calibration biases. More critically, their simplified treatment of extreme events tends to underestimate both the magnitude of climate risks and their spatial heterogeneity [18,19].

Among the available approaches to estimating drought, the Standardized Precipitation Index (SPI) quantifies precipitation over a given time scale as a standardized deviation from the long-term mean, and the Standardized Precipitation Evapotranspiration Index (SPEI) incorporates potential evapotranspiration alongside precipitation to assess available water resources. These indices provide valuable tools for characterizing drought [20,21]. They have been widely applied worldwide to assess drought trends, variability and impacts, including applications across Brazil [22,23,24].

In addition to water deficit indicators, the Warm Spell Duration Index (WSDI) offers a complementary perspective by capturing the frequency and duration of heat events, making it a valuable tool for detecting trends in heat waves and their intensification over time [25]. The WSDI is designed to measure prolonged periods of unusually warm days relative to a baseline climatology. It quantifies the number of days in a year that are part of warm spells lasting at least six consecutive days, with daily maximum temperatures exceeding the 95th percentile of a reference period. Although brief, these events can overlap with droughts and occur during critical stages of plant growth, causing damage to yield and other crop phenotypes. Together, SPI, SPEI, and WSDI provide bias-free insights into the frequency, severity, and spatial extent of drought and heatwave events across multiple timescales.

This study addresses the question: How do droughts and heatwaves affect soybean yields across Brazil’s soybean production regions? The objectives are: a) Integrate temporally variable and spatially distributed climate and soybean yield data at the municipal scale across Brazil’s soybean-producing areas; b) Identify DHE in Brazil’s soybean landscapes using relevant climate indices; c) Quantify the spatiotemporal contributions of DHE to soybean yield losses, analyzing potential drivers and underlying causes in each municipality. The study focuses on Brazil’s key soybean-producing regions, combining high-resolution climate data with SPI, SPEI, and WSDI analyses to characterize regional patterns of drought and heatwaves.

2. Materials and Methods

2.1. Study Area

A wide range of climatic conditions is observed across Brazil’s extensive territory (Figure 1a). It spans from the subtropical southern states within the Pampa and Atlantic biomes, characterized by Cfa and Cfb Köppen-Geiger climate classes, to the low-latitude regions of the Cerrado biome, the country's main agricultural frontier. In these northern regions, soybean “tropicalization” has enabled its expansion into areas with Am and Aw climate classes [26].

In this study, nine soybean-producing municipalities were selected for analyses of drought and heat indices (Figure 1; Table 1). Each municipality represents a homogeneous climatic zone (CZ) as defined by [27] and [28]. This workflow is based on the climate and agricultural representativeness of characteristic soybean production regions at the municipal level across Brazil. Additionally, the workflow and associated analytics can be expanded to include new crops and their corresponding municipalities.

To ensure consistency with the period when the crop is biologically active and water use is agronomically relevant, we used only observations occurring between the local planting start and harvest end dates, as in [29]. The soybean growing season for each municipality was defined using on the official crop calendar published by [30], which provides region specific information on planting, vegetative development, and harvest periods. For each location, the planting window, the onset of vegetative growth, and harvest dates were converted into daily time series, allowing all climate and yield observations to be filtered to the months encompassed by the respective growing season.

2.2. Soybean Yield Dataset

We analyzed annual soybean yield (Yr) data for the selected Brazilian municipalities (Table 1) from the 1996/97 to the 2019/20 crop seasons. For each municipality, the dataset included harvested area (ha) and observed yield (Yr, t ha⁻¹), obtained from the Brazilian Institute of Geography and Statistics [31]. Drought years, used for subsequent analysis of yield reductions, were defined as those with SPEI-12 < –0.5, consistent with previous studies assessing soybean yield losses in Brazil [11,21].

This section focuses on the $Yr$ which captures both interannual variability and long-term trends driven by technological improvements, management practices, and cultivar adoption. To isolate interannual deviations attributable to climate, $Yr$ was detrended using a five-year moving window, which effectively removes short- to medium-term trends while retaining year-to-year variation [11,32]. For periods at the start or end of the series, the window was adjusted to use all available data, with at least one year of observations.

For municipality $i$ and year $t$, the detrended yield ${Ydetrended}_{,i,t}$ was computed as the centered moving average (Equation 1):

$${Ydetrended}_{,i,t}={\displaystyle \frac{1}{5}}\sum _{k=-2}^2{Yr}_{i,t+k}^{obs}$$

(1)

where ${Yr}_{i,t+k}^{obs}$ represents the observed yield for years within the 5-year window.

The yield anomaly (ΔYi,t) quantifies the relative deviation of $Yr$ from the $Ydetrended$(Equation 2):

$${△Y}_{i,t}={\displaystyle \frac{{Yr}_{i,t}^{obs}-{Ydetrended}_{i,t}}{{Ydetrended}_{i,t}}}\times 100$$

(2)

Positive anomalies indicate yields above the expected trend, while negative anomalies indicate yield reductions relative to the local trend.

For each crop season, we calculated the negative yield anomalies or losses ($Yloss$) individually, estimating that for the soybean area at the municipal level. It was calculated based on the ${△Y}_{i,t}$and the $Ydetrended$calculated by five-year moving windows, according to Equation 3:

$${Yloss}_{i,t}=\left\{\begin{array}{c}-{\displaystyle \frac{{△Y}_{i,t}}{100}}\times {Ydetrended}_{i,t}\times A_{i,t}\\ 0\end{array}\right.if{△Y}_{i,t}<0,otherwise$$

(3)

where $A_{i,t}$is the soybean harvested area for each municipality $i$ in year $t$. This calculation ensures that only negative (less than zero) deviations contribute to losses, consistent with previous Brazilian soybean studies [11].

2.3. Weathear Dataset

Daily precipitation (P), maximum temperature (Tmax), minimum temperature (Tmin), relative humidity (RH), solar radiation (Rs), wind speed at 2 m height (u2), and reference evapotranspiration (ETo) were obtained from the Brazilian Daily Weather Gridded Data (BR-DWGD) dataset [33]. The dataset provides daily temporal resolution, a spatial resolution of 0.1° x 0.1°, and covers the period from 1989 to 2022 for the entire Brazilian territory. Z-score normalization was applied to standardize the scales of all climatic variables, facilitate multivariate analysis, and identify statistical outliers. This normalization reduced within-cluster variation which was essential for subsequent analytical procedures [34]. The normalized dataset was used to compute the climate indices and DHE metrics. Additionally, standardized climate indices, such as SPI, SPEI, and WSDI, were derived from the normalized BR-DWGD datasets.

2.4. Crop Evapotranspiration (ETc)

Crop coefficient (Kc) values were determined using the methodology proposed by [35] in FAO Irrigation and Drainage Paper 56. The FAO 56 procedure defines Kc as a function of the crop’s phenological stages, subdivided into initial, mid, and late-season phases. For each municipality, phase boundaries were aligned with the region-specific crop calendar, and fixed Kc values were assigned to each phase based on the mean coefficients recommended for soybean in FAO 56 and validated in recent Brazilian studies using lysimetry and remote sensing data [29,36,37].

Daily Crop Evapotranspiration (ETc) was calculated using the standard FAO 56 formulation [35] (Equation 4):

$${ET}_c=K_c\times {ET}_o$$

(4)

where ETo is the reference evapotranspiration, estimated using BR-DWGD for each municipality, and Kc is the stage-dependent crop coefficient assigned to each day based on the local phenological calendar [33]. This approach aligns with experimental and modeling studies that quantify soybean water use across Brazil using daily Kc curves and FAO 56 methodology [9,38].

2.5. Drought Indices

Meteorological drought was also characterized using the Standardized Precipitation Index (SPI), which is based solely on precipitation and quantifies precipitation anomalies across multiple accumulation time scales [20]. Monthly precipitation totals were aggregated over 1-, 3-, 6-, and 12-month periods to represent short- to long-term water deficits relevant to hydrometeorological and agricultural processes. For each municipality, SPI was computed independently to preserve local precipitation climatology. Accumulated precipitation series were fitted to a two-parameter Gamma distribution using maximum-likelihood estimation, following the original definition proposed by [20]. The Gamma distribution is appropriate for modeling nonnegative, right-skewed precipitation data and has been extensively validated for SPI applications across diverse climatic regimes [39]. Zero precipitation values were handled in accordance with standard SPI conventions using a mixed-probability adjustment. Cumulative probabilities derived from the fitted Gamma distribution were transformed into standardized normal variates, yielding SPI values with a mean of zero and unit variance. Negative SPI values indicate drier-than-normal conditions, while positive values indicate wetter- than-normal conditions, with larger magnitudes corresponding to greater anomaly intensity. SPI calculations were implemented in Python using a reference implementation consistent with the original methodology.

In addition to SPI, we used the Standardized Precipitation Evapotranspiration Index (SPEI), which extends the SPI framework by incorporating atmospheric evaporative demand through the climatic water balance [40]. The monthly climatic water balance was calculated as the difference between precipitation and potential evapotranspiration, and then accumulated across multiple temporal scales. SPEI was also calculated at two temporal scales, 3 and 6 months, to capture both short-term agricultural impacts and longer-term hydrological variability. The 6-month scale reflects water deficits during critical crop growth stages and is widely used to assess the effects of droughts on agricultural productivity, as documented in previous studies [40,41].

For each municipality and accumulation scale, the accumulated climatic water balance was fitted separately for each calendar month using a three-parameter Pearson Type III distribution. Monthly fitting ensures stationarity of the distributional parameters and enables consistent standardization across seasons, as recommended for drought index calculation [21,42]. The Pearson Type III distribution was selected to accommodate both negative and positive values of the water balance, preserving the physical interpretation of dry and wet periods. Cumulative probabilities from the fitted Pearson Type III distribution were transformed into standardized normal variates, yielding SPEI series with zero mean and unit variance. Negative SPEI values indicate drier-than-normal conditions, reflecting combined precipitation deficits and enhanced evaporative demand, whereas positive values indicate wetter-than-normal conditions. All SPEI calculations were performed using a custom Python implementation based on the standard SPEI methodology, adapted to use the Pearson Type III distribution for robust handling of the full range of water balance values [21,43]. Drought conditions were classified using SPEI values and standardized categories (Table 2), ranging from extreme drought (SPEI < -2.0) to non-drought or wet conditions (SPEI > 1.00).

Droughts were characterized using four key metrics derived from the SPEI time series, namely severity (cumulative SPEI during the event), duration (number of days the event persists), intensity (severity divided by duration), and peak (lowest SPEI value). These metrics provide a multidimensional understanding of drought behavior, encompassing magnitude, persistence, and extremity, which are essential for agricultural and climatological assessments [21,44].

2.6. Heatwave Index

The first step in calculating the WSDI was to determine the 95th percentile of daily maximum temperatures (Tmax) for a reference period, often using a 30-year baseline. For each day of the year, the 95th percentile (T95) was calculated from the historical temperature record (Equation 5) [45].

$$T_{95}\left(d\right)={Percentile}_{90}\left(T_{max}\right(d,referenceperiod\left)\right)$$

(5)

where T95(d) is the 95th percentile of the daily maximum temperature for day d, and Tmax (d, reference period) represents the historical maximum temperatures on day d over the reference period [48].

$$WSDI={\sum }_{i=1}^{Nspells}D{urationofwarmspell}_i$$

(6)

where Nspells is the number of warm spells in the year, and the duration of warm spelli is the number of consecutive days in that spell.

2.7. Spatial Processing

Spatial processing and extraction of SPI, SPEI and WSDI indices at the municipal level were conducted using Google Earth Engine (GEE) and Python. Multi-polygon geometries representing Brazilian municipal boundaries were used as masks to aggregate and summarize weather variables and indices. This approach enabled high-resolution, municipality-specific assessments of climate extremes and variability. The data underwent quality control and preprocessing, including handling missing values, treating outliers, and standardizing variables where necessary [46]. The methodological steps for drought detection and classification are summarized in the flowchart (Figure 2), which integrates threshold selection, temporal-scale specification, and event parameterization.

3. Results

The central thesis states that Brazil’s soybean production, the largest worldwide, has been affected by chronic and interannual drought and heatwaves, challenging its expansion beyond its current geographic limits. The implications of exposure to changing extreme events begin at the municipal level and extend across vast Brazilian biomes and agroclimatic regions, including the Cerrado, Caatinga, Pampa, Pantanal, Atlantic Forest, and Amazon. For example, the Cerrado region accounts for almost 50% of Brazil’s soybean production, and the Atlantic Forest and Amazon represent the geographic limits of soybean production in this study and the frontiers of its expansion and intensification in the years to come. The following sections provide evidence of the spatiotemporal attributions of soybean yield responses to droughts and heatwaves.

3.1. Climate Characterization

We grouped municipalities using Principal Component Analysis (PCA) based on estimated temperature and precipitation anomalies across Brazil to capture regional patterns in agroclimatic conditions and soybean productivity. As shown in Figure 3, the PCA identified three climate clusters: Central-West, Northeast, and South (Primavera de Leste, MT, and Cristalina, GO, were included in the Central-West and the South, respectively). These clusters coincide with the agroclimatic regions reported in [47].

Brazil’s tropical dry winter covers ~53% of MT’s and all TO’s surface in the Central-West region (Figure 1). Annual precipitation (mean ~1600 mm) and chronic winter dryness are evidenced by the fact that 88% of the annual precipitation falls during the growing season (Table 3). This hydroclimatic regime contributes to reduced yield variation in an agroclimatic region that also overlaps with the Cerrado, an area also considered the largest soybean agricultural expansion in Brazil [48,49]. Agricultural expansion frontiers in the region include the state of Matopiba, accounting for approximately 73 million hectares, consistent with recent findings by [50] and [51]. The conspicuous expansion of agriculture in the region evidences the growing climate risks to crops due to increasing hydroclimatic instability. Yet agricultural practices may have secured soybean yields as the consistency of ETc spatiotemporal variability across Central-West Brazil illustrates (Table 3 and Figure 4). The ET changes also shed light on how cropping systems use the limiting water and energy availability to support plants’ water demands as the growing season evolves. To support this process [52] suggests that energy is a limiting factor during the early stages of the growing season when moisture is present. Then, as water vapor deficits increase in later stages, water released to the atmosphere gradually limits soil water content and, in turn, ET. While the mechanisms described above are relevant to other regions with similar rainfall patterns, [53] found that vapor pressure and solar radiation are among the most important drivers of ETo. As noted by [54] and [24], high Tmax and low precipitation often coincide during the terminal stages of the crop cycle, creating conditions of elevated atmospheric demand (high ETc), exacerbating soil moisture deficits, and increasing the risk of crop failure in non-irrigated systems [55,56]. In the Central-West and Northeast regions, ETo values are similar, but during long dry spells, when water becomes a limiting factor, radiation offsets the vapor pressure deficit, reducing ET. Thus, during wet and dry periods, chronic or interannual precipitation and temperature maxima and minima alternate as drivers of ETc, highlighting crops’ exposure to local-to-regional hydrometeorological conditions associated with drought [55,56].

The municipalities in the Hot-Dry Northeast cluster in BA, PI, and MA are characterized by high mean annual temperatures (mean ~32°C) and low annual precipitation (mean ~ 1100 mm), along with an elevated annual evapotranspiration (mean ~1050 mm), reflecting semi-arid to sub-humid conditions [59]. The municipalities within the Hot-Dry Northeast cluster, which includes regions in Bahia (BA), Piauí (PI), and Maranhão (MA), experience high mean annual temperatures of around 32°C and low annual precipitation of approximately 1100 mm. Additionally, these areas have elevated annual evapotranspiration, averaging about 1050 mm, which indicates semi-arid to sub-humid conditions. Despite these conditions, similar intraseasonal variations in crop evapotranspiration between the Northeast and Central-West Brazil indicate their proximity to the tropics. Consequently, solar radiation levels are a major factor influencing patterns. During the growing season in the Northeast, summer precipitation accounts for 91-98% of total annual rainfall, the highest proportion observed in our study. In times of drought, summer precipitation becomes the key factor influencing agricultural outcomes in these arid areas, particularly due to the reliance on rainfed agriculture [60,61]. Further, similar intraseasonal variations between the Northeast and Central-West Brazil suggest that locations near tropical Brazil and, consequently, the associated radiation levels are the key drivers of these interdependent hydroclimate intraseasonal variations.

Overlapping with the South cluster, the Paraná is the second-largest area for soybean production in Brazil. This PCA cluster includes the states of Rio Grande do Sul (RS) and Paraná (PR), along with the Cristalina municipality, part of Brazil's humid subtropical region, characterized by hot to temperate summers. This area enjoys favorable agroclimatic conditions for soybean cultivation. In RS and PR, summer rainfall accounts for 64-69% of the annual precipitation, which totals approximately 1400 mm. These humid subtropical conditions, coupled with summer temperatures averaging around 23.5°C, contribute to soybean yields that often exceed 2400 kg/ha. Subject to low temperatures, frost risk in the South becomes a factor in August, when the expected warmer and wet conditions help successful sowing, the onset, and the evolution of the growing season across this area [9,46,54,62]. Thus, average precipitation and temperature magnitudes and fluctuations indicate that water deficits are not the primary constraint, since the autumn-winter months and mid-spring are typically wet and have low evaporative demand (Figure 4 and Figure 5). Temperature and precipitation values in the South show limited differences compared to those in the Central-West and Northeast. However, at municipal level, differences and similarities emerge suggesting that the sensitivity of municipalities such as Cristalina in the upper lands of CO and Primavera do Leste and Sorriso in MT to moisture or energy limited conditions can be due to changes in atmospheric evaporative demands and the inherent water pressure deficits driven by natural climate variability (Table 3) [60,61,62,63].

3.2. Interannual and Spatial Variability of Drought and Heatwave Events

The spatiotemporal variability of rainfall is a key determinant of rainfed soybean productivity in Brazil. The drivers of rainfall distribution over the soybean production areas suggest that a series of intertwined processes fueled by the contrasting sea surface temperatures (SST) between the Pacific and Atlantic Oceans, land surface-atmosphere interactions, and other climate-driven teleconnections, regulate the precipitation deficits and surpluses in Central-West, Northeast, and South Brazil [12,51]. Syntheses of recent research studies [64,65,66,67,68] supporting some of the variations of precipitation in soybean production regions indicate that the coupling between SSTs in the tropical Atlantic circulation patterns across the equator, which consequently drive contrasting wet and dry conditions in South America’s North and Northeast and Southeast. Furthermore, [69] indicate that contrasting SSTs now between the Atlantic and Pacific Oceans, during the warm and cold phases of the El Niño Southern Oscillation, lead to contrasting droughts, extended wet periods, and floods between northern and southern Brazil.

Analysis of rainfall and temperature event frequencies in selected municipalities within Brazil’s primary soybean-producing regions, as illustrated in Figure 6, reveals a “climate divide” between the North (MT, TO, BA, PI, and MA) and South (RS, PR, and GO). This divide corresponds to drought conditions in the North and Northeast clusters, which are closely linked to El Niño events, and to wet conditions in the Southeast that occur during the opposite phase of the El Niño Southern Oscillation (ENSO), specifically La Niña. Although these interannual climate variability modes offer important insights for crop and water management, intraseasonal changes remain uncertain, particularly when soybean planting, flowering, or grain filling are impacted by reduced residual soil moisture during the austral winter or by prolonged water deficits throughout the growing season [23,24,54,64,70,71].

3.2.1. SPI and SPEI Relationship to Soybean Yield Losses

A balance between precipitation and crop water requirements underscores how hydrometeorological and climate variability shape agricultural outputs. Even slight deviations from normal seasonal rainfall can stress soil moisture, directly influencing drought indices such as SPI and SPEI. These and other drought indices capture both precipitation anomalies and atmospheric evaporative demand, integrating climatic and water availability [20,21], and helping diagnose agricultural drought conditions and assess the sensitivity of yields to climate extremes [72]. Examining the temporal variability and spatial distribution of the SPI and SPEI in relation to soybean yield losses is essential for understanding how climatic fluctuations translate into agricultural vulnerability across the study regions.

Figure 7 shows the spatial distribution of SPI-6 for each growing season, with March as the reference month, coinciding with the end of the soybean cycle across most regions. The spatial patterns indicate a persistent drought signal in the North and Northeast, especially in the Northeast, where SPI-6 values frequently fall below -1.0. As mentioned in the previous section, ENSO has contrasting effects between the Northeast and Central-West regions and the South. During El Niño events (e.g., 2001–2003, 2015–2016), the North and Northeast typically experience negative precipitation anomalies, increasing drought risk, while the South generally receives above-average rainfall [73].

The cold phase (La Niña) produces the opposite effect of El Niño, resulting in reduced precipitation in the South and increased rainfall in the North and Northeast. These distinct drought drivers establish ENSO as a critical factor in risk management for southern soybean production. During El Niño years, more abundant and regular rainfall supports higher yields and mitigates the irregularity of frontal rains in the South [74]. In contrast, La Niña events significantly elevate the risk of yield losses in the South, as reduced and unpredictable precipitation may coincide with critical stages of the soybean growth cycle, intensifying water deficits and severely impacting potential yields [62]. Additionally, municipalities in the South and Southeast exhibit moderate SPI-6 variability, indicating alternating wet and dry years. Although these regions have historically experienced relatively stable rainfall patterns, there has been an increase in moderate-to-severe drought episodes since the early 2000s. This trend is indicative of regional drying and increased evaporative demand, consistent with a shift in the South Atlantic Convergence Zone (SACZ) and a reduction in convective rainfall [70,75,76].

Figure 8 and Figure 9 present the frequency distributions of drought classes derived from the 12-month SPI and SPEI, respectively, along with their relationships with observed soybean yield losses across the studied municipalities. The analysis demonstrates a consistent association between severe-to-extreme drought conditions (SPI/SPEI ≤ −1.5) and significant yield reductions, particularly in municipalities in the Northeast. Under these conditions, yield losses often exceed 20% relative to the municipal trend-adjusted mean, emphasizing the strong dependence of soybean productivity on sufficient hydrologic balance throughout the crop cycle.

In contrast, municipalities in the South experienced a higher frequency of near-normal and moderately wet classes (−0.5 < SPI/SPEI < 0.5), resulting in lower yield losses (<10%). This spatial differentiation highlights the gradient in drought exposure across the soybean frontier and is consistent with previous findings of northward intensification of hydroclimatic risk [73,75]. Southern municipalities such as Cascavel and Palmeira das Missões exhibited robust yield trends but also persistent drought signals and rising temperatures, demonstrating that technological advances alone are insufficient to counteract intensifying climatic extremes. This observation supports earlier conclusions by [77], which found that technology can reduce but not fully offset the impacts of increasing climate variability on crop yields.

Comparison of SPI and SPEI distributions indicates that SPEI identified a greater proportion of drought-affected months in most municipalities, especially during the 2015–2016 and 2019–2020 growing seasons. The discrepancy between indices suggests that evaporative demand, rather than precipitation alone, was the primary driver of water deficits during these years, consistent with regional warming trends and increased atmospheric aridity [76]. The stronger correspondence between SPEI-defined droughts and yield anomalies further supports the consensus that SPEI provides a more comprehensive representation of agricultural drought under current climatic conditions [21].

3.2.2. The Warm Spell Duration Index (WSDI)

The Warm Spell Duration Index (WSDI), defined as the annual count of days with at least six consecutive days exceeding the 95th percentile of daily maximum temperature, serves as a robust indicator of heat events affecting agroecosystems [25]. WSDI is essential for assessing thermal stress episodes that directly influence crop phenology, water demand, and yield potential [70,78]. Analysis indicates a statistically significant upward trend in WSDI across all surveyed municipalities from 1996 to 2020 (p < 0.05), consistent with broader regional warming trends reported for Central-West Brazil [70]. Northeast municipalities exhibited the largest increases in warm spell duration, with average annual increases exceeding one day per decade (Figure 10). These prolonged warm spells intensify drought conditions by increasing crop evapotranspiration demand and accelerating soil moisture depletion, thereby compounding the effects of heat and water deficits on soybean productivity [49,73]. The observed WSDI trends are consistent with projections of increased frequency and intensity of heatwaves in tropical agricultural zones under future climate scenarios [4].

Compound drought and heat events are increasingly responsible for yield losses, and the alignment of the Standardized Precipitation Evapotranspiration Index (SPEI) and WSDI indicates that combined stress effects exceed those of individual factors [49,79]. The 2015–2016 El Niño event exemplified this phenomenon, resulting in yield reductions exceeding 30% across municipalities in Northeast and Central-West Brazil [73,80].

3.3. Quantification of Soybean Yield Losses at the Municipal Level

Soybean yield losses at the municipal level in Brazil demonstrated significant spatial and temporal variability over the 30 crop seasons analyzed. During the first decade (1989/1990–1999/2000), total yield losses amounted to 252.2 thousand tonnes (kt), representing an average reduction of 3.0% across municipalities. The most affected municipalities in this period included Sorriso (74.9 kt; 2.6%) in the Central-West and Barreiras (66.0 kt; 4.3%) in the Northeast. While losses were concentrated in the northern portion of the study area, several municipalities, such as Bom Jesus and Lagoa da Confusão in the Northeast and mountainous regions in Central Brazil, did not experience drought or yield losses during the 1990s (Table 4). The spatial distribution of these municipalities suggests that local convective events and other teleconnections play a critical role in modulating the regional impacts of ENSO and the ITCZ (Costa et al., 2024). This pattern is also observed in regions where both sides of the “climate divide” previously experienced drought and associated yield losses, including Palmeira das Missões (278.1 kt; 13.8%) in the South and Barreiras (227.0 kt; 6.7%) in the Northeast. As shown in Table 4, subsequent decades reveal both contrasting and consistent yield responses to drought.

In a broader international context, these losses are comparable to drought-related soybean yield reductions reported in major producing regions, including the US Midwest and the Argentine Pampas, where decadal-scale yield penalties of 3 to 7 percent during periods of climate stress have been documented. These results suggest that Brazilian municipalities are exposed to climate phenomena of similar intensity as those affecting other global soybean-producing regions [81,82].

During the most recent decade, total soybean yield loss reached 1618.1 kt (3.7%), maintaining a similar relative magnitude to the preceding period but exhibiting different spatial distribution patterns. Over the entire 30-year period, municipalities experienced cumulative losses of 3292.3 kt, corresponding to an average of 3.7% of potential production.

Combined Climate Drivers of Droughts in Soybean Regions

Table 5 and Figure 11 demonstrate a clear trend of increasing soybean production stress across Brazil’s principal producing regions, primarily due to climate-induced water deficits. Yields have risen markedly in municipalities within the eastern states included in the study (BA, MA, and GO), reflecting advancements in technology and management practices [73,83]. Nevertheless, negative SPEI-6 anomalies and positive WSDI trends reveal that hydroclimatic stressors are intensifying, especially in municipalities in the Northeast (MA, TO, PI, and BI). In the latter half of the 2000s, negative SPEI-6 values became more persistent, coinciding with higher mean temperature anomalies and extended warm spells. This pattern indicates a transition from precipitation-driven to temperature- and evapotranspiration-driven droughts [75,76].

Positive WSDI trends, particularly in the Northeast, indicate longer heatwave durations during the soybean season, consistent with projections of increased heat stress under warming climate scenarios [75]. The inconsistent spatial and temporal patterns of SPI-6 across most sites underscore the complexity of precipitation variability and its implications for agricultural productivity. These findings suggest that, despite improvements in yield, climate stressors such as drought and heat may increasingly limit soybean production in rapidly expanding agricultural regions, highlighting the necessity for climate-adaptive strategies. In municipalities located in the Central-West and along the Cerrado region, robust and consistent trends imply that adaptation measures, including irrigation, enhanced genetics, and improved management practices, have achieved relative success (Table 5). Nonetheless, the significant WSDI and SPI trends observed in these areas indicate that even well-adapted systems remain susceptible to heat-driven evapotranspiration losses, especially under dry-season planting schedules [84].

In the Northeast, the concurrence of negative SPEI-6 values with prolonged heat spells (high WSDI) intensifies crop water stress beyond the effects of precipitation anomalies alone. Elevated temperatures accelerate evapotranspiration and reduce soil moisture availability during critical phenological phases [49,70]. The widening gap observed in these municipalities indicates increasing climatic constraints on productivity, particularly where irrigation capacity is limited and evaporative demand is high. In contrast, municipalities in the Cerrado region exhibit smaller negative SPEI-6 trends and only moderate yield reductions despite rising WSDI. This pattern suggests that adaptive capacity, including cultivar selection, management practices, and technological investment, can partially mitigate heat-related impacts [85,86].

Table 5 and Figure 11 further demonstrate that rising heat stress and water deficits interact with technological yield gains, resulting in heterogeneous effects across regions. In areas experiencing rapid agricultural expansion, such as the Northeast, the increasing frequency of compound drought–heat events may constrain future yield potential. This trend highlights the necessity for climate-adaptive strategies, especially in rainfed systems where thermal extremes and evapotranspiration-driven drought are becoming primary limitations.

4. Discussion

The findings demonstrate that high-resolution, spatially disaggregated gridded climate datasets, when validated at the municipal scale, accurately represent the spatial and temporal variability of soybean yield responses to DHE stress across Brazil’s principal soybean macroregions. These results substantiate the assertion that standardized hydroclimatic indices (SPI, SPEI, and WSDI) effectively characterize regional yield sensitivities to DHE. Integrated trend analyses of these indices and detrended yield data reveal that all studied municipalities experienced significant declines in water availability and increased heat stress during the soybean growing season from 1989 to 2020. The Northeast region exhibits the most pronounced yield losses, corresponding with recurrent rainfall deficits and intensifying temperature anomalies, while the Central-West displays emerging negative yield trends associated with prolonged heatwave durations despite technological advancements. In the South and Central-West highlands, previously marked by greater climatic stability, precipitation variability is now increasing, disrupting crop development, particularly during reproductive stages. These results further indicate that the compound effects of DHE are undermining yield stability across all CZ, aligning with projections of tropical and subtropical agricultural vulnerability [4,76].

The recent intensification of DHE has exerted a greater influence on soybean production trends over the past decade, as indicated by precipitation anomalies. Additionally, detrended yield analyses demonstrate that negative yield deviations increasingly coincide with periods of simultaneous negative SPI and SPEI anomalies and positive WSDI trends, confirming that concurrent heat and moisture stress dominate interannual variability [49,73]. Seasonal analyses identify water deficits during the reproductive and grain-filling stages as primary drivers of yield loss, resulting from insufficient recovery following early-season rainfall deficits. This seasonal mismatch between water demand and availability amplifies reproductive stress and accelerates phenological progression, ultimately reducing grain mass and final yield. The integration of high-resolution climate data, standardized indices, and yield detrending represents a key strength of this study, facilitating the identification of underlying yield–climate mechanisms independent of technological progress. Nevertheless, the absence of management and socioeconomic variables constrains the ability to attribute yield trends exclusively to climatic factors. Technological and agronomic adaptations, such as irrigation scheduling and optimized planting calendars, partially mitigate climatic impacts but remain inadequate to reverse the persistent downward yield trend under intensifying compound extremes.

Yields across the three clusters remain similar (Figure 3e), indicating that soybean varieties and agricultural management practices have adapted to relatively stable intraseasonal climate variability in Brazil’s soybean production areas (Figure 4 and Figure 5). However, the sufficiency of these adaptations in addressing the increasing frequency, duration, and intensity of droughts remains uncertain [87]. The compounded nature of drought events presents significant challenges to water and food security in regions where soybean agriculture is expected to expand further. Soybean responses to anomalously dry seasons or years are often triggered by heatwaves and deficiencies in agronomic and water management. [70] propose several strategies to mitigate drought impacts in Brazil, including diversification of agriculture, selection of drought-resistant varieties, and the use of native species for forage and wood. These strategies are closely linked to plant exposure to climate variations during the growing season. Thus, the exploration of changes in atmospheric evapotranspirative demand (AED), vapor pressure deficits, and ETc is critical for adaptive soybean cropping systems, particularly during unexpected dry years, and involves complex interactions between natural and human environments [63,88,89]. AED influences evapotranspiration and the processes governing available energy in areas where land use, vegetation, and soil moisture are continuously changing. Therefore, scientific and technological advancements aimed at mitigating rising temperatures and associated droughts and heatwaves in regions such as Paraná and Cerrado should prioritize constraining the positive feedback loop, in which crop evapotranspiration and soil moisture decline as temperatures continue to increase.

In summary, the analysis confirms that DHE exerts an increasing influence on soybean yield trajectories across Brazil’s soybean regions. The combined application of SPI, SPEI, and WSDI demonstrates their complementary diagnostic capabilities in identifying both precipitation- and temperature-driven components of agricultural variability. These indices establish an empirical foundation for operational climate-risk surveillance and adaptive management within compound-risk frameworks, particularly in tropical cropping systems where warming intensifies atmospheric water demand. Future agricultural expansion or intensification will benefit from studies that integrate process-based crop models with high-frequency meteorological data, management records, and remote sensing to refine yield response functions under DHE, as well as socioeconomic indicators to enhance vulnerability diagnostics and inform resilient agricultural planning [90,91,92]. Furthermore, integrating genomic, climate, and management data to model phenotypes will improve the capacity to nowcast, forecast, and predict phenotypes across regions, accounting for inherent variability and exposure to extreme events [93,94,95]. These advancements are essential for developing predictive models that capture the multidimensional and interactive nature of climate risks affecting regional soybean systems in Brazil.

5. Conclusion

This study provides a spatiotemporal assessment of compound drought and heat extremes affecting Brazil’s main soybean producing regions from 1989 to 2020. By integrating high resolution meteorological data, yield records, and standardized hydroclimatic indices SPI, SPEI, and WSDI, the analysis demonstrates that climate stressors are intensifying and becoming increasingly compound in nature. Distinct CZ exhibit contrasting sensitivities, with southern regions showing relatively higher yield stability, while northeastern and frontier regions experience greater concurrence of rainfall deficits and heatwaves during critical phenological stages, leading to greater yield variability.

The results show that DHE has increased in frequency and intensity, particularly during reproductive phases, amplifying stress on rainfed soybean systems. Drought related indices, especially SPI-6, exhibit strong correlations with yield variability in vulnerable regions such as Maranhão, Piauí, and Bahia, underscoring their relevance for early warning systems and climate risk monitoring. Temperature driven increases in evapotranspiration and prolonged heatwaves further exacerbate drought severity, placing the Northeast as the most climate sensitive soybean frontier, while the Central-West and South regions face emerging risks linked to heat extremes and precipitation variability despite technological advances.

Overall, the findings underscore the need to revise agroclimatic risk frameworks to explicitly account for compound hydroclimatic extremes. Regionally differentiated adaptation strategies, supported by climate informed planning, improved monitoring using compound indices, and investments in resilient cultivars and infrastructure, are essential to sustain soybean production amid rising climatic volatility. Future research should integrate finer scale agronomic, genetics and phenological, and socioeconomic data, along with remote sensing and crop modeling, to better capture multidimensional vulnerability and support targeted adaptation policies.