Natural variability vs forced signal in the 2015–2019 Central American drought

The recent multi-year 2015–2019 drought after a multi-decadal drying trend over Central America raises the question of whether anthropogenic climate change (ACC) played a role in exacerbating these events. While the occurrence of the 2015–2019 drought in Central America has been asserted to be associated with ACC, we lack an assessment of natural vs anthropogenic contributions. Here, we use five different large ensembles—including high-resolution ensembles (i.e., 0.5∘ horizontally)—to estimate the contribution of ACC to the probability of occurrence of the 2015–2019 event and the recent multi-decadal trend. The comparison of ensembles forced with natural and natural plus anthropogenic forcing suggests that the recent 40-year trend is likely associated with internal climate variability. However, the 2015–2019 rainfall deficit has been made more likely by ACC. The synthesis of the results from model ensembles supports the notion of a significant increase, by a factor of four, over the last century for the 2015–2019 meteorological drought to occur because of ACC. All the model results further suggest that, under intermediate and high emission scenarios, the likelihood of similar drought events will continue to increase substantially over the next decades.


Introduction
Most of Central America-here loosely defined as the continental land between southern México and Panamá-experienced below average rainfall between 2015 and 2019 (Fig. 1a,  Fig. S1). This prolonged rainfall deficit resulted in severe and continued drought conditions separated by the midsummer drought (MSD; Magaña et al. 1999;Gamble et al. 2008, see also Fig. S1). In areas of Central America on the side of the Caribbean Sea, the precipitation distribution is more uniform during the wet season though, with higher annual accumulations compared to the rest of Central America (Alfaro 2002;Martinez et al. 2019). Because of the highly seasonal character of Central American rainfall, even small changes in rainfall can have a substantial impact on regional water resources and local rain-fed agriculture. Indeed, crop shortages due to extreme drought conditions over the last decade and the subsequent impoverishment have been identified as the main driver of mass migrations from Central America in recent years (World Food Programme 2020). Climate migration from Central America is expected to continue and ramp up in the next decades due to stronger impacts of global warming combined with rapid population growth (Rigaud et al. 2018;Hidalgo and Alfaro 2012).
Observational evidence of more prolonged MSDs along with more dry days has been recently discussed in the literature (Maurer et al. 2017;Anderson et al. 2019), but there has been little progress so far in attributing these changes to ACC. Herrera et al. (2018) examined the 2013-2016 Pan-Caribbean drought, which affected mostly the Caribbean islands, and concluded that ACC played a major role in increasing the drought's severity through the effects of higher temperature on potential evapotranspiration. This made its attribution to ACC somewhat less problematic, given the observed, long-term warming of the region (Aguilar et al. 2005;Pachauri and Meyer 2014). Attribution of a meteorological drought to ACC is instead more troublesome, since this is determined by local and large-scale circulation changes which, in turn, are generally less directly controlled by temperature (NAS 2016). We have no evidence of specific attribution studies connecting ACC to the recent 2015-2019 meteorological drought in Central America. Neelin et al. (2006) analyzed the recent multi-decadal summer drying trend over Central America and concluded that the observed trend cannot be unambiguously attributed to ACC, that is, it cannot be excluded that it is due to natural interdecadal variability.
This study has two objectives. First, to determine if the recent drying trend observed from the late 1970s is attributable to ACC. Second, to evaluate whether the 2015-2019 Central American meteorological drought has been influenced, i.e., made more likely, by ACC. These two questions will be addressed from a meteorological point of view by looking at the impact of climate change on the precipitation deficit. We realize that this is only the first step in analyzing the drivers of a complex phenomenon like a drought, which is not only caused by large scale circulation anomalies but also by complex interactions of local precipitation, wind, temperature, land use, and water management (Mishra and Singh 2010;AghaKouchak et al. 2021).
We will use five different large ensembles-including two high-resolution (i.e., 0.5 • horizontally) ensembles-to estimate the contribution of ACC and distinguish it from that of internal variability. Large ensembles of climate models (Deser et al. 2020) are very useful for climate risk and attribution studies (Otto et al. 2018;Swain et al. 2018;Pascale et al. 2020) as they provide thousands of years of data and thus allow for a direct reconstruction of the underlying probability distribution of hydroclimatic extremes without relying on a hypothesized statistical model of extremes (Van der Wiel et al. 2019). Central American summer rainfall has large magnitude interannual and decadal variability (Giannini et al. 2000(Giannini et al. , 2001Wang 2007;Hastenrath and Polzin 2013;Muñoz-Jimènez et al. 2019;Anderson et al. 2019) and a large ensemble is, thus, also a powerful method to isolate, at the decadal timescale, internal variability from the forced signal (Kay et al. 2015;Deser et al. 2020).

Observation data
To take into account uncertainty in observed precipitation over Central America, four different observational precipitation datasets are used in this study: (1) the Global Precipitation Climatology Centre (GPCC) dataset version 7, at 0.5 • horizontal resolution (Schneider et al. 2013), (2) the Climate Research Unit high-horizontal resolution grids of monthly rainfall at the University of East Anglia, version 3.24, at 0.5 • horizontal resolution (Harris et al. 2013), (3) the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS), at 0.05 • horizontal resolution (Funk et al. 2015), and (4) the CPC Merged Analysis of Precipitation (CMAP), at 2.5 • (Xie and Arkin 1997) horizontal resolution.
For atmospheric variables, like, e.g., the 925 hPa winds, we employ the ERA5 reanalysis (Hersbach et al. 2019). ERA5 is the most recent reanalysis product from the European Centre for Medium-Range Weather Forecast and it features major improvements such as higher horizontal and vertical resolution (0.28 • and 139 vertical levels), an updated and enhanced version of the Integrated Forecast System (IFS) Earth system model and associated observational assimilation system (Hersbach and Dee 2016).
Sea surface temperatures (SSTs) are taken from the Extended Reconstructed Sea Surface Temperature version 5 (ERSSTv5) dataset (Huang et al. 2017), which is a global monthly SST dataset from 1854 to present on a 2 • × 2 • grid derived from the International Comprehensive Ocean-Atmosphere Dataset.

Modelling data
We use a suite of large ensemble simulations from the the Seamless System for Prediction and EArth System Research (SPEAR; Delworth et al., 2020). SPEAR is the newest modeling system for seasonal to multi-decadal prediction developed at NOAA Geophysical Fluid Dynamics Laboratory (GFDL) and shares underlying component models with the CM4 climate model (Held et al. 2020). The SPEAR atmospheric model is run at different horizontal resolutions (atmosphere, land) in this paper: 0.5 • (SPEAR MED) and 1 • (SPEAR LO), and it has 33 atmospheric levels in the vertical. More details about SPEAR and the SPEAR large ensemble can be found in Delworth et al. (2020), Lu et al. (2020), Pascale et al. (2020), andMurakami et al. (2020). The SPEAR MED large ensemble is characterized by a horizontal grid-spacing that is finer than those of most other available large ensembles, which makes the SPEAR MED ensemble an unprecedented and unique tool to study regional climates.
To evaluate forced vs. natural variability, we use four different numerical experiments: (1) CTRL, a long-term control run with constant preindustrial (1850) forcing; (2) NATU-RAL, an ensemble driven only by natural forcing (i.e., volcanic eruptions and solar cycles); (3) ALLFORC4.5, an ensemble driven by observed anthropogenic and natural forcing up to 2014 (HIST), and then according to the Shared Socioeconomic Pathway (SSP2-4.5) developed for the Coupled Model Intercomparison Project Phase 6 (CMIP6) Eyring et al. 2016;Riahi et al. 2017); and (4) ALLFORC8.5, an ensemble driven by observed anthropogenic and natural forcing up to 2014 (HIST) and then according to the Shared Socioeconomic Pathway (SSP5-8.5). The SSP5-8.5 ("high emissions") represents the high end of the range of future scenarios (radiative forcing by 2100 of 8.5 W m 2 ) and it updates CMIP5 RCP8.5 ). SSP5-8.5 is compatible with a future development heavily based on fossil fuels. The SSP2-4.5 ("middle of the road emissions") instead sits in the middle of the range of future forcing pathways and updates the CMP5 RCP4.5 .
We also analyze four additional large ensembles to assess model uncertainty: SPEAR LO, the Forecast-Oriented Low Ocean Resolution model with flux adjustment, FLOR FA (Zhang and Delworth 2018), the Community EARTH System Model Large Ensemble, CESM-LENS (Kay et al. 2015), and the Max Planck Institute Grand Ensemble, MPI-GE (Maher et al. 2019). While SPEAR LO forced runs follow the CMIP6 SSP5-8.5, the last three ensembles are available with various CMIP5 scenarios. Table 1 provides additional details of the SPEAR and other models' experiments used in this study. Although developed following different protocols (CMIP5 vs CMIP6) and thus being different in the details of the future anthropogenic forcings, SSP5-8.5 and SSP2-4.5 are comparable to RCP8.5 and RCP4.5, respectively, when evaluated on the base of carbon dioxide emissions and radiative forcing Meinshausen et al. 2500).

Model evaluation
An important step in multi-model analysis of regional hydroclimates is to evaluate each model's performance to determine which models more reliably represent the hydroclimate of a certain region. We assessed each model's performance in terms of (1) the annual cycle of precipitation over Central America   Table 2), (2) the amplitude of the interannual, multi-annual, and decadal variability of the Central American summer rainfall (Fig. S5) and ability to simulate the tail of the probability distribution of the 5-year MJJAS rainfall anomalies (that is, multi-annual negative rainfall deficits leading to droughts like the 2015-2019 one, Fig. S6), and (3) the ability to capture the remote SST drivers of Central American rainfall in the Atlantic and Pacific Ocean (Fig. S7).
SPEAR MED has the most realistic representation of the Central American summer rainfall patterns (Fig. S3) and seasonal cycle ( Fig. S4) according to six different precipitation metrics (Table 2), followed by SPEAR LO and FLOR FA. All models reproduce quite well the magnitude of rainfall variability, with the exception of CESM that underestimates it (Fig. S5). The SPEAR models best reproduce the left tail of the 5-year MJJAS rainfall anomaly cumulative probability distribution (Fig. S6). All models capture-though exaggerating it-the remote influence of the Atlantic and Pacific oceans on Central American rainfall (Fig. S7), with MPI-GE and FLOR FA exaggerating the influence of the tropical Pacific and partially missing that of the eastern tropical Atlantic. Collectively, these results suggest that all models used in this study are appropriate tools for use in characterizing future changes in regional precipitation over Central America, with the SPEAR MED and SPEAR LO the best performing models.
Finally, we also evaluate if the models' historical trends of MJJAS rainfall are consistent with observations over Central America over the period 1979-2019 (Fig. S8). To do so, we compute rainfall trends over the period 1979-2019 in GPCC and compared them with individual members of the ALLFORC8.5 ensembles over the same time period. If the observed trend at one grid point was within the range of those simulated by the 30 ensemble members, then we said that the model is consistent with observations in that grid box. We found that models are consistent with observations over most of the grid points (Fig. S8).

Climatological indexes and event definition
As precipitation in Central America is mostly concentrated from May to September, we use a Central American Summer Rainfall Index (SRI hereafter) to monitor the interannual and The Central American Summer Rainfall Index (SRI), the SRI-CLLJ and SRI-TNATNP relationship in the models is compared to observations. RMS, root-mean-square error; MSD, midsummer drought. The CLLJ-SRI correlation is −0.58 in observations (ERA5, GPCC). The TNATNP-SRI correlation is 0.62 in observations (ERSSTv5,GPCC).
STDDEV in the bottom two lines is the ensemble standard deviation monthly rainfall anomalies. The SRI is defined as the area-averaged MJJAS precipitation over all land grid points between 11 • N-18 • N and 95 • W-84 • W (Fig. 1a). This is the area most affected by the 2015-2019 meteorological drought (Fig. 1a), and will therefore be our study region. A more intense Caribbean Low-Level Jet (CLLJ) tends overall to reduce summertime precipitation, especially over the west side of Central America (Wang 2007;Amador 1998;Mo et al. 2005;Cook and Vizy 2010).
To monitor the correlation between the CLLJ and SRI, we introduce, following the definition of (Wang 2007), the CLLJ index as the area-averaged MJJAS 925 hPa zonal wind in the region represented in Fig. 2a (12.5 • N-17.5 • N, 80 • W-70 • W) multiplied by −1. This is justified by the fact that easterly winds feature a maximum (larger than 13 m/s) in the lower troposphere at about 925 hPa (Fig. 2a). The 5-year running mean of the CLLJ index is shown in Fig. 2b. The correlation between CLLJ and SRI in summer is −0.58 in observations, with models generally reproducing such relationship well (Table 2).
Interannual variability of precipitation over Central America is affected by interactions with Pacific and Atlantic Ocean SSTs. Warm Atlantic-cool Pacific (cool Atlantic-warm Pacific) conditions cause a weakening (strengthening) of the CLLJ and favor increased (decreased) precipitation over Central America (Taylor et al. 2002(Taylor et al. , 2011Fuentes-Franco et al. 2015, see also  To quantitatively define the relationship between SRI and the Tropical North Atlantic (TNA)-Tropical North Pacific (TNP) SST difference, we calculate the TNATNP SST index (TNATNP hereafter) as defined by Fuentes-Franco et al. (2015). The TNATNP is the difference between the TNA SST averaged in the region between 5 to 22 • N and from 85 to 35 • W, and the TNP SST averaged in the region between 9 to 27 • N and from 110 to 90 • W (Fig. 2c). As expected, the TNATNP index is positively correlated to SRI (0.6, Table 2), with models generally featuring an even larger correlation (between 0.7 and 0.8). The 5year running mean of the TNATNP anomaly index indicates the persistence of conditions characterized by TNP warmer than the TNA (Fig. 2d).
Hereafter, we refer to the 2015-2019 Central American SRI negative anomaly (Fig. 1d) as "prec event 1519". Averaging all four observation datasets, we obtain a 5-year anomaly of −37 mm/month, which is the largest (in magnitude) measured over the observational record (i.e., 1921-2019). Similarly, for the CLLJ and the TNATNP, we refer to the mean 2015-2019 MJJAS CCLJ and 2015-2019 MJJAS TNATNP anomalies as the "cllj event 1519" and 'tnatnp event 1519". The values of these anomalies are 14.8 m/s and −0.23 K, respectively ( Fig. 2b and c). The CLLJ has been persistently intense during 2014-2019, resulting in record high values of the 5-year running mean values of the CLLJ index (Fig. 2b). Likewise, the mean 2015-2019 SSTs in the tropical Pacific were generally higher than the SSTs in the tropical Atlantic by about 1 K (Fig. 2c), resulting in large negative values of TNATNP (Fig. 2d).

Estimate of probability and risk ratios
The probability of occurrence of prec event 1519 is estimated according to the empirical probability distribution of the 5-year SRI anomalies in the conterfactual climate, i.e., a climate undisturbed by anthropogenic influence (e.g., NAS 2016; Hauser et al. 2017). As noted by Hauser et al. (2017), choices of different counterfactual climates (e.g., pre-industrial vs. early twentieth century) can lead to contradicting results. We account (when possible, see Table 1) for three counterfactual climates: (1) CTRL, (2) NATURAL, and (3) PAST. PAST employs historical simulations forced with observed boundary conditions (i.e., HIST) over time period  when the anthropogenic impact on climate was substantially smaller than today.
As in Pascale et al. (2020), to build the CTRL empirical probability distribution, we randomly selected a 50-year and 5-year sequence (nonoverlapping) and then calculate the anomaly of the 5-year period relative to the 50-year climatology. This choice mimics the 2015-2019 mean minus the 1921-1970 mean. We then repeat this process 10,000 times to form a distribution of the 5-year SRI anomalies. We deal with the NATURAL runs in the same way after concatenating all the ensemble members. For the HIST runs, we concatenate the ensemble members after taking years 1921-1970. To evaluate the decadal change in the probability of occurrence during the historical, present (factual climate) and future projected climate, we empirically estimate a decadalvarying probability distribution of the 5-year MJJAS rainfall anomalies (relative to the 1921-1970 mean) using the ALLFORC experiments (HIST and SSP5-8.5 or SSP2-4.5 in the case of SPEAR MED). Following the method used in Pascale et al. (2020), we estimate the empirical probability distribution for a 20-year time window centered around the referenced year (so, for example, decade 2010 is built from all years from 2001 to 2020). This choice is a trade-off between a time period not too wide (in order to assume the stationarity of the probability distribution) and a number of instances large enough to allow for sufficiently accurate estimates of probabilities of rare events (e.g., 100-y return time). Once we have the decadal probability distribution, we can estimate the probability of occurrence, for each 20-year period, of prec event 1519 for any random 5-year segment within the 20-year time window. We repeat this calculation every 5 years, that is, 1921-1940, 1926-1945, ..., 2076-2095, 2081-2100, to obtain a time-varying estimate of the probability.
We estimate the risk ratio of prec event 1519 by computing the ratio of the probability that the MJJAS precipitation deficit is below the chosen threshold in the factual world and the probability that the MJJAS precipitation deficit is below the same threshold in the counterfactual world. For the probability of the factual world (at the time of the event), we take a weighed mean of the probabilities for the 2015 (i.e., 2006-2025) and 2020 (2011-2030) decades. The 95% confidence intervals in these probabilities were estimated by applying bootstrap-with-replacement resampling. The same approach described so far is applied to the cllj event 1519 and tnatnp event 1519 events.

Recent 40-year trend
The recent negative rainfall trends are common to most of Central America ( Fig. 3a and Fig. S2). Previous studies (Neelin et al. 2006) found that simulations from climate models of the CMIP3 archive (Meehl et al. 2007) generally show a good intermodel agreement on this pattern. This led Neelin et al. (2006) to wonder whether this regional precipitation anomaly is a consequence of ACC. After analyzing CMIP3 (Meehl et al. 2007) models preindustrial runs, they concluded that attribution to ACC of the observed negative trend is plausible but inconclusive.
Here, we first evaluate how likely the 1979-2019 SRI trend is in a counterfactual climate with no anthropogenic influence, i.e., in the long term control runs (CTRL, Table 1). The CTRL 40-year trend distributions (Fig. 3b) suggest that this is indeed a rare yet possible event in four out of five model ensembles analyzed in this study (0.9% in SPEAR MED, 0.6% in SPEAR LO, 0.07% in FLOR FA, 0% in CESM and 1.1% in MPI-GE). Differences in these probability estimates might depend on the model's ability to reproduce the variety of trends due to internal variability. For example, CESM largely underestimates natural variability of Central American summer rainfall compared to observations (standard deviation of 25 mm/month vs 33 mm/month). Nevertheless, we can state that a 40-year trend like that observed recently is likely to be possible even without any anthropogenic influence on the climate.
As a second step of our analysis, we then analyze the 1979-2019 trends from the ALL-FORC ensembles (Table 1), whose forcing includes both natural and anthropogenic effects. Fig. 3c shows the distribution of SRI 1979-2019 trends for each model ensemble. These distributions are centered around their forced value, with their widths determined by internal climate variability. For each model, the trend's ensemble mean, which represents the forced signal, is close to zero (between −0.3 and 0 mm/month/year) and substantially smaller than observations (−1.3 mm/month/year). This indicates a small forced signal relative to the noise generated by internal variability. Large negative trends comparable to the observed trend are achieved only by a few ensemble members (Fig. 3c).
Finally, for models for which we have the NATURAL experiments (SPEAR MED and FLOR FA), we compare the trend distributions across the ensemble members for the ALL-FORC and NATURAL experiments (Fig. 3d). This can be very informative if we recall that ALLFORC and NATURAL differ only by the presence of anthropogenic forcing. Thus, any  Table 1). The observed SRI negative trend, with its 95% significance confidence interval, is shown by the black vertical line and grey bar. c Ensemble distribution of the 1979-2019 SRI trends. The trend ensemble mean, for each ensemble, is shown by the color marks on the upper x-axis. The pale blue denotes the range of the forced (i.e., ensemble mean) trends. d Ensemble distribution of the 1979-2019 SRI trends from the ALLFORC and NATURAL (Table 1) ensembles. Precipitation trends are estimated using least-squares linear trends for each grid point. Significance of trends at the 5% level is assessed through a t-test significant difference in the two distributions is to be attributed to anthropogenic forcing. We find that, for both SPEAR MED and FLOR FA, the two distributions are qualitatively very similar (Fig. 3d). A Kolmogorov-Smirnov test performed on the two arrays of trends, for each model, reveals that the probability that the two arrays are drawn from the same distribution is 94% for both models. This test so does not pass the 5% level, although it falls close to it (i.e., 95%). This results so suggests, at least for these two models for which we have NATURAL, that the addition of anthropogenic forcing to natural forcing does not lead to a substantial difference in the 40-year trends distribution across the ensemble.
The overall evidence presented so far suggests that the observed recent drying trend may be mostly due to internal, multi-decadal variability rather than ACC. However, we admit that we cannot totally exclude the role of ACC. This is because we are unable to make a more precise statement (of the form "ACC has increased the likelihood of the 1979-2019 SRI trend by X times") given the scarcity of points, thirty, used to reconstruct the probability distribution of the 1979-2019 SRI trends (Fig. 3d). Instead, we would have needed hundreds to thousands of equivalent (i.e., under the same 1979-2019 natural only and natural plus anthropogenic forcing) estimates of the 1979-2019 SRI trends. These would have then allowed us to reconstruct reliably the probability distribution of the 1979-2019 SRI trends and so estimate the probability of having a 40-year trend smaller (i.e., more negative) than the one observed in 1979-2019 in both the factual and counterfactual climate. Figure 4 shows, for each model ensemble, the risk ratio of prec event 1519. These values are obtained by the ratio of the probabilities of occurrence of prec event 1519 in the factual world and the probability of occurrence of prec event 1519 in three different counterfactual climates (CTRL, NATURAL and PAST, see Section 2.5). For SPEAR MED, we have two different sets of values, depending on whether years 2015-2030 are drawn from SSP2-4.5 or SSP5-8.5 (Fig. 5a).

Attribution of the 2015-2019 drought event
All ensembles, with the exception of SPEAR MED ALLFORC4.5, indicate a risk ratio significantly larger than one for all three counterfactual climates. The difference between SPEAR MED ALLFORC4.5 risk ratio and SPEAR MED ALLFORC8.5 risk ratio originates from years 2015-2030, which lead to significantly different values of the probability of occurrence of prec event 1519 (Fig. 5a). These differences may be in part explained by the differences in radiative forcings between SSP2-4.5 and SSP5-8.5, which are minimal until 2020-2030 and start differing substantially only afterwards ). Differences may be also attributable to the limited number of ensemble members (30), which does not allow a complete quantification of internal variability (Fig. 5a).
For the second part of the twenty-first century, the probabilities of occurrence of prec event 1519 in SSP2-4.5 and SSP5-8.5 diverge substantially (Fig. 5b). In SSP2-4.5 (middle of the road pathway) it remains below 20% by the end of this century. Given a probability of around 2% of prec event 1519 in the counterfactual climates, this implies a risk ratio which remains below 10. In SSP5-8.5 (high emission scenario), this probabilities  increases swiftly, reaching 90% by the end of the century, that is a risk ratio of nearly 40. SPEAR LO and FLOR FA share similar values (median values between 2 and 5). For MPI (100 ensemble members), we use RCP8.5 for years 2005-2030, but we find no significant differences when RCP4.5 and RPC2.6 are instead used (not shown).
Combining all different risk ratios for different models, using a simple average to synthesize the results, the probability of an event like the observed 2015-2019 meteorological drought has increased by a factor of 4 relative to the two counterfactual climates PAST and CTRL, with confidence intervals (3, 7) and (3.5, 5) respectively, and by a factor 1.5 (0.95, 2) relative to NATURAL. Let us note, however, that the value for NATURAL is obtained averaging only three cases, so we should bear in mind that is a less comprehensive and somewhat biased result. Overall, results summarized in Fig. 4 suggest that ACC has significantly increased the likelihood of a multi-year drought like that occurred in Central America between 2015 and 2019.

The role of the CLLJ and Tropical Atlantic-Tropical Pacific temperature SST differences.
The 2015-2019 meteorological drought unfolded in association with a very strong CLLJ (Fig. 2a) and a negative TNATNP (Fig. 2c). In particular, the 5-year mean of the CLLJ index features the largest positive anomalies over the entire ERA5 time span (Fig. 2b) while the 5-year mean of TNATNP features a negative value equalled only other two times in the last hundred years (Fig. 2d).
The observed mean 2015-2019 925 hPa zonal wind anomalies and SST anomalies have a strong resemblance to those projected for the end of this century under the highest emission scenario (Fig. 6). Strengthening of the CLLJ and faster warming up of the tropical Pacific than the tropical north Atlantic is indeed a robust feature across climate models (Rauscher et al. 2008Fuentes-Franco et al. 2015). It is therefore natural to ask whether the likelihood of the unusual 2015-2019 CLLJ and TNATNP anomaly was increased by ACC. To address this question, we estimate the probability of occurrence of cllj event 1519 and tnatnp event 1519 (Section 2.4) in the factual and counterfactual climates. We perform this  -8.5. b Projected (2071-8.5. b Projected ( -2100-8.5. b Projected ( vs. 1921-8.5. b Projected ( -1970 MJJAS SST anomalies (black contours show the climatological mean) under  analysis only for SPEAR MED and SPEAR LO, for which we have data readily available. However, given the similar future projections in terms of CLLJ and SSTs, we speculate that similar results would be obtained from other ensembles too. The CLLJ index features a long term significant trend in both NATURAL and HIST (SPEAR MED). Since this trend is not seen in reanalyses, and there is no long-term forcing in NATURAL which can physically justify it, it must be a model drift unrelated to anthropogenic forcing. We therefore subtract from the CLLJ index the trend estimated from NATURAL to both NATURAL and HIST before estimating decadal probabilities of occurrence and risk ratios. Figure 7 shows the risk ratio of cllj event 1519. We have contrasting results across different counterfactual climates (e.g., in SPEAR LO) and for the same model (i.e.,  . Overall, there is not strong evidence for the 2015-2019 CLLJ mean anomaly to be attributable to ACC. In SPEAR MED (SSP5-8.5), risk ratios of cllj event 1519 are close to one (indicating no attribution to ACC) even though those of prec event 1519 are clearly larger than one (Fig. 4). Similarly, in SPEAR LO, risk ratios for prec event 1519 are clearly larger than one, but risk ratios for cllj event 1519 are larger than one only when PAST is chosen as the counterfactual climate. This means that for almost all cases in which prec event 1519 is attributable to ACC, the 2015-2019 CLLJ anomaly is not attributable to ACC. While the CLLJ has not shown fluctuations unambiguously attributable to ACC during 2015-2019, the decadal evolution of the probability of cllj event 1519 in SPEAR MED suggests that this will be the case in the next few decades ( Fig. 9a and b). Figure 8 shows the risk ratios associated with tnatnp event 1519. Results here are much more consistent across the SPEAR models and with those about the prec event 1519 event, and they indicate that ACC has increased the likelihood of the tnatnp event 1519 by a factor 2 (with a 95% confidence interval between 1.5 and 5). The NATURAL data are lacking for this case, although those do not generally lead to very differ results from those of, e.g., CTRL (Fig. 7) and therefore that does not alter the validity of our conclusions.
The emergence of the ACC impact on TNATNP is faster than that on CLJJ, as evident from the time evolution of the probability of occurrence (e.g., Fig. 9a vs. c, b vs. d). We conjecture that may be due to the larger level of background noise of a dynamical, atmospheric field as the CLLJ as compared the more thermodynamically constrained SST field. This leads us to speculate that-differently from SST, for which ACC is already detectable (Chan and Wu 2015)-the signal associated with ACC on the regional circulation, specifically the CLLJ, has not yet emerged from the noise of internal variability. For example, previous work on the strengthening of the North Atlantic Subtropical High-which is connected to the CLLJ, although it is not the only driver of the CLLJ-is also inconclusive about the time of emergence of the ACC forced signal (Li et al. 2011, 2013, Diem 2013.

Discussion and conclusions
In this study, we have analyzed the recent drying trend observed from the late 1970s and the 2015-2019 Central American meteorological drought to determine whether these two events have been made more likely by ACC.
From the comparison of large ensembles driven by natural forcing only and natural plus anthropogenic forcing, we conclude that the observed 1979-2019 drying trend is consistent with internal, multi-decadal variability, although our analysis on this was not conclusive and therefore we cannot totally exclude the role of ACC. Instead, we found strong evidence that the 2015-2019 rainfall deficit was made four times more likely by ACC (Fig. 4). Closely related to that, we also found that ACC increased the probability (1.5 to 3), for the period 2015-2019, of large, negative SST differences between the Tropical North Atlantic and Tropical North Pacific, which are known to drive meteorological droughts in Central America.
At this point it is worth speculating on why we managed to attribute the 2015-2019 rainfall deficit to ACC while we did not for the 1979-2019 trend. First, as already discussed in Secttion 3.1, we are unable to quantify a risk ratio for the 1979-2019 trend event-and thus quantify by how many times ACC made it more or less likely-because we cannot empirically reconstruct the probability distribution for this event due to scarcity of data. Second, we note that the effect of ACC starts becoming noticeable on the 5year SRI distribution's tail only after 2010 in SPEAR MED (Fig. 5a) and after 2000 in FLOR FA. Therefore, while this allows us to see the effects of this climatic change on a short-term event like prec event 1519, it may not directly lead to a detectable effect on the 1979-2019 trend, whose estimate is based, for the most part, on years during which the effect of ACC was not yet detectable (at least according to these two models).
The extent of the precipitation anomaly observed between 2015 and 2019 revealed a key role for rainfall in this multi-year drought, leading us to focus on precipitation as the dominant variable to characterize this event. Indeed, recent studies highlight the projected precipitation reduction over Central America as the key factor in the increase of drought risk in the coming decades (Giorgi 2006;Rauscher et al. 2008;Fuentes-Franco et al. 2015;Almazroui et al. 2021;Depsky and Pons 2020). However, it is also important to mention that regional ongoing warming trends (Aguilar et al. 2005;Gourdji et al. 2015) play an important role in driving enhanced evapotranspiration and thus further exacerbating future droughts Alfaro-Cordoba et al. 2020). An example of that is the 2013-2016 Pan-Caribbean drought, which was not driven by a precipitation deficit but it was a direct consequence of exceptionally high evapotranspiration due to higher temperatures (Herrera et al. 2018). Furthermore, higher surface temperature from surrounding oceans along with drier land conditions can induce stronger lower-to-middle tropospheric stability, especially in the early summer, by increasing the midtropospheric dry static energy relative to near-surface moist static energy (Seth et al. 2011(Seth et al. , 2013Giannini 2010). Increased lower-to-middle tropospheric stability suppresses convection, reducing precipitation in monsoonal regimes where summertime precipitation is mostly associated with convective systems (Pascale et al. 2017(Pascale et al. , 2019. Higher temperatures therefore may exacerbate the consequences of a meteorological drought in multiple ways. Future work should include analysis of the risk associated with compound events in Central America driven by exceptionally high temperatures and prolonged rainfall deficits. Although a multi-year meteorological drought like the one that hit Central America between 2015 and 2019 is a rare event, anthropogenic climate change has significantly increased the likelihood for such a prolonged rainfall deficit to occur. All of the models analyzed in this study suggest that the probability of occurrence of this event will continue to increase into the future, with a rate depending on the specific socioeconomic scenario (Eyring et al. 2016;O'Neill et al. 2017;Riahi et al. 2017). As a consequence, we expect that prolonged multi-year periods of below-average precipitation in the region will be more severe in the coming decades.
The more frequent occurrence of multi-year meteorological droughts in Central America shown in this study, along with increased aridity due to higher temperatures and the likelihood of increased water demand due to a growing population (Hidalgo 2021), suggests that managing water resources will likely be increasingly challenging in the decades to come. Advanced planning for drought years and measures to increase water efficiency for agricultural and urban use is therefore urgently needed at the national and regional level to better adapt to the new hydroclimatic conditions.
Funding Open access funding provided by Alma Mater Studiorum -Università di Bologna within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.