Streamflow Intensification Driven by the Atlantic Multidecadal Oscillation (AMO) in the Atrato River Basin, Northwestern Colombia

The impact of the Atlantic Multidecadal Oscillation (AMO) on the variations in the streamflow in the Atrato River Basin (ARB) during the 1965–2016 period was analyzed here by considering the cold (1965–1994) and warm (1995–2015) phases of this oscillation. The mean streamflow increased after 1994 (AMO phase change). This increase is related to the strengthening of the zonal gradients of the sea surface temperature (SST) and sea level pressure (SLP) between the tropical central Pacific and the tropical Atlantic after 1994 (warm AMO phase). These gradients contributed to strengthen the Walker cell related upward movement over northern and northwestern South America, in particular during November-December (ND). Consistently, the frequency (R20 mm) and intensity (SDII) of extreme daily rainfall events increased during the 1995–2015 period. Our results show a connection between the AMO and the increase in the streamflow in the ARB during the last five decades. These results contribute to the studies of resilience and climate adaptation in the region.

The El Niño Southern Oscillation (ENSO) in the interannual time scales is the main ocean-atmospheric coupled phenomenon that affects the rainfall distribution in tropical regions and causes droughts and/or floods in the watersheds involved [15,18,19]. Several studies have documented the ENSO effects on rainfall and streamflow in northern and northwestern South America (SA) [1,18,[20][21][22][23][24][25] with consequent negative impacts on the economy, biodiversity and human health in Colombia [6,25,26]. The ENSO warm phase or El Niño (EN) is associated with negative anomalies of the rainfall, soil moisture, evapotranspiration, and streamflow, mainly in central, northern, and western Colombia [1,20,25,27,28], which are caused by variations in both the Walker and Hadley cells [29][30][31]. The La Niña (LN) events are related to the opposite anomaly patterns of these variables. At interdecadal scales, the average discharges of Colombian rivers and precipitation correlate significantly and negatively with the North Atlantic Oscillation (NAO) index and the Pacific Decadal Oscillation (PDO) index and the rainfall in the North Caribbean stations correlate significantly and positively with the sea surface temperature (SST) in the tropical Atlantic [13,15,32].
However, the ARB streamflow might present low-frequency variability which is modulated by the SST low-frequency variability mode. One of these modes in known as the Atlantic Multidecadal Oscillation (AMO) in the North Atlantic [33,34]. This mode is characterized by anomalies of the same signal in the North Atlantic with two centers of maximum anomalies at 15 • N and 55 • N and a periodicity between 65 and 80 years. The warm phase (positive) occurred during the 1860-1885 and 1930-1964 periods and the cold phase (negative) during the 1895-1924 and 1965-1994 periods [35].
Several studies have provided evidence of the AMO influence on the interannual climate variability in many regions in SA [36][37][38][39][40][41]. Chiessi et al. [40] indicated the significant effect of the AMO on the SA monsoon system, such that during the warm AMO (WAMO) phase the main area with summer rainfall is displaced northward, affecting the rainfall and runoff in the region. Dong et al. [41] showed that the amplitude of the Niño 3 index during WAMO is reduced in relation to the CAMO, so that the AMO can modulate ENSO variability. Additionally, Kayano and Capistrano [36] demonstrated that EN events during the cold AMO (CAMO) phase were generally stronger than those during WAMO phase and that the EN-related rainfall anomalies were more organized and extensive during the CAMO than in the other AMO phase. They also showed that LN events were stronger during the WAMO than during the CAMO. Timmermann et al. [39] found that the strong (weak) ENSO variability during the CAMO (WAMO) phase reduces (increases) the magnitude of the SST annual cycle in the tropical Eastern Pacific, and modulates the ENSO effect on the South American rainfall. In addition, Kayano et al. [38] indicated that the negative multidecadal relation between AMO and the PDO, during the WAMO reinforces the anomalous cooling of the tropical Pacific and favors the settlement of LN events. This leads to a configuration of the dipolar pattern with positive rainfall anomalies in northern SA and negative anomalies in central and eastern SA, especially from winter to summer, which is modulated by the anomalous Walker and Hadley cells. In turn, during the CAMO and warm PDO phase, the warming in the tropical Pacific is reinforced and the LN pattern and its effects on SA rainfall are weakened.
Floods, droughts, and other extreme hydrometeorological events have become recurrent natural disasters over the past several decades, and are expected that these types of hydrometeorological threats will intensify in the future [42,43]. The 2010-2012 LN event caused extreme rainfall and flooding in Colombia with disasters in the Magdalena and Cauca river valleys [44,45] and in the Atrato River valley, a region associated with a high degree of social vulnerability in the country [46]. The ARB is located in extreme northwestern Colombia known as the Chocó Biogeographic ( Figure 1) and is one of the richest regions in the world in terms of natural and cultural diversity [47][48][49]. It is a region with some of the highest rainfall levels on earth, with an annual mean rainfall greater than 12,500 millimeters [7,50,51]. Therefore, understanding this region's extreme hydrometeorological events is a priority for the conservation and planning of water resources of a region that shows strong contrasts.
Water 2020, 12, x FOR PEER REVIEW 3 of 23 and is one of the richest regions in the world in terms of natural and cultural diversity [47][48][49]. It is a region with some of the highest rainfall levels on earth, with an annual mean rainfall greater than 12,500 millimeters [7,50,51]. Therefore, understanding this region's extreme hydrometeorological events is a priority for the conservation and planning of water resources of a region that shows strong contrasts. In this sense, a better understanding of the influence of the low-frequency changes associated with the AMO in the ARB streamflow is aimed in this study, with the main focus on the peaks during May-June (MJ) and November-December (ND). The years with extreme ARB streamflow were stratified according to the AMO phases during the last five decades . For these years, the related anomalous oceanic and atmospheric patterns and the anomalies in the total rainfall and extreme rainfall indices were examined. Section 2 describes the data and methodology. The mean anomaly fields for SSTs, sea level pressures (SLP), vertical velocities at 500 hPa and rainfall are presented and discussed in Section 3. Finally, the conclusions are presented in Section 4. In this sense, a better understanding of the influence of the low-frequency changes associated with the AMO in the ARB streamflow is aimed in this study, with the main focus on the peaks during May-June (MJ) and November-December (ND). The years with extreme ARB streamflow were stratified according to the AMO phases during the last five decades . For these years, the related anomalous oceanic and atmospheric patterns and the anomalies in the total rainfall and extreme rainfall indices were examined. Section 2 describes the data and methodology. The mean anomaly fields for SSTs, sea level pressures (SLP), vertical velocities at 500 hPa and rainfall are presented and discussed in Section 3. Finally, the conclusions are presented in Section 4.

Study Area
The ARB is located between 5 • and 9 • N and between 75 • and 80 • W covering an area of approximately 37,730 km 2 and an altitude range of 1-3820 meters above the sea level (Figure 1a). The Atrato river is born in the Western Cordillera, running 600 km before delivering nearly 4140 m 3 ·s −1 each year to the Urabá Gulf in the Caribbean Sea [52][53][54][55]; one of the longest navigable rivers in Colombia, only exceeded by the Magdalena and the Cauca Rivers [56]. Moreover, its specific streamflow of 0.109 m 3 ·s −1 km −2 , is among the world's highest runoff rates [7,25,53,54].
The precipitation annual cycle in the ARB (Figure 1b) is dominated by a uniform regime with rainfall throughout the year, especially between April and December, that exceeds the 450 mm·month −1 [57,58]. According to Velasquez-Restrepo et al. [53] and the Institute of Hydrology, Meteorology and Environmental Studies (IDEAM) [54], the upper basin (i.e., upstream Bellavista stream gauge station) has a drainage area of 15,358 km 2 , a mean annual streamflow of 2500 m 3 ·s −1 for the 1965-2015 period, two notable peaks during May-June (MJ, 2640 m 3 ·s −1 ) and November-December (ND, 2820 m 3 ·s −1 ), and the lowest streamflow (1770 m 3 ·s −1 ) in March ( Figure 1c). So, considering the streamflow annual cycle the analyses focus on peaks during MJ and ND.

Hydrometeorological Data
The stream flow data of the 1965-2015 period at the stream gauge stations of Belén, San Antonio, and Bellavista (at 26, 18 and 12 m a.s.l., respectively) from the IDEAM, located in upper ARB were used in preliminary analyses (Appendix A) based on the principal component analysis (PCA) and correlations among the stations. These analyses indicated that the Bellavista streamflow time series reproduced most of the variations shown in the other two streamflow time series. So, the streamflow data at Bellavista was used here as representative of the study domain.
Monthly precipitation data were obtained from the observed meteorological network (from IDEAM). Only the most extensive records at Tagachi, Buchado, and Bellavista were considered ( Figure 1a). Hence, the average precipitation of these three stations represent the average precipitation in the ARB. Furthermore, daily precipitation data from the Global Meteorological Forcing Dataset for Land Surface Modeling (GMFD) [59] at a horizontal resolution of 0.25 • × 0.25 • were also used. The GMFD data result from a combination of the global datasets and reanalysis from NCEP-NCAR reanalysis, Climatic Research Unit monthly climate variables, Global Precipitation Climatology Project (GPCP) daily precipitation, Tropical Rainfall Measuring Mission three-hourly precipitation and the NASA Langley monthly surface radiation budget. This dataset provides a significant improvement over the original reanalysis variables and has been used for a wide variety of applications and diagnostic studies in hydroclimatology [24,60,61]. This data was obtained in the area limited at 65 • W, 85 • W, 20 • N and 5 • S for the 1965-2015 period.

Atmospheric and Oceanic Data
Monthly reconstructed SST dataset (National Oceanic and Atmospheric Administration/Extended Reconstructed Sea Surface Temperature-NOAA/ERSST) [62], version 5, derived from the International Comprehensive Ocean-Atmosphere Dataset [63] was used. The monthly SLP and 500 hPa vertical velocity in pressure coordinate data were obtained from the National Center for Environmental Prediction and National Center for Atmospheric Research (NCEP/NCAR) [64] Reanalysis 1 project [65]. The SST data have a latitude-longitude resolution of 2 • × 2 • and the SLP and vertical velocity data, of 2.5 • × 2.5 • . These atmospheric and oceanic data were obtained in the domain limited at 120 • E, Greenwich longitude, 40 • N and 40 • S during the 1965-2015 period. The monthly standardized anomalies of these variables were based on the 1965-2015 period. Additionally, the AMO index used in this study was accessed at (https://www.esrl.noaa.gov/psd/data/timeseries/AMO/) [35]. This index is based on the Kaplan SST data.

Extreme Rainfall Index
Daily rainfall data from GMFD were used to calculate the extreme rainfall indices proposed by the Expert Team on Climate Change Detection and Indices (ETCCDI) [66,67]. According to Folland et al. [68], changes in these extremes often have greater impacts than the changes in the mean values and can be climate risk indicators. Two indices related to extreme rainfall were calculated at the monthly scale. The frequency index (R20 mm) characterizes heavy rains and is the number of days in a month for which the daily rainfall is greater than or equal to 20 mm [69]. The simple daily intensity index (SDII) is the ratio between the total annual rainfall and the number of wet days (≥1 mm). These indices have been widely used to quantify severe rainfall and to detect possible relationships between hydrological threats and areas that are potentially vulnerable to catastrophic events [24,[70][71][72][73][74][75].

Mean SST Anomaly Patterns During the AMO Phases
Following the methodology proposed by Kayano et al. [76,77], the reconstructed SST data of the 1870-2015 period were used only to support that the SST anomaly patterns during the 1965-1994 and 1995-2015 periods indeed contain, respectively, the CAMO and WAMO features [35,78]. In order to remove the climate change effects, the linear trend for the 1870-2015 of the time series at each grid point was eliminated using the method of least squares. Monthly standardized anomalies based on the 1965-2015 period were obtained in the domain limited at 60 • N, 60 • S, 120 • E and the Greenwich longitude. Then, the bimonthly SST anomaly mean states were obtained for MJ and ND during the 1965-1994 and 1995-2015 periods. The statistical significance was assessed using the Student's t-test with a confidence level of 95% and the degrees of freedom are the number of years in each period.

Composite Analysis
The extreme streamflow events in MJ and ND of the 1965-2015 period were selected and stratified according to AMO phases. Bimonthly streamflow time series for MJ and ND at Bellavista stream gauge station were obtained. These bimonthly series were then divided into two time series, one covering 30 years of the 1965-1994 period and the other, 21 years of the 1995-2015 period. These periods are referred to as the periods before and after 1994. Subsequently, the percentile time-series (R) were constructed by ranging the years from one to 30 (one to 21) and dividing them by 30 (21). Then, for each period, the extreme events above the percentile 75 were selected.
First, monthly standardized anomalies of precipitation, SST, SLP, and 500 hPa vertical velocity in pressure coordinate during the 1965-2015 period were obtained and used in the composite analyses. In this case, the means and standard deviations were based on the 1965-2015 period. Bimonthly composites of the anomalies of the rainfall, SST, SLP, and 500 hPa vertical velocity in pressure coordinate were obtained for MJ and ND. Furthermore, the composites for the extreme rainfall indices R20 mm and SDII were also obtained for MJ and ND. Composite analyses of the extreme events before and after 1994 were performed separately. The difference composites between these two periods are also presented.
The statistical significance of the composites was evaluated using the Student's t-test at 95% confidence level, and the degrees of freedom, the number of events for each period. The significance test used for comparing two means [79] considers that a variable X with n values, standard deviation S and Student's t distribution has a significant mean when the absolute value of the mean exceeds t α(n−1) S/ (n − 1); in which t α(n−1) corresponds to the cumulative probability for (n − 1) degrees of freedom in Student's table and α = 0.05. The two-sample t-test of Levine and Wilks was used to evaluate the statistical significance of differences in the composites before and after 1994 [80]. This parametric test considers two variables, X 1 and X 2 , with values n 1 and n 2 , standard deviations S 1 and S 2 and composites indicated by X 1 and X 2 and assumes that the difference X 1 − X 2 has a Student t distribution with n 1 + n 2 − 2 degrees of freedom and that the absolute values of X 1 − X 2 greater than t α(n 1 +n 2 −2) (n 1 − 1)S 2 1 + (n 2 − 1)S 2 2 n 1 +n 2 n 1 n 2 (n 1 +n 2 −2) are statistically significant. Figure 2 shows the bimonthly SST anomaly patterns during the two selected periods. During 1965-1994, both seasons show significant negative anomalies in most of the North Atlantic and in the tropical West Pacific, and the negative ones in the extratropical South Atlantic and Southeast Pacific, with the negative anomalies in Tropical North Atlantic (TNA) being significant during ND and non-significant during MJ (Figure 2a,c). The maps during the 1995-2015 period show inverted signal patterns (Figure 2b,d). These maps show previously documented typical AMO features [35,39,76,81]. Thus, the two selected periods present contrasting SST mean states, under which the streamflow and rainfall in the ARB may have distinct characteristics which may be related to the low-frequency background associated with the AMO.

SST Mean States
Water 2020, 12, x FOR PEER REVIEW 6 of 23 Figure 2 shows the bimonthly SST anomaly patterns during the two selected periods. During 1965-1994, both seasons show significant negative anomalies in most of the North Atlantic and in the tropical West Pacific, and the negative ones in the extratropical South Atlantic and Southeast Pacific, with the negative anomalies in Tropical North Atlantic (TNA) being significant during ND and nonsignificant during MJ (Figure 2a,c). The maps during the 1995-2015 period show inverted signal patterns (Figure 2b,d). These maps show previously documented typical AMO features [35,39,76,81]. Thus, the two selected periods present contrasting SST mean states, under which the streamflow and rainfall in the ARB may have distinct characteristics which may be related to the low-frequency background associated with the AMO.

Difference of the Average Precipitation in the ARB and Streamflow in Bellavista
The differences between the two AMO phases of the bimonthly streamflow at Bellavista and precipitation in the ARB in MJ and ND were evaluated using the Student's t-test (Table 1)). The t-test has been used to compare the differences between climate and environmental time series [82,83]. The normality of the series was evaluated using the Shapiro and Wilk [84] and Anderson and Darling [85] tests. All time series satisfied the normality test. The differences in the means were analyzed at a significance level of 5%. Table 1 shows the t-test comparisons of streamflow and precipitation time

Difference of the Average Precipitation in the ARB and Streamflow in Bellavista
The differences between the two AMO phases of the bimonthly streamflow at Bellavista and precipitation in the ARB in MJ and ND were evaluated using the Student's t-test (Table 1)). The t-test has been used to compare the differences between climate and environmental time series [82,83]. The normality of the series was evaluated using the Shapiro and Wilk [84] and Anderson and Darling [85] tests. All time series satisfied the normality test. The differences in the means were analyzed at a significance level of 5%. Table 1 shows the t-test comparisons of streamflow and precipitation time series before and after 1994 in MJ and ND. In both seasons, the mean streamflow at Bellavista shows a significant difference before and after 1994 and a statistically significant increase (p < 0.05) in the mean discharge from the period before to after 1994. Also, there was a statistically significant increase in the average precipitation in the ARB from the period before to after 1994 in ND.  Figure 3 shows the ARB streamflow and AMO index, both filtered with a 10-year moving mean filter. In both seasons, the streamflow presents a significant (p < 0.05) increase in the mean values of the period before and after 1994, which accompanies the change of the AMO phase. The mean streamflow values before 1994 and after were of 2530 m 3 ·s −1 and 2800 m 3 ·s −1 in MJ, and 2615 m 3 ·s −1 and 3106 m 3 ·s −1 in ND (Figure 3). These results are consistent with the significant correlations of 0.68 (p < 0.05) between the filtered MJ streamflow and AMO index time series, and 0.75 between the filtered ND streamflow and AMO index time series. Similar graphics were constructed for the average precipitation in the ARB (Figure 4). The ND precipitation time series shows a significant (p < 0.05) increase with the means values before and after 1994 of 472 mm and 579 mm. The filtered streamflow and precipitation time series show higher correlation in ND (r = 0.50) than in MJ (r = 0.43). Therefore, there were estimated increases from the period before to after 1994 of 12% and 19% in the streamflow in MJ and ND, and of 21% in the precipitation in ND.

Composite Analysis
Composites considering the streamflow extreme events were calculated in the periods before and after 1994. The extreme events selected for the CAMO and WAMO phases are listed in Table 2.

Composite Analysis
Composites considering the streamflow extreme events were calculated in the periods before and after 1994. The extreme events selected for the CAMO and WAMO phases are listed in Table 2.

Composite Analysis
Composites considering the streamflow extreme events were calculated in the periods before and after 1994. The extreme events selected for the CAMO and WAMO phases are listed in Table 2. Figures 5 and 6 show the SST, SLP, and 500 hPa vertical velocity anomaly patterns during the periods before and after 1994 and their differences. The SST composite during the CAMO shows significant negative anomalies in the subtropical North Atlantic, the Tropical South Atlantic (TSA), the equatorial Western and subtropical South Pacific during MJ and in TNA, and most of the tropical Pacific during ND (Figure 5a,g). Additionally, small areas of anomalous warming occurred in the tropical central Pacific in MJ and in the subtropical areas of the Pacific and TSA in ND. During the WAMO, the SST anomaly patterns were characterized by an anomalous warming in a boomerang-shaped area in the Australasian region in both bimonthly periods and in most of the subtropical North Atlantic in MJ and in the TNA in ND, while negative anomalies were depicted in the central tropical Pacific in MJ (not significant) and ND (significant) (Figure 5b,h). This analysis therefore shows that the greatest differences in SST composites between the two AMO phases occurred in the TNA, subtropical North Atlantic, and central tropical Pacific, which defined the SST variations associated with the low-frequency background shown in Section 3.1 (Figure 5c,i). Furthermore, the SST anomalies exhibited higher magnitudes and greater spatial extents in the tropical sector in ND than in MJ.   (Figure 5b,h), which creates a positive feedback modulating the climate variability in these Oceanic basins. Previous studies argued that the WAMO/CPDO low-frequency mean state strengthens the cooling in the tropical Pacific and, therefore, leads to a higher percentage of La Niña (LN) events [38,86]. This process explains why the LN-related SST anomaly patterns in the tropical Pacific are stronger after 1994 (Figure 5h). These results reinforce the fundamental role of the Atlantic in determining the intensities of the LN events and, consequently, the rainfall in SA [36,38,[87][88][89]. Coherently, all events selected in ND after 1994 coincided with the LN event occurrences ( Table 2). Choi et al. [90] showed that the ENSO decadal oscillation presents a stronger relationship with the SST mean state during its mature stage (November-January). Li et al. [91] and Barichivich et al. [92] observed, since the beginning of the 1990s, a strong warming trend in the tropical Atlantic and Indo-Western Pacific Oceans and a cooling trend in the Eastern Pacific, which intensifies the LN response in the tropical Pacific. The results in the present study confirm their findings (Figure 5c,i). Consistent with the SST composites, the SLP anomaly patterns after 1994 show an abnormally high pressure center in the central and eastern tropical Pacific, and an anomalous low pressure center in the TNA, both centers more intense in ND than in MJ (Figure 5e,k). These patterns lead to a stronger inter-basin SLP gradient, especially during ND (Figure 5k).   Consistent with the SST and SLP anomalies, the anomalous Walker cells before 1994 featured downward motion in the tropical Eastern Pacific and northern and northwestern SA and upward motion in the central tropical Pacific in MJ, and downward motion over most of the tropical Pacific between the South American coast and 180° W and upward motion on the northern coast of SA and northeastern Brazil in ND (Figure 6a,d). On the other hand, after 1994, the SST and SLP gradients between the tropical central Pacific and the tropical Atlantic strengthened the Walker cell with an upward motion to northern and northwestern SA and downward in the tropical central Pacific and the Caribbean Sea in both bimonthly periods (Figure 6b,e). The composite differences represent the intensification (weakening) of the SLP in the Pacific (tropical Atlantic and northwestern SA) associated with anomalous cooling (warming) in the tropical Pacific (TNA and subtropical North Atlantic) and the intensification of Walker circulation with downward (upward) motion in the tropical central Pacific (northern and northwestern SA) after 1994 (Figure 6c,f). Since the Walker circulation connects the tropical oceans, and its upward branch extends to northern and northwestern SA, intensification of the atmospheric circulation and the convergence of humidity on these regions are expected, increasing the rainfall and the streamflow in the ARB, with the highest intensity in ND, when the pattern is better established. The results show a physical connection between the ARB streamflow and the recent SST and SLP inter-basin gradients, boosted by the Atlantic warming during the WAMO, as observed in other studies [91,[93][94][95]. Consistent with the SST and SLP anomalies, the anomalous Walker cells before 1994 featured downward motion in the tropical Eastern Pacific and northern and northwestern SA and upward motion in the central tropical Pacific in MJ, and downward motion over most of the tropical Pacific between the South American coast and 180 • W and upward motion on the northern coast of SA and northeastern Brazil in ND (Figure 6a,d). On the other hand, after 1994, the SST and SLP gradients between the tropical central Pacific and the tropical Atlantic strengthened the Walker cell with an upward motion to northern and northwestern SA and downward in the tropical central Pacific and the Caribbean Sea in both bimonthly periods (Figure 6b,e). The composite differences represent the intensification (weakening) of the SLP in the Pacific (tropical Atlantic and northwestern SA) associated with anomalous cooling (warming) in the tropical Pacific (TNA and subtropical North Atlantic) and the intensification of Walker circulation with downward (upward) motion in the tropical central Pacific (northern and northwestern SA) after 1994 (Figure 6c,f). Since the Walker circulation connects the tropical oceans, and its upward branch extends to northern and northwestern SA, intensification of the atmospheric circulation and the convergence of humidity on these regions are expected, increasing the rainfall and the streamflow in the ARB, with the highest intensity in ND, when the pattern is better established. The results show a physical connection between the ARB streamflow and the recent SST and SLP inter-basin gradients, boosted by the Atlantic warming during the WAMO, as observed in other studies [91,[93][94][95].
Our results are in agreement with the finding that the atmospheric response to the tropical Atlantic warming can lead to an intensification of the upward branch of the Walker circulation in northern SA [96][97][98]. This process is accompanied by a downward motion on the central equatorial Pacific and upward in the Indo-Pacific region. These anomalous motions increase the zonal pressure gradient between the tropical Pacific and Atlantic and increase the divergence at the surface, generating a cooling in the central equatorial Pacific, thus favoring the development of LN events. This is in agreement with the above analyzed SST anomaly patterns with large negative SST anomalies in the tropical eastern Pacific after 1994 in ND (Figure 5h).
Consistent with the anomalous SST and SLP patterns with cooling in the tropical Pacific, especially during WAMO, positive rainfall anomalies were observed in Colombia in both bimonthly periods (left column in Figures 7 and 8). After 1994, the positive rainfall anomalies increased mainly in the Caribbean region and in the Andes region adjacent to the ARB in MJ (left column in Figure 7). After 1994, significant positive anomalies with larger magnitudes in ND than in MJ covered a large part of the Colombian territory, including the ARB (left column in Figure 8). The rainfall increases are associated with the strengthened zonal SST and SLP gradients between the tropical central Pacific and Atlantic in ND (Figure 5h,k). These gradients during ND strengthen the anomalous upward motion in northern and northwestern SA, increasing the rainfall in the ARB and the Colombian Caribbean. An examination of the anomalous rainfall characteristics using the monthly database of the GPCC [99,100] is consistent with the GMFD results, with the greatest rainfall increase in the Colombian Caribbean during MJ (see Appendix B, left column in Figure A2) and more extensive anomalies in ND, covering the northern, central, and western parts of the country (right column in Figure A2).
Severe local rainfall events and environmental changes can often increase the incidence of hydrological risks, such as floods and landslides, especially in tropical regions, where an increase in the intensity and frequency of rainfall extremes in the coming decades is expected [24,75,[101][102][103][104].
Here, rainfall extremes were analyzed using the R20 mm index and SDII.
The differences in composite in MJ for the R20 mm show no significant increases in the number of very rainy days in the ARB (middle column in Figure 7), or increases associated with the SDII (Figure 7). This result agrees with the lack of significant differences found in the t-test of the means before and after 1994 for the mean precipitation in the ARB shown in Section 3.2. However, in ND, both extreme rainfall indices describe significant increases in most of Colombia, especially in the Caribbean, Pacific, and Andean regions (middle and right columns in Figure 8). These results show that the upward change of the mean streamflow in the ARB at Bellavista in ND coincides with the increases in intensity and frequency of extreme rainfall events during the WAMO.

Conclusions
We evaluated how the low-frequency background associated with the AMO modifies the anomaly patterns of some atmospheric and oceanic variables associated with extreme streamflow events in the Atrato River basin (ARB), particularly during the two bimonthly streamflow peaks in May-June (MJ) and November-December (ND). Through a Student' t-test, the mean rainfall and mean streamflow were compared considering the cold AMO (CAMO) and warm AMO (WAMO), which refer to the 1965-1994 and 1995-2015 periods, respectively. Our results showed a significant increase in the streamflow means at Bellavista in MJ and ND and in the rainfall means in ND from the first to the second period (Table 1, Figure 4).

Conclusions
We evaluated how the low-frequency background associated with the AMO modifies the anomaly patterns of some atmospheric and oceanic variables associated with extreme streamflow events in the Atrato River basin (ARB), particularly during the two bimonthly streamflow peaks in May-June (MJ) and November-December (ND). Through a Student' t-test, the mean rainfall and mean streamflow were compared considering the cold AMO (CAMO) and warm AMO (WAMO), which refer to the 1965-1994 and 1995-2015 periods, respectively. Our results showed a significant increase in the streamflow means at Bellavista in MJ and ND and in the rainfall means in ND from the first to the second period (Table 1, Figure 4).
Composite analyses of the SST and SLP anomalies based on the streamflow extreme events in MJ and ND during each AMO phase reproduce the main anomaly patterns of these variables previously documented to be related to the AMO [38,39,41,86]. The rainfall anomaly composites show a more defined pattern, with significant positive anomalies in extensive areas during the WAMO phase, which are indeed modulated by the LN events in the tropical Pacific, especially during ND. In other words, rainfall anomalies related to the cooling of the tropical Pacific during WAMO (1995WAMO ( -2015 are more pronounced and occupy more extensive areas than those during the CAMO phase . Our findings show that the Atlantic Ocean plays an important role in modulating the rainfall and streamflow variability in the ARB. At the low frequency time scale, the SST and SLP anomalies in this Oceanic sector create the east-west SST and SLP gradients between the tropical Atlantic and eastern Pacific Oceans, which during the WAMO phase favors the La Niña event establishment [77]. Physical explanation is that the changes in the tropical Atlantic SST modify the inter-basin (tropical Atlantic and tropical eastern Pacific) pressure gradient, creating an atmospheric bridge that transmits the SST information between the basins [94,105]. The warm tropical Atlantic during the WAMO intensifies the Walker circulation [106] and confirms the previously discussed AMO forcing during the last two decades (1995-2015) [93,[107][108][109]. McGregor et al. [94] indicated that the tropical Pacific Ocean shows changes in the average state at a multidecadal scale resulting from the AMO, particularly in the tropics, along with interannual and decadal scales [94,105,109].
The strengthening of the Walker circulation after 1994 played a fundamental role by causing upward motion in northern and northwestern South America (SA) (Figure 6), which contributed to the positive rainfall anomalies in the ARB. Significant increases in intensities (SDII) and frequencies (R20 mm) of extreme rainfall events in ND were found in the ARB during the WAMO (Figure 8). However, in MJ there was no significant increase in the extreme rainfall in the basin (middle and right columns in Figure 7).
Our analyses suggest that the low-frequency background associated with the AMO phases modulated the rainfall and streamflow in the ARB during the last five decades. However, the changes in the streamflow and rainfall in the ARB might also be partially related to the global warming in the 1965-2015 period. The results here improved our understanding on the hydrometeorological variability in the ARB. We believe that our results might contribute to the assessment of climate resilience and the monitoring of hydrometeorological threats. Changes in land use may alter the basin responses to rainfall, however, this issue is outside the scope of this study, and future studies are needed to evaluate possible changes.

Conflicts of Interest:
The authors declare no conflict of interest. The founding sponsors had no role in the design, analysis and interpretation of data, in the writing manuscript, or in the decision to publish the results.

Appendix A
The streamflow data in Bellavista, San Antonio and Belén for the 1965-2015 period were provided by the Institute of Hydrology, Meteorology and Environmental Studies (IDEAM) (Figure 1a). These data were subjected to the principal component analysis (PCA) and the correlations among them. For PCA, the methodology proposed by Jolliffe and Cadima [110] was used. For the monthly analysis, the first mode explains 92% of the total monthly streamflow variance. For the bimonthly analyses, the MJ first mode and ND first mode explain 79% and 89% of the corresponding total bimonthly streamflow variances ( Figure A1). These results indicate that the three stations are highly correlated with each other and that the first mode of each analysis incorporates the largest amount of information from the original variables. In addition, the monthly correlations between Bellavista and the other stations were greater than 0.76 (p < 0.05) at monthly scale (Table A1). The analyses here ensure that the streamflow time series at Bellavista can be used as representative of the other two stations in the ARB.

Appendix A
The streamflow data in Bellavista, San Antonio and Belén for the 1965-2015 period were provided by the Institute of Hydrology, Meteorology and Environmental Studies (IDEAM) ( Figure  1a). These data were subjected to the principal component analysis (PCA) and the correlations among them. For PCA, the methodology proposed by Jolliffe and Cadima [110] was used. For the monthly analysis, the first mode explains 92% of the total monthly streamflow variance. For the bimonthly analyses, the MJ first mode and ND first mode explain 79% and 89% of the corresponding total bimonthly streamflow variances ( Figure A1). These results indicate that the three stations are highly correlated with each other and that the first mode of each analysis incorporates the largest amount of information from the original variables. In addition, the monthly correlations between Bellavista and the other stations were greater than 0.76 (p < 0.05) at monthly scale (Table A1). The analyses here ensure that the streamflow time series at Bellavista can be used as representative of the other two stations in the ARB.