The Pantanal under Siege—On the Origin, Dynamics and Forecast of the Megadrought Severely Affecting the Largest Wetland in the World

The Pantanal is the largest wetland of the world and one of the most important biodiversity hotspots in South America. An unprecedented ongoing megadrought is severely affecting its ecological functioning, flood pulse dynamics, and fire regime. Regarding this problematic, the present study generates reliable information about the following key issues: 1—Evolution and dynamics, 2—Origin and determinants, and 3—Forecast based on identified determinants and current trends. Results show that the evolution of the megadrought has been differentiable in both, space and time. As for its origin and determinants, Climate Change was ratified as one of the most important threats to the Pantanal, and to vast areas of South America, since a strong correlation was identified between megadrought’s dynamics and the occurrence of intense marine heatwaves at Northern Hemisphere oceanic waters, and more specifically, at the Northeast Pacific. Results also show that the megadrought is expected to continue at both the Pantanal and the surrounding Highlands, at least until December 2023. Thus, an intensification of fires risk, extending now to areas historically flooded or perhumid should be expected, concomitantly to a very negative impact on non-fire-resistant vegetation cover, as well as ecosystem functioning and biodiversity, perhaps even worse than those from 2020, widely covered by the international media.


Introduction
The Pantanal is the largest wetland in the world [1][2][3]. It was designated a World Heritage Site by UNESCO in 2000. Besides a high number of resident animal species, the Pantanal is also very important biome for the migratory ones, making this wetland one of the most important biodiversity hotspot in South America [4,5]. The Pantanal functions as a large reservoir that collects water from the surrounding Highlands during the rainy season and then gradually delivers it to the lower sections of the Paraguay River resulting in characteristic and historically predictable seasonal flooding or flood pulse [6,7]. Here, human activities mostly rely on ecosystem services, and include cattle farming, agriculture, Although wide reference was achieved about the direct effects on environmental and societal issues, besides that provided by Thielen et al. [13], little evidence was available, to the media, to stakeholders and to policy-makers, regarding key issues such as the dynamics, origin, and forecast of the extreme drought that already is considered to be the most prolonged and most severe in the last 60 years [21], and that has already severely affected the Pantanal [20]. Which, as forecasted by Thielen et al. [13], and according to preliminary observations from Thielen et al. [26], may be part of an ongoing global extreme oceanic event (i.e., an extreme MHW) generating an unprecedented climatic response (i.e., a megadrought) that might be far from over.
Thus, the main objective of this work was to generate reliable information about three key issues regarding the megadrought severely affecting the largest wetland in the world. First, to characterize its spatial and temporal dynamics along the different regions and subregions from the Pantanal and surrounding areas, and since its very beginning until present. Second, to evaluate recent oceanic SST dynamics and identify the oceanic region that best explains current megadrought in terms of its duration and intensity. And third, to generate a reliable forecast based on identified determinants and current trends. The combined use of remote sensed data from different sources, analyzed through sophisticated spatial analysis tools, allowed this study to generate reliable information about the evolution of current megadrought has been differentiable in both, space and time. An evolution that was not explained by the SST dynamics from neighboring oceanic regions, but rather from much far distanced regions such as the Northeast Pacific. From which, the presence of a very intense and prolonged warm SST gives solid hints to expect as to forecast for the megadrought to last much longer and to be even more pervasive.

Study Area
The study area consisted of regions and sub-regions of the Upper Paraguay River Basin (UPRB) which forms part of the upper La Plata River Basin [28] within west-central Brazil. With a total extension of 361,338 km 2 , it includes the Pantanal (138,183 km 2 ), a run-on drainage area occupying the center of the UPRB, and the Highlands (223,155 km 2 , see Figure 1), which consists of the surrounding watershed areas at (or above) elevations of 200 m asl.

Data
Precipitation data was obtained from the Climate Hazards group Infrared Precipitation with Stations (CHIRPS), freely available at https://iridl.ldeo.columbia.edu/ SOURCES/.UCSB/.CHIRPS/.v2p0/.monthly/.global/ (accessed on 10 October 2021). CHIRPS V2.0 is a quasi-global gridded rainfall time series dataset, spanning 50 • S-50 • N, from 1981 to near-present, 0.05 • resolution satellite imagery with in situ station data [30], with great applications monitoring precipitation extremes [13,26]. According to Beck et al. [31], in a global-scale evaluation of 23 precipitation datasets, CHIRPS V2.0 tends to perform the best in hydrological modeling of tropical regions, specifically in Central America, and in central and eastern Brazil. Thielen et al. [13], as well as Marengo et al. [20], tested its applicability and potentiality in the study of precipitation spatial-temporal dynamics at the UPRB. For the present study, monthly data for the time series January 1981/June 2021 was obtained from 486 rasters. Monthly and annual means, as well as some other basic precipitation parameters, were obtained by using GIS applications.
Quantifying global trends and variability in SST is of fundamental importance to understanding changes in the Earth's climate. One approach to observing SST is via remote sensing [32]. In this sense, the monthly mean SST remote sensed data was provided by the NOAA Extended Reconstructed Sea Surface Temperature Version 5 dataset (ERSSTv5), which is a global dataset derived from the International Comprehensive Ocean-Atmosphere Dataset (ICOADS). It is produced on a 2 × 2-degree grid with spatial completeness enhanced using statistical methods (http://iridl.ldeo.columbia.edu/maproom/Global/ Ocean_Temp/Monthly_Temp.html, (accessed on 10 October 2021)). For the present study, monthly data for the time series January 1982/June 2021 was obtained from 474 rasters.

Calculation of the Standardized Pluviometric Drought Index-SPDI
In this study, precipitation spatial-temporal dynamics was analyzed by the Standardized Pluviometric Drought Index (SPDI) developed by Pita-López [33]. The SPDI is a monthly rainfall index that is based on the calculation of cumulative monthly rainfall anomalies, similar to the well-known Standardized Precipitation Index (SPI) of McKee et al. [34], more specifically, the 12-month SPI. As in this index, values ranging from +1 to +1.5 and +1.5 to +2.0 are associated with moderately humid and very humid episodes, respectively, and values exceeding +2 are representative of extremely humid episodes. Moderately dry, very dry and extremely dry spells are characterized by the same ranges with a negative sign (see Table 2). The SPDI is calculated as follows: First stage: where APi is the monthly precipitation anomaly, Pi is the monthly precipitation, and P MED is the median precipitation of the month for the series. As for this study: 1981-2010. Second stage: where APAi is the accumulated precipitation anomaly of the month. Third Stage: where APA is the average value of accumulated precipitation anomalies of all the months of the series, and σAPA is the standard deviation of accumulated precipitation anomalies of all the months of the series. GIS applications allowed to implement these equations to the 486 aforementioned CHIRPS V2.0 rasters and generate SPDI products such as images of 0.05 • resolution or monthly zonal values, and at different space and/or time criteria. In the present study, monthly values of SPDI were estimated based on the climatological series 1981-2010. From this information, most recent spatial-temporal precipitation anomaly dynamics regarding current megadrought was identified for the UPRB, as well as for the Pantanal and the Highlands, and their respective sub-regions (see Figure 1). Significance of statistical difference between sets of monthly precipitation and/or SPDI values was identified by two-tailed Paired t-Test (α = 0.01). While similarities in spatial-temporal dynamics of SPDI between the different sub-regions were identified by Cluster Analysis (K-means clustering using Euclidean distance). Gong and Richman [35] noted that nonhierarchical methods, such as the K-means algorithm, outperformed hierarchical methods (the Ward's method and the average linkage method) when tested with precipitation data, as well as for SPI series [36]. Table 2. Categories resulting from SPDI estimation, adapted from McKee et al. [34] and WMO [37].

SPDI Values Category
≥2.0 Extremely wet 1.5-2.0 Severely wet 1.0-1. 5 Moderately wet −1.0-1.  Figure 2). According to Thielen et al. [13], these specific oceanic regions show the strongest (positive or negative) correlation with historical precipitation dynamics at the UPRB (series 1981-2018). Additionally, traditional SST indexes such as TNA-Tropical North Atlantic and TSA-Tropical South Atlantic, as well as NEP-Northeast Pacific, and the ENSO regions: NIÑO 1+2, NIÑO 3 and NIÑO 3.4, described by PSL/NOAA [38] and CPC/NCEP/NOAA [39], were considered in the present study. Mann-Kendall tests were applied to determine the significance of the time series trends in the monthly SSTs for each oceanic zone.  A Marine Heatwave (MHW), which is an anomalously warm water event that may last many months and extend over thousands of square kilometers, was identified and categorized following Hobday et al. [17]. According to these authors, the different categories of MHWs are based on multiples of the difference between the climatological mean (January 1982/December 2011) and the climatological 90th percentile, where: Moderate (1-2×, Category I), Strong (2-3×, Category II), Severe (3-4×, Category III), and Extreme (>4×, Category IV).

Calculation of Correlation of Precipitation with the Oceanic SSTs
GIS analytical tools were used to investigate correlations between mean monthly SST and mean monthly precipitation in the Pantanal and the surrounding Highlands. The application of correlation analysis permits the identification of regions where the interactions between hydrological droughts and teleconnection are strong [40]. Results from Section 2.3, that is from the precipitation spatial-temporal dynamics, were used to identify the series that best describes current megadrought. While results from Section 2.4 provided hints about the oceanic region with most significant SST variability. Analysis assumed that the observed SST preceded precipitation by lead times of 0 to 6 months. The Pearson correlation coefficient was used to estimate the degree of explained variance (α = 0.01) between the monthly SST and precipitation. A paired, two-tailed t-Test was used to assess the significance of statistical differences. Confidence intervals of 99% were used to identify correlation coefficients (r) ranges. Close attention was paid to resulting correlations that would best explain recent precipitation dynamics (2019-present) at the Pantanal and the surrounding Highlands.
A scheme of the workflow for the methodology of the present research is provided in Figure 3.  shows that the Pantanal is significantly drier than surrounding Highlands (p < 0.01). Even though the Highlands showed annual precipitation amounts to be about 19% higher than in the Pantanal (1454 vs. 1176 mm/year, respectively), rainfall here was less variable: over 90% of months during this series had SPDI near-normal ((−1.0 ≤ SPDI < 1.0). Until the beginning of 2020, precipitations events at the Highlands capable of generating abnormal SPDI values, either for dry or a wet condition, were rather short-lived and moderately intense (see Figure 4). Since June 2019 onwards, the monthly precipitation amounts suffered a reduction of about 30% from historical mean. Right from January 2020, precipitation anomaly in the Highlands passed from moderate drought (−1.5 ≤ SPDI < −1.0) to a severe drought condition (−2.0 ≤ SPDI < −1.5), reaching a sustained extreme drought condition (SPDI ≤ −2.0) from September 2020 onwards.  [34] and WMO [37] (see Table 2). The grey area represents the near-normal condition (i.e., −1.0 ≤ SPDI < 1.0).
As for the Pantanal, during the series January 1981-June 2021, about 78% of the months showed a near-normal precipitation condition. Here, anomalies in precipitation, dry or wet, occurred in longer pulses than in surrounding Highlands. In the Pantanal, dry spells have been significantly more frequent than wet ones (p < 0.01; 15.9 vs. 6.4% of the months). Three main dry spells have occurred in the Pantanal in the last 40 years (see Figure 4). The first one is a moderate drought lasting 15 months (September 1993/November 1994) with a mean SPDI value of −1.48. A second dry spell of 14 months (January 2013/February 2014) generated a severe dry anomaly in the Pantanal (mean SPDI = −1.68), and during three months (from March/May 2013), precipitation anomalies generated an extreme drought condition (SPDI ≤ −2.0). None of these two dry pulses affecting the Pantanal generated precipitation anomalies in the surrounding Highlands. As for the third dry spell, and like it had occurred in the neighboring higher lands, since June 2019 onwards, the monthly precipitation amounts at the lowlands of the Pantanal underwent an important reduction from the historical mean (of about 32%). As for resulting SPDI values, the drought condition was reached four months' sooner in the Pantanal than in the Highlands (September 2019 and January 2020, respectively). The evolution of this particular dry spell has been most dramatic in the Pantanal-in the term of only six months, the precipitation anomalies passed from a near-normal condition to a sustained extreme drought condition (SPDI ≤ −2.0).
Such evolution is also appreciated in Figure 5 which represents the temporal dynamics of the SPDI monthly values from January 2019 to June 2021, together with mean historic and observed precipitations, at the Pantanal and the surrounding Highlands. As for the Pantanal, from June 2019 onward, 83.3% of the months had precipitations below the historical values. This value was 87.5% in the case of the Highlands. Such persistent precipitation anomaly generated an extreme drought that passed to affect initially around 5.7% of the Pantanal area in June 2019, to 41.5% in January 2020 ( Figure 5a). By March 2020, over 90% of the Pantanal was affected by a drought condition (SPDI < −1.0), and over 60% by an extreme drought (SPDI ≤ −2.0). The extreme drought condition not only has persisted for a long time in the Pantanal, but it has expanded very importantly: 82.2% by June 2021. Even though the number of abnormally dry months generating the extreme and extended drought in the Pantanal was not significantly different (p > 0.05) from that occurred in the Highlands, here, the drought had a more discrete dynamic (Figure 5b). It is only by September 2020, six months later than in the lower lands, that 50.7% of the Highlands area was affected by the extreme drought condition. As for most recent update, June 2021, 70.5% of the Highlands was affected by a drought condition, from which 56.1% by an extreme drought.

Spatial Dynamics
Spatial dynamics of the drought (January 2019 to June 2021) in the UPRB is given in Figure 6. Here, differences between the Pantanal and the surrounding Highlands, as well as in their different sub-regions, are better observed. According to this figure, the drought dynamics along the series is spatially differentiable. First evidence of a drought developing was observed in the sub-region of Paraguai (P-08) in central section of the Pantanal, showing a severe drought right from the beginning of 2019. Along the first year, and specifically by the second half, an easterly spread of the drought occurred. First as moderate (−1.5 ≤ SPDI < −1.0) and reaching the severe condition (−2 ≤ SPDI < −1.5) by the end of 2019. By the end of the first trimester of 2020, and from there on, the extreme drought condition (SPDI ≤ −2.0) was generalized in the middle and most of the northern section of the Pantanal and corresponding Highlands. Only the southernmost sectors of the UPRB, and a small section of the Highlands (H-11) draining to Cáceres (P-11), showed near-normal precipitations during the entire series. Accurate quantification of such dynamics is given in Figure 7. For instance, the central Pantanal sub-region of Paraguai (P-08), from the beginning of the series until October 2019 reached a mean SPDI of −1.84 (named Group 1, according to Figure 7a [13]. Cluster analysis (K-means clustering using Euclidean distance) was performed to both rows and columns, and with the statistical tool ClustVis [41].
During the entire series, that is from June 2019 until June 2021, there are three subregions that did not show any precipitation anomaly (mean SPDI = 0.08): Porto Murtinho (P-01) in the Pantanal, and in the neighboring Highland the sub-regions of H-01-02 and H-03; all three located in the southern most section of the UPRB (Group 3, according to Figure 7a,b). There is also a fourth group of sub-regions regarding their latter and/or more spatial heterogeneous response in precipitation anomaly dynamics (Group 4, Figure 7a,b), comprising the Pantanal sub-regions of Nabileque (P-02), Miranda (P-03), and Abobral (P-05), and by the two Highland sub-regions of H-04 and H-11). Here, drought conditions did not start until March or April 2020. Although for Abobral (P-05) and H-11, the drought was from severe to extreme for the rest of the series (mean SPDI = −1.94 and −3.82, respectively), drought conditions in the other three sub-regions were less intense (moderate) and less spread (heterogeneous), and even intermittent, as in the case of the southern highland sub-region of H-04.  Figure 2). It was established that most significant SST dynamics (p < 0.01) resulted from the warming of oceanic waters in the Northern Hemisphere, more specifically in the three regions already identified by Thielen et al. [13], such as: PAC-NE, PAC-NW and ATL-N, with resulting mean SSTA of 1.28, 0.77 and 0.50 • C, respectively. In absolute terms, resulting SSTA at PAC-NE was from 3-6× higher than that for any of the ENSO regions, and 6-9× higher than for any referential SST indexes such as TNA and TSA. As for NEP, and even tough 81% of PAC-NE area is comprised in this northeast oceanic region, resulting NEP's mean SSTA is significantly lower than that of PAC-NE (p < 0.01; 1.04 vs. 1.28 • C, respectively). As for PAC-NE, during the series June 2019/June 2021, 88% of the area presented monthly mean SSTA between 1.0 and 2.0 • C, with a maximum of 2.21 • C in August 2019. While for NEP, anomalies of such intensities involved a lesser area: 51%, reaching a maximum of 2.02 • C also by August 2019. PAC-NW, on the other side, had a lesser extended and warm SST: only 23% of the area involved SSTA >1.0 • C, and 59% involved SSTA between 0.5 and 1.0 • C. On another ocean, the Atlantic, such mild warming involved 50% of ATL-N. The rest of the area of this oceanic region had a slight warmth in its SST of less than 0.5 • C.  Figure 2). Now, besides such spatial SST dynamics, it is precisely in the northern oceanic regions of PAC-NE, NEP, PAC-NW and ATL-N, where the most significant warming trend was observed (p < 0.01, series January 1982/June 2021). Among these, the oceanic regions from the Northeast Pacific, PAC-NE and NEP, showed the strongest trend ( Figure 9). An important contribution for such warming trend comes from two very distinctive SSTA positive pulses. The first one is a very extended SST warm pulse occurring from April 2013 until October 2016, lasting 43 months. The warming in this first pulse was less variable and significantly higher in NEP than in PAC-NE (1.02 vs. 0.84 • C, respectively; p < 0.01). As for the second Northeastern Pacific SST warming pulse, it started around April 2019 and, as for June 2021, still goes on (SSTA June 2021 = 1.35 and 0.84 • C, for PAC-NE and NEP, respectively), lasting until present, over 27 months. Differently than the first warm pulse, current pulse is significantly warmer at PAC-NE than at NEP (1.19 vs. 0.99 • C, respectively; p < 0.01). As for the Northwest Pacific, PAC-NW, this oceanic region has shown a very extended but variable warm SST pulse. With a mean of 0.69 • C, it was significantly less warm (p < 0.05) than Northeastern counterparts. The warm pulse in PAC-NW started in June 2017, way before current megadrought, and has lasted almost continuously until present (49 months). A brief pause in the anomalously positive condition in PAC-NW occurred in May 2019, little after PAC-NE and NEP started their second warming pulse (i.e., by April 2019). As for the Atlantic Ocean, the warming did not seem to respond to very extended intense SST warm pulses. For instance, most recent SST dynamic for ATL-N shows a rather mild pulse starting in January 2019 with a mean SST of only 0.48 • C. From Figure 9 it is also evident that SST warming dynamics in the long term has been generally significant for TNA (p < 0.01). But, such dynamics, differently than the other aforementioned oceanic regions, is rather discrete regarding most recent SST data: a warmth of only 0.36 • C was evident from January 2020 until March 2021.

Marine Heatwave (MHW) Dynamics at the Northeast Pacific
Another important output from the temporal analysis of SST dynamics at the Northeast Pacific, besides its intensity, is regarding its duration. About this, Figure 10 shows that, for over two continuous years now (since June 2019/June 2021), SSTA at PAC-NE has been, not only well above the historic mean (climatological series January 1982/December 2011), but almost constantly above the 90th percentile. According to Hobday et al. [17], the combination of SSTA of such intensity and such persistence defines very well a Marine Heatwave (MHW) condition. In the case of PAC-NE current warming, the MHW condition has oscillated from the categories of moderate to about severe (Figure 10a). The MHW condition was also reached at NEP about the same time than in PAC-NE (by June 2019), but in a significantly (p < 0.01) lesser extent. As a matter of fact, while the MHW condition has persisted in PAC-NE until present (June 2021), at NEP it experienced a three-month break after February 2021 (Figure 10b).

Marine Heatwave (MHW) Dynamics at the Northeast Pacific
Another important output from the temporal analysis of SST dynamics at the Northeast Pacific, besides its intensity, is regarding its duration. About this, Figure 10 shows that, for over two continuous years now (since June 2019/June 2021), SSTA at PAC-NE has been, not only well above the historic mean (climatological series January 1982/December 2011), but almost constantly above the 90th percentile. According to Hobday et al. [17], the combination of SSTA of such intensity and such persistence defines very well a Marine Heatwave (MHW) condition. In the case of PAC-NE current warming, the MHW condition has oscillated from the categories of moderate to about severe (Figure 10a). The MHW condition was also reached at NEP about the same time than in PAC-NE (by June 2019), but in a significantly (p < 0.01) lesser extent. As a matter of fact, while the MHW condition has persisted in PAC-NE until present (June 2021), at NEP it experienced a three-month break after February 2021 (Figure 10b).  Figure 2). MHW categorization is taken from Hobday et al. [17].
As mentioned, two major warm SST pulses have occurred in the Northeast Pacific since January 1982 (see Figure 9). Such pulses generated two important and extended MHWs in both, PAC-NE and NEP (see Figure 11). No other oceanic regions, from those considered in the present study, showed MHW conditions of such intensity and/or duration. As for NEP, these two extreme warm pulses were identified in the present study as Blob 1 and Blob 2. Because 80% of PAC-NE is encompassed in NEP, for the purpose of the present study, the Blob 1 and Blob 2 are SST phenomena simultaneously occurring in PAC-NE and NEP (Figure 11a,b). gorization is taken from Hobday et al. [17].
As mentioned, two major warm SST pulses have occurred in the Northeast Pacific since January 1982 (see Figure 9). Such pulses generated two important and extended MHWs in both, PAC-NE and NEP (see Figure 11). No other oceanic regions, from those considered in the present study, showed MHW conditions of such intensity and/or duration. As for NEP, these two extreme warm pulses were identified in the present study as Blob 1 and Blob 2. Because 80% of PAC-NE is encompassed in NEP, for the purpose of the present study, the Blob 1 and Blob 2 are SST phenomena simultaneously occurring in PAC-NE and NEP (Figure 11a,b).
Regarding Blob 1 dynamics, MHW condition during this extreme oceanic event was more persistent in NEP than in PAC-NE: Four breaks of two or three months at PAC-NE, versus only one two-month break occurred at NEP. As for Blob 2, although it is still an ongoing phenomenon, its MHW dynamics has been inversely that of Blob 1: more persistent in PAC-NE than in NEP, showing only a brief one-month break from the MHW condition.  [13]) and (b) the region NEP-Northeast Pacific SST index (from CPC/NCEP/NOAA [39]) (see Figure 2). Solid blue line represents the historic mean SST for the climatological series January 1982-December 2011. As for the solid black line, it represents the 90th percentile, while the dashed red line represents the observed monthly SST variability. Red areas represent the moments when MHW condition has been reached (i.e., SSTA > 90th percentile). As for intensity, the thicker the resulting (red) area, the more intense the MHW. Grey areas are relative to the occurrence of Blob 1 and Blob 2.   Figure 2). Solid blue line represents the historic mean SST for the climatological series January 1982-December 2011. As for the solid black line, it represents the 90th percentile, while the dashed red line represents the observed monthly SST variability. Red areas represent the moments when MHW condition has been reached (i.e., SSTA > 90th percentile). As for intensity, the thicker the resulting (red) area, the more intense the MHW. Grey areas are relative to the occurrence of Blob 1 and Blob 2.

Blob
Regarding Blob 1 dynamics, MHW condition during this extreme oceanic event was more persistent in NEP than in PAC-NE: Four breaks of two or three months at PAC-NE, versus only one two-month break occurred at NEP. As for Blob 2, although it is still an ongoing phenomenon, its MHW dynamics has been inversely that of Blob 1: more persistent in PAC-NE than in NEP, showing only a brief one-month break from the MHW condition.

Cross Correlating Monthly Precipitations with SST Dynamics
Besides expressing the highest SSTA, as well as showing the most persistent MHW condition, SST dynamics from PAC-NE correlates the best with precipitations anomalies generating current megadrought at the UPRB. Results from cross correlating for the series June 2019/June 2021 show that, with a lead time of two months, SST at the PAC-NE significantly correlate (p < 0.01) with precipitations, explaining the ongoing megadrought condition: about 73.5% (for an r = −0.8574) in the case of the Pantanal, and 84.5% (for an r = −0.9194) for surrounding Highlands. Therefore, the warming trend observed in SST at PAC-NE implies a decrease in UPRB precipitation.

Megadrought's Forecast
From historical precipitation anomalies (series January 1981/June 2021) at the Pantanal and the surrounding Highlands (Figure 12a) and the application of time series forecasting techniques such as the Exponential Triple Smoothing (ETS), it is possible to make a projection of the megadrought. As shown in Figure 12b, negative precipitation anomalies (as SPDI) are expected at least until December 2023. Thus, current megadrought is expected to persist in the UPRB, severely affecting both, the Pantanal and the Highlands. In this time series forecasting, no significant difference was identified between observed vs. forecasted precipitation anomalies, in neither the Pantanal nor the Highlands (t-Test, p = 0.3367 and 0.6101, respectively).

Megadrought's Evolution
The evolution of current megadrought has been differentiable in both, space and time. First evidence of a drought developing was observed right from the beginning of 2019 in the central section of the Pantanal. By mid-2019, an easterly spread of the drought occurred. And, in the term of nine months (by April 2020), over 60% of the area of the Pantanal, as well as 30% of the surrounding Highlands, succumbed to an extreme dry condition. Such extreme condition has persisted as a megadrought and, as for most recent update (June 2021), spreading over 80% of the Pantanal and about 60% of the Highlands. It is appropriate to classify the extreme event as a megadrought for this being very lengthy and pervasive, as well as lasting much longer than any previous drought in the area [42] (see Figures 5-7 and 12).
In their analysis, which involves exclusively precipitations occurring at the Pantanal and only a segment of the megadrought (up to November 2020), Marengo et al. [20], referring to the initial dynamics of the drought, suggest a late onset of the rainy season (Oct-Mar) of the hydrological year 2019-2020 with a rainfall deficit of 40%. In the present study, such deficit was prominently lower (23%). A deficit which was not sensibly different from that generated during the next rainy season of 2020/2021 (17%). Now, the origin of the differences between the results from both studies may lay in the extension of the series involved.
At the Pantanal, dry spells have been significantly more frequent than wet ones. Three main dry spells have occurred in the Pantanal in the last 40 years. The first one is a moderate drought lasting 15 uninterrupted months (September 1993/November 1994), while the second dry spell lasted 14 continuous months (January 2013/February 2014) and generated a severe drought in the Pantanal. Although historically dry spells have not affected the surrounding Highlands [13], is it during the third and last dry spellthe current megadrought-that negative precipitation anomalies not only affected the Pantanal, but the surrounding Highlands as well. Nevertheless, the megadrought temporal dynamics along the UPRB was differentiable-the drought condition was reached four months' sooner in the lowlands of the Pantanal than in the higher lands of the surrounding Cerrado (September 2019 and January 2020, respectively).

SSTA Dynamics and Resulting MHWs
Since mid-2019, and about simultaneously to the dramatic spread of the megadrought in the UPRB (with a lead of two months), a very important warming of Northern Hemisphere oceanic waters has been in process. More specifically, in the oceanic regions of PAC-NE, PAC-NW and ATL-N, previously identified by Thielen et al. [13]. Such warming has been most significantly intense and extended in waters of the Northeast Pacific, involving both, PAC-NE, and NEP. No other oceanic region, from those considered in the present study, showed SSTA dynamics, for warming (or even for cooling) as intense and/or as enduring as that observed in the Northeast Pacific since mid-2019.
Such dynamics of SSTA is precisely what describes a MHW-A discrete period of prolonged anomalously warm water at a particular location, exceeding a fixed [43], seasonally varying [44], or cumulative [45] threshold. The MHW condition is specifically reached when SSTA not only surpasses the historical mean but above the 90th percentile [17]. Short lasting MHWs (over few days) may disturb the marine environment, it is the long-lasting events that have the most significant impact [46]. With only a short one-month break at the end of 2019, PAC-NE has constantly reached the MHW condition since June 2019. The MHW condition has been significantly higher and more persistent in PAC-NE than in NEP. Such persistence of a MHW at the Northeast Pacific was also referred by the reports from Global Ocean Monitoring and Observing Program (GOMO), ratifying the MHW condition for the Northeast Pacific in its report of 12 July 2021 [39].
Results also show that observed SSTA dynamics and resulting MHWs at the Northeast Pacific is part of a generalized warming trend of the Northern Hemisphere oceanic waters. While for the Atlantic Ocean such trend has been rather gradual since 1982, for the Pacific this has been more accentuated from the beginning of the 21st century (as for PAC-NW), and most extreme this last decade for eastern side of the North Pacific (as for PAC-NE and NEP). In the case of Northeast Pacific, an important contribution for such warming trend came from two very distinctive MHWs. The first one, is a very extended SST warm pulse occurring from April 2013 until October 2016, lasting 43 months. Originally known by the moniker of "The Blob" in the general press [47]-as Blob 1, in the present research. As for the second northeastern Pacific MHW, Blob 2 (after Amaya et al. [48]), it started around April 2019 and still goes on lasting, as for June 2021, over 27 months. Up to this point, the warming dynamics during Blob 1 and Blob 2 has been distinctive, not only between these two important MHWs, but between the oceanic regions of PAC-NE and NEP. For instance, during Blob 1, SSTA was less variable and significantly higher in NEP than in PAC-NE, while the opposite has been true during the ongoing Blob 2. Blob 1 is the largest and longest MHW ever recorded between 1981 to 2017, not only for NEP, but in general [49][50][51]. Now, because Blob 2 is an ongoing process and, as for June 2021, there is no sign of a soon end, much from this important MHW event is still be expected.
The participation of the different climate modes modulating the frequency, intensity and duration of the different MHWs is still under great debate [44,52]. The Pacific Decadal Oscillation (PDO), the North Pacific Gyre Oscillation (NPGO), El Niño-Southern-Oscillation (ENSO), among others, influence ocean temperatures and thus the development of MHWs, directly or remotely via atmospheric or oceanic teleconnections, which reverberate the effects globally [53]. As for potential drivers of Blob 1, there has been identified a link to the dynamics of the two dominant modes of winter sea surface temperature variability in the North Pacific, the Pacific Decadal Oscillation (PDO), and the North Pacific Gyre Oscillation (NPGO) [52,54]. Specifically, between winter warm NPGO anomalies and the following winter PDO arising from extratropical/tropical teleconnections. But, unlike Blob 1, Blob 2 peaked in the summer, a season when little is known about the physical drivers of such events [48]. According to these authors, Blob 2 primarily results from a prolonged weakening of the North Pacific High-Pressure System, which reduces surface winds and decreases evaporative cooling and wind-driven upper ocean mixing. Warmer ocean conditions then reduce low-cloud fraction, reinforcing the MHW through a positive low-cloud feedback. Amaya et al. [48] also found that remote SST forcing from the central equatorial and, surprisingly, the subtropical North Pacific Ocean contribute to the weakened North Pacific High. This study also highlighted the possibility that, due to resulting Northeast Pacific mixed layer depths extremely shallow, Blob 2 would not last much longer than winter and spring 2020 seasons. Results from present study show that, even though the MHW condition at PAC-NE and NEP did weaken and even disappeared for a couple of months during that time, as for PAC-NE, the MHW condition was not only reestablished and has persisted until present, with no evidence of weakening. This result reinforces the importance not only to identify the drivers capable of initiating MHWs at different seasons [48], but also about the necessity of evaluate potential alternating dynamics, intra e inter annually, among the different drivers capable of determining the life-span and intensity of important foreseeable MHWs occurring at the Northeast Pacific.

Cross Correlations between Monthly Precipitations and SST Dynamics
Even though the exact mechanisms for how SST may influence precipitation at a local level was not considered in this study, the application of correlation analysis on different hydrological indices and SST data permits the identification of regions where the interactions between hydrological droughts and teleconnection are strong [40]. As a consequence of aforementioned Blob 2 (i.e., an extended SST warm pulse occurring at the Northeast Pacific), a very persistent and intense MHW condition was generated at the oceanic region PAC-NE and the associated NEP right from June 2019, which is the same moment identified in the present study as the beginning of current megadrought. The results from lagged correlation between the monthly mean precipitations at the Pantanal and the surrounding Highlands with monthly SST at the different oceanic regions analysis show that, for the series June 2019/June 2021 and with a two-month offset (lead time), it is PAC-NE's MHW condition that best correlates and explains, right from the beginning of the megadrought, the precipitation anomalies at the UPRB: 73.5% in the case of the Pantanal, and 84.5% for surrounding Highlands.
As for the rest of the oceanic regions considered in the present study: from the tropics (i.e., PAC-SAC-Pacific South American Coast, and the ENSO regions: NIÑO 1+2, NIÑO 3 and NIÑO 3.4), as well as from the Southern Hemisphere (i.e., ATL-SC-Atlantic South Central, ATL-S-Atlantic South, PAC-S-Pacific South, and IND-S-Indian South), and even referential SST indexes such as TNA-Tropical North Atlantic and TSA-Tropical South Atlantic; all these failed explain the megadrought dynamics as significantly as did the Northeast Pacific region of PAC-NE. In these terms, a warming trend in SST at PAC-NE and the concomitant presence of a MHW implies a decrease in UPRB precipitations, with the possibility of generating a megadrought condition. Nevertheless, because the existence of a significant correlation between the SST of the other two Northern Hemisphere oceanic regions-PAC-NW in the West Pacific and ATL-N in North Atlantic-and the historical dynamics in the precipitations of the Pantanal [13], there is also the possibility of future megadroughts resulting from current warming trend observed in these two specific oceanic areas.
Evaluating the potential causes of the precipitation deficit from the austral summer (December 2019/February 2020) at the Pantanal, Marengo et al. [20] propose that it was caused by a reduction in the transport of warm and humid (austral) summer air from Amazonia into the Pantanal. These authors suggest that such event occurred in association with an anomalous warming of the Atlantic Ocean captured by the AMO and TNA indices, inducing an anomalously northward position of the Atlantic ITCZ, producing less rainfall in southern Amazonia and the Pantanal. But, Marengo et al. [20] also found that other drought episodes in the Pantanal not necessarily follow this pattern. Regarding this possibility, the outcomes from the present study, which include analysis from the second and most recent austral summer from current megadrought (December 2020/February 2021), show that, due to the generalized but significant warming trend observed in the North Atlantic, as for oceanic regions such as ATL-N and TNA, SSTA for last austral summer (December 2020/February 2021) did was higher than for the preceding (December 2019/February 2020). But, regardless of such warmth difference (+0.18 • C and +0.38 • C, in the case of TNA), no significant concomitant response in precipitation anomaly was detected among the two austral summers, as would be expected.
Regarding major SST based indices such as TNA, AMO, PDO, and ONI, Marengo et al. [20] identified a very low correlation between rainfall in the Pantanal, as well as Paraguay River level at Ladário and SST dynamics. The results from present research, not only agrees with such findings, but provides information about the specific areas, encompassed in the very extended ocean sections where these and other macroclimate indices express their SST dynamics at different time scales, from interannual to interdecadal, and expressing the highest correlation with precipitation dynamics at the UPRB, not only at different time scales, but also at different space scales. About this last issue, Batista Silva et al. [55] established that, with a lead of 3 months and explaining 66% of the monthly variance, PDO negative phases concomitantly occurred with drier periods in the Pantanal. Now, according to CPC/NCEP/NOAA [39], as for June 2021, the PDO phase is negative and has persisted so since January 2020. Thus, current negative phase of PDO fails, by itself, to fully explain current megadrought dynamics, right from its beginning in June 2019. Now, results from present study show that, when using in this kind of correlation analysis a more specific and targeted oceanic region, as resulted to be PAC-NE explaining precipitation variability at the UPRB, the lead time for response shortens (2 months), as well as sensibly improves the capacity to explain (74-85%) any resulting drought, both in time and in space, and all along the evolution, not only of current megadrought, but of other recent important drought events.
One of such events occurred during 2013 and part of 2014. South America was experiencing an important drought at about the same time of the so-called "second dry spell" of the Pantanal, described earlier in this study (see Section 4.1). Now, according to Rodrigues et al. [56,57], such dry condition was generated by an unprecedented MHW developed in a rather small and specific area of the western South Atlantic that lasted for almost all of the summer (<4 months), which apparently suppressed the South Atlantic convergence zone and its associated rainfall. Regarding this specific MHW, from the three South Atlantic related oceanic regions considered in the present study (i.e., ATL-SC, ATL-S and TSA), only ATL-S showed a moderate MHW condition, lasting only 3 months (February/April 2014). As for the specific South Atlantic oceanic region indicated by Rodrigues et al. [56], resulting SSTA and MHW dynamics do not explain satisfactorily the drought in the Pantanal since, by the time the warming of SST occurs and the short-lived MHW reaches its maximum (by February 2014), the severe dry condition was already affecting the Pantanal for over a year (since February 2013). Results from the present study regarding the presence of warm SST at the Northeast Pacific and the presence of Blob 2 and the concomitant generation and dynamics of a MHW at PAC-NE, provide the solid arguments as to explain this particular 2013/2014 dry spell from the concurrent presence of Blob 1 at the Northeast Pacific and the generation of a persistent and intense MHW at NEP. An explanation that points out to the most probable origin of such specific dry spell in the Pantanal, as might explain also, the origin of a more extended climatic phenomenon, such as the generalized drought concomitantly affecting vast areas of South America. Results from present study leave solid open grounds as for the (re)evaluation the origin and dynamics, not only extreme droughts that have affected extended regions along history, but also more recent and even ongoing extreme dry events such as the drought that is currently spreading and intensifying across South America [58].

Megadrought's Forecast
The application of time series forecasting techniques to historical precipitation allowed to make a projection of the ongoing megadrought affecting the UPRB. Results from present study show that the dry condition (as negative precipitation anomalies) is expected at both the Pantanal and the surrounding Highlands at least until December 2023. Even if precipitation amounts normalize from this point on, it would take over 30 continuous months as for any drought condition to disappear in the Pantanal area. Now, the sooner the current extreme drought condition can become about normal is if, at both the Pantanal and the surrounding Highlands, each month of next rainy season (i.e., October 2021/March 2022) reaches a relative positive anomaly of 43% or above from historical mean (series 1981-2020). A positive anomaly of such magnitude has never been reached in the Upper Paraguay river basin at least for the last 40 years. For instance, the wettest rainy season (i.e., October 2016/March 2017) presented an anomaly of 21% above historical mean, that is less than half the amount of rain currently demanded.
Even though fire is considered a natural disturbance in the Neotropical savannas, in the case of the Pantanal, due to its wetter condition, fires here have been historically less frequent and spatially restricted [26,59]. However, given the ongoing megadrought, reductions in precipitation are still expected to cause significant disturbances in its ecological functioning, affecting hydrological, floodplain inundation dynamics, as well as fire regime. The intensification of fire risk, extending now to areas historically flooded or perhumid, is expected. Concomitantly, and as predicted by Thielen et al. [26], it is likely a negative impact on non-fire-resistant vegetation cover, as well as ecosystem functioning and biodiversity, perhaps even worse than those from 2020. Now, according to the latest report from the Global Water Monitor and Forecast Watch List (15 July 2021, [58]) which provides a 12-month prediction through March 2022 indicate that exceptional deficits are forecasted not only along the Paraguay River through its namesake, but also along the Paraná River through Brazil, Paraguay, and Argentina. As mentioned in the previous section, the exceptional drought is becoming a generalized condition in many parts of South America, a condition that is expected to keep spreading. From preliminary estimations, we found that 22% of South America is currently under a drought as severe as the one currently affecting the Pantanal. With 82% of its area, Pantanal is the biome most affected by the megadrought, followed by the Mata Atlântica (55%), the Cerrado (37%), Caatinga (16%) and Amazônia (9%). As for Brazil, 35% of its territory is currently under different degrees of droughts, from which 23% is extreme. According to Parris [58], widespread water deficits ranging from severe to exceptional are expected in the states of Mato Grosso, Mato Grosso do Sul, São Paulo, and Paraná. Deficits reaching exceptional intensity are forecast for parts of western Amazonas and Acre, Pará and Amapá, and Maranhão. Regarding SSTA forecasts, the Climate Forecast System of the National Centers for Environmental Predictions [39] recognized that initial projections underestimated the strength of current Northeast Pacific MHW and suggested that the MHW condition in this oceanic region should continue for the next 6 months.
As mentioned, the projection of the present study is founded on consideration to historical variability, status and trends of precipitation and the SSTs from several oceanic regions, specific to the series June 2019-June 2021 and based in the climatological January 1981-December 2011. Any case, further drought projections should be reevaluated in terms of the newly identified causal connections between the Pantanal dry pulses and the increasing Northern Hemisphere SST and foreseeable increase in occurrence MHWs [13,26]. This adjustment would allow to determine whether predicted conditions are consistent and supported by current data and models or whether these conditions indicate larger scales and more complex climate dynamics. The correlations between precipitations at the UPRB with more than one oceanic region stresses the importance of future studies evaluating combined influence [13,20].
Regarding MHWs, even though knowledge about their functioning, drivers, impacts and future evolution, is still in its infancy [44,52,60], consensus is emerging that anthropogenic climate change has significantly increased the likelihood of recent MHWs, especially for the more extreme category events [15,16,18,51]. Frölicher et al. [43] found that MHW duration doubled from 1982 to 2016 and attributed 87% of MHWs to human-induced warming. Significant increases in MHW intensity and count of annual MHW days are projected to accelerate, with many parts of the ocean reaching a near-permanent MHW state by the late 21st century [18].
Improved skill to predict oceanographic conditions in the North Pacific is highly desirable [61]. Blob 1's warmth penetrated hundreds of meters below the surface, but Blob 2's has not reached such depth yet [48]. To answer how long will Blob 2 last is a matter of knowing how deep its warmth will penetrate below the surface. This kind of information would improve drought forecasting capacity for the Pantanal, and even for extended regions of South America, since both Blob 1 and Blob 2 have proven capacity of generating extreme MHW conditions in PAC-NE as well as in NEP, an oceanic phenomenon that is strongly correlated with the occurrence and dynamics of extreme droughts in the Pantanal and surrounding areas. Now, in the last quarter of this decade, global ocean heat gain has increased in the upper 700 m of the ocean and heat has been sequestered in deeper ocean layers at depths down to more than 2000 m [62]. Ocean Heat Content refers to the heat absorbed by the ocean. More than 90% of the excess heat is stored within the world's oceans, where it accumulates and causes increases in ocean temperature [63,64]. Knowing how much heat energy is stored in the ocean-and where it is stored and released-is essential for understanding the state, variability, and changes of Earth's climate system [62,65]. The gradual buildup of subsurface heat content throughout 2020 in the region suggests the potential for persistent ecological impacts [66]. As warmer temperatures gradually mix downward, they can persist long after the surface MHW disappears, suggesting that the ocean can provide memory for long-lived MHWs [66]. The stronger ocean warming within upper layers versus deep water has caused an increase of ocean stratification in the past half century [67]. With increased stratification, heat from climate warming less effectively penetrates into the deep ocean, which contributes to further surface warming [68].
As for 2019, the world's oceans (especially the upper 2000 m) were the warmest in recorded human history [69]. Due to resulting teleconnections, this is one of the key reasons why the Earth experienced such catastrophic fires in the Amazon, California, and Australia in 2019 [65]. An update for this list would have to include the unprecedented fires in the Pantanal, for 2020, and upcoming years until the precipitation siege from MHW at PAC-NE ceases. It is important to note that ocean warming will continue even if the global mean surface air temperature can be stabilized at or below 2 • C (the key policy target of the Paris Agreement) in the 21st century [64,65,70]. The ocean, and some other components in the Earth system, are slow to respond and equilibrate, and will continue to change even after radiative forcing stabilizes [64]. However, the rates and magnitudes of ocean warming, and the associated risks will be smaller with lower levels of anthropogenic GHG emissions [64,65,70].

Conclusions
The present study generates reliable information about three key issues: 1-Megadrought's evolution and dynamics, 2-Megadrought's origin and determinants, and 3-A forecast based on identified determinants and current trends.

1.
The evolution of the megadrought currently affecting the Pantanal has been differentiable in both, space and time. First evidence of a drought developing was observed in the central section of the Pantanal, right from the beginning of 2019. By mid-2019, an easterly spread of the drought occurred, and as for June 2021, over 80% of the Pantanal area and about 60% of the Highlands is affected by the megadrought. As for temporal dynamics, the drought condition was reached four months' sooner in the lowlands of the Pantanal than in the surrounding higher lands.

2.
The study ratifies that Climate Change is one of the most important threats to the Pantanal. The evolution and dynamics of current megadrought is correlated to a significant warming trend of the Northeast Pacific sea surface waters (SST) which resulted, since June 2019, in one of the most intense and persistent marine heatwave (MHW) from recent times, identified by other authors as Blob 2. The present study also provides evidence that Blob 1, a preceding MHW also from the Northeast Pacific strongly correlates to the second most important drought (2013-2016) occurred in the Pantanal area in the last 40 years. Encompassed in the Northeast Pacific, the SST warming of PAC-NE region explains, right from the beginning of the megadrought, 73.5% of the precipitation anomalies in the case of the Pantanal, and 84.5% for surrounding Highlands. As for South American neighboring oceanic regions (e.g., TNA, TSA and the ENSO regions), their SST dynamics was rather discrete, and did not explain Pantanal's megadrought in the terms as did the presence of a MHW in the Northeast Pacific. Now, the SST warming trend at the Northeast Pacific is part of a generalized warming trend of the Northern Hemisphere oceanic waters. There is a historical correlation between the presence of warm SST at PAC-NW in the Northwest Pacific and ATL-N in the North Atlantic and the occurrence of dry pulses in the Pantanal area.

3.
Pervasive MHWs from warming oceanic waters of the Northern Hemisphere-The "new climatic reality" that the Pantanal and other vast regions of South America have to overcome. Results from this study as well as from others provide arguments as to consider the problematic of warmer SST as a new climatic reality. Anthropogenic climate change has significantly increased the likelihood of recent MHWs, especially for the more extreme category events. Blob 2 is an ongoing process with no signs to end any soon. Blob 1's warmth penetrated hundreds of meters below the surface, but Blob 2's has not reach such depth yet. To answer how long will Blob 2 last, ergo the megadrought, is a matter of knowing how deep will its warmth penetrate below the surface. In addition, and because the existence of a significant correlation between the SST dynamics of the other two Northern Hemisphere oceanic regions (i.e., PAC-NW and ATL-N) and the historical dynamics in the precipitations of the Pantanal, there is also the possibility of future megadroughts resulting from current warming trend and potential MHWs generating in these two oceanic areas. Regarding a forecast for the megadrought, results from present study show that it is expected to continue at both the Pantanal and the surrounding Highlands at least until December 2023. Even if precipitation amounts normalize from this point on, it would take over 30 continuous months as for any drought condition to disappear in the Pantanal area. The sooner the current extreme drought condition can become about normal is if, at both the Pantanal and the surrounding Highlands, each month of next rainy season (i.e., October 2021/March 2022) reaches a relative positive anomaly of 43% or above from historical mean (series 1981-2020). But a positive anomaly of such magnitude has never been reached in the Upper Paraguay river basin at least in the last 40 years. Due to the continuity of the megadrought, the intensification of fires with unprecedented duration and intensity, extending now to areas historically flooded or perhumid are expected. Concomitantly, we predict a most definite negative impact on non-fire-resistant vegetation cover, as well as ecosystem functioning and biodiversity, perhaps even worse than those from 2020. Now, the exceptional drought is becoming a generalized condition in many parts of South America. From preliminary own estimations, about 22% of South America is currently under a drought as severe as the megadrought affecting the Pantanal. With 82% of its area, Pantanal is the biome most affected by the megadrought, followed by the Mata Atlântica (55%), the Cerrado (37%), Caatinga (16%) and Amazônia (9%). Here, the exceptional drought condition is expected to keep spreading.
The present study provides valuable and reliable skills as to deal with key issues about such abrupt and extreme drought, providing to decision-makers, stakeholders and society most useful and opportune information regarding its origin, evolution and a forecast based on identified determinants and current trends. Any case, further drought projections should be reevaluated in terms of the newly identified causal connections between the dry pulse's dynamics and the increasing Northern Hemisphere SST and foreseeable increase in occurrence MHWs.  Data Availability Statement: Precipitation data was freely provided by the Climate Hazards group Infrared Precipitation with Stations at https://iridl.ldeo.columbia.edu/SOURCES/.UCSB/.CHIRPS/ .v2p0/.monthly/.global/ (accessed on 10 October 2021). As for the monthly SST data, this was freely provided by the NOAA at http://iridl.ldeo.columbia.edu/maproom/Global/Ocean_Temp/ Monthly_Temp.html (accessed on 10 October 2021).