Teleconnections between Monthly Rainfall Variability and Large-Scale Climate Indices in Southwestern Colombia

: Given that the analysis of past monthly rainfall variability is highly relevant for the adequate management of water resources, the relationship between the climate-oceanographic indices, and the variability of monthly rainfall in Southwestern Colombia at di ﬀ erent time scales was chosen as the research topic. It should also be noted that little-to-no research has been carried out on this topic before. For the purpose of conducting this research, we identiﬁed homogeneous rainfall regions while using Non-Linear Principal Component Analysis (NLPCA) and Self-Organizing Maps (SOM). The rainfall variability modes were obtained from the NLPCA, while their teleconnection in relation to the climate indices was obtained from Pearson’s Correlations and Wavelet Transform. The regionalization process clariﬁed that Nariño has two regions: the Andean Region (AR) and the Paciﬁc Region (PR). The NLPCA showed two modes for the AR, and one for the PR, with an explained variance of 75% and 48%, respectively. The correlation analyses between the ﬁrst nonlinear components of AR and PR regarding climate indices showed AR high signiﬁcant positive correlations with Southern Oscillation Index (SOI) index and negative correlations with El Niño / Southern Oscillation (ENSO) indices. PR showed positive ones with Niño1 + 2, and Niño3, and negative correlations with Niño3.4 and Niño4, although their synchronous relationships were not statistically signiﬁcant. The Wavelet Coherence analysis showed that the variability of the AR rainfall was inﬂuenced principally by the Niño3.4 index on the 3–7-year inter-annual scale, while PR rainfall were inﬂuenced by the Niño3 index on the 1.5–3-year inter-annual scale. The El Niño (EN) events lead to a decrease and increase in the monthly rainfall on AR and PR, respectively, while, in the La Niña (LN) events, the opposite occurred. These results that are not documented in previous studies are useful for the forecasting of monthly rainfall and the planning of water resources in the area of study. W.A.-M.; preparation—T.C., W.A.-M., W.L.C., Y.C.-E.; W.A.-M., visualization—T.C., W.A.-M., Y.C.-E., E.C.-B.; supervision—W.A.-M.,

. The geographic location of the study area and the distribution of the rainfall gauge stations.

Rainfall Dataset
We used monthly rainfall time series, from 1983 to 2016, which included the forty-four rainfallgauge stations located in different zones in Nariño, as provided by Instituto de Hidrología, Meteorología y Estudios Ambientales (IDEAM) of Colombia. Table 1 presents the main features and descriptive statistical attributes of each rainfall gauge station. The monthly rainfall series presented missing data in the percentage, also specified in Table 1. Here, we estimated this data using nonlinear principal component analysis (NLPCA), as suggested by Scholz et al. [29]. Refer to Canchala et al. [28] for further details regarding how missing data in Nariño using this non-linear approach was estimated. After completing the procedure for each of the series, we calculated the accuracy using the Root Mean Square Error (RMSE)

Rainfall Dataset
We used monthly rainfall time series, from 1983 to 2016, which included the forty-four rainfall-gauge stations located in different zones in Nariño, as provided by Instituto de Hidrología, Meteorología y Estudios Ambientales (IDEAM) of Colombia. Table 1 presents the main features and descriptive statistical attributes of each rainfall gauge station. The monthly rainfall series presented missing data in the percentage, also specified in Table 1. Here, we estimated this data using non-linear principal component analysis (NLPCA), as suggested by Scholz et al. [29]. Refer to Canchala et al. [28] for further details regarding how missing data in Nariño using this non-linear approach was estimated.
After completing the procedure for each of the series, we calculated the accuracy using the Root Mean Square Error (RMSE) where r is the reconstructed value, o indicates the observed value, and n refers to sample size.

Climate-Oceanographic Indices
Additionally, we selected eight monthly atmospheric and ocean time-series. This selection was performed while taking into account the results of previous research studies on teleconnections between the hydro-climatology of Colombia, and large-scale climate indices that have identified relationships with Southern Oscillation Index (SOI), Multivariate ENSO Index (MEI), Oceanic Niño Index (ONI), and SST in regions of the TPO (Niño1 + 2, Niño3, Niño3.4, Niño4) and Pacific Decadal Oscillation (PDO) [9][10][11][12][22][23][24]27,[30][31][32]. These indices are available in the National Oceanic and Atmospheric Administration (NOAA) (https://www.esrl.noaa.gov/psd/data/climateindices/list/). SOI is a measure of the large-scale oscillations in air pressure that take place between the west and east of the TPO, and is estimated as the difference in the Sea Level Pressure (SLP) anomalies between Tahiti and Darwin. MEI is an index obtained from six variables over the TPO:SLP, winds, SST, air temperature, and cloudiness [33]. ONI corresponds to the three-month moving average of SST anomalies on the Niño3.4 region [34]. Niño1 + 2, Niño3, Niño3.4, and Niño4 were calculated with the average of SST in the boxes 0-10 • S, 90

Regionalization of Monthly Rainfall
In this study, we established homogeneous rainfall regions in Nariño through the combination of two non-linear methods: NLPCA, and Self-Organized Feature Map (SOM), as efficient approaches to identify clusters with similar features of rainfall. NLPCA is a multi-layer feedforward network with a self-associative topology, also known as auto-encoder [29], which reduces the dimensionality of the dataset extracting its main features [36]. SOM is an Artificial Neural Network (ANN) based unsupervised clustering algorithm, introduced by Kohonen [37], which approximates the probability density function of the input data. It preserves the neighborhood properties and local resolution of the input space proportional to the data distribution [38,39]. Due to the robustness of SOM, its applications are focused on clustering and the modeling of hydrological datasets, such as precipitation, streamflow, runoff, and others [40]. Lin and Chen [41] and Chen et al. [42] have evidenced the advantage of SOM over conventional regionalization methods. Thus, the input data with 44-gauge stations into a temporal series of 408 months (34 years), required the use of NLPCA to reduce the dimensionality of the dataset to five non-linear components using a network with a (408-200-25-5) topology. These five nonlinear components were then used as input of SOM to capture a profile of the homogeneous areas and obtain a classification of gauge stations. The grid obtained from SOM in the output layer corresponds to 25 nodes (5 × 5 cells), which was enough to split the data-gauge stations. Finally, we validated the identified clusters through Discriminant Analysis (DA), as developed by Fisher [43]; more details about this test are described in Phalm et al. [44]. We used the Kolmogorov-Smirnov (K-S) test to measure the homogeneity of the identified sub-regions. This non-parametric test allows us to infer that the subjacent distribution of two datasets has significant differences [45,46].

Nonlinear Principal Component Analysis
The NLPCA is a non-linear generalization of principal component analysis (PCA) [47]. Given that the NLPCA generalizes the components, including lines and curves, it, therefore, depicts the structure of the data in curved subspaces [48]. NLPCA, as developed by Hsieh. [47] and Scholz. [48], uses the multi-layer perceptron of an auto-associative topology, known as a bottleneck (or auto-encoder). Figure 2 shows the neural network (NN) model for calculating NLPCA. topology. This diagram shows 3 "hidden" layers of "neurons" (indicated by circumferences) inserted between the input layer on the left and the output layer ̂ on the right. h(x) is the encoding layer, u is the bottleneck layer and h(u) is the decoding layer. 1, 2, 3, and 4 are the transfer functions. Adapted from Kenfack et al. [49].
The dimensions of x and h(x) are n and m, respectively. Where is the input column vector of length n and m is the number of hidden neurons in the encoding and decoding layers for u. More details about NLPCA are available in Scholz [50] and Hsieh [47].
NLPCA is a non-linear statistical method that is employed to analyze hydro-climatological, meteorological, and oceanographic data, such as rainfall, streamflow, or oceanographic variables [28,36,47,49,[51][52][53]. It has been widely used in recent years due to its advantages (over PCA) in the maximization of the explained global variance [47], most notably when non-linear processes, such as the rainfall distribution, rainfall clustering, atmospheric circulation regimes, among others, are involved [53,54]. Its objective is to obtain some non-linear components that explain the main modes of the variability of the original dataset that offer the least loss of information. This research uses the Non-linear PCA toolbox (http://www.nlpca.org/matlab.html) to obtain the hierarchically-ordered features by sequentially training and calculating the explained variance of each of the NLPCs. The NLPCA method helps to establish the dominant modes of variability of the monthly rainfall of Nariño in each of the clusters previously identified. The best architecture is based on the best performance, in terms of the highest percentage of explained variance.
Before performing monthly NLPCA for each group, we turned the rainfall values into anomalies by subtracting each month's mean.

Teleconnections
We applied the cross-correlations through Pearson Correlations, using a 12-month maximum lag, and taking into account that the climate indices precede the monthly rainfall, in order to measure the degree of correlation between the main modes of variability of monthly rainfall and the climateoceanographic indices, as well as the persistence of the relationship. The student t-test helped to define the statistical significance with a confidence value of 95% ( = 0.05).
Once the climatic indices with the highest incidence on the NLPCs were determined, we estimated the energy spectrum by wavelet transform to verify the spatio-temporal variability of rainfall, and the cross-wavelet analysis to compare two-time series. The hydro-climate field widely used the Wavelet Transform Methodology for different purposes, such as to: analyze the ENSO influence over the spatio-temporal variability of droughts [12]; study rainfall extremes and its relationship with climate indices [55]; identify teleconnections between monthly streamflow and climate indices [11,27]; and, describe the variability in annual streamflows [56]. In this regard, we used Morlet's Wavelet to obtain the time-frequency representation of NLPCs following the methodology presented by Grinsted et al. [57] and the computational procedure proposed by topology. This diagram shows 3 "hidden" layers of "neurons" (indicated by circumferences) inserted between the input layer x on the left and the output layerx on the right. h(x) is the encoding layer, u is the bottleneck layer and h(u) is the decoding layer. φ1, φ2, φ3, and φ4 are the transfer functions. Adapted from Kenfack et al. [49].
The dimensions of x and h(x) are n and m, respectively. Where x is the input column vector of length n and m is the number of hidden neurons in the encoding and decoding layers for u. More details about NLPCA are available in Scholz [50] and Hsieh [47].
NLPCA is a non-linear statistical method that is employed to analyze hydro-climatological, meteorological, and oceanographic data, such as rainfall, streamflow, or oceanographic variables [28,36,47,49,[51][52][53]. It has been widely used in recent years due to its advantages (over PCA) in the maximization of the explained global variance [47], most notably when non-linear processes, such as the rainfall distribution, rainfall clustering, atmospheric circulation regimes, among others, are involved [53,54]. Its objective is to obtain some non-linear components that explain the main modes of the variability of the original dataset that offer the least loss of information. This research uses the Non-linear PCA toolbox (http://www.nlpca.org/matlab.html) to obtain the hierarchically-ordered features by sequentially training and calculating the explained variance of each of the NLPCs. The NLPCA method helps to establish the dominant modes of variability of the monthly rainfall of Nariño in each of the clusters previously identified. The best architecture is based on the best performance, in terms of the highest percentage of explained variance.
Before performing monthly NLPCA for each group, we turned the rainfall values into anomalies by subtracting each month's mean.

Teleconnections
We applied the cross-correlations through Pearson Correlations, using a 12-month maximum lag, and taking into account that the climate indices precede the monthly rainfall, in order to measure the degree of correlation between the main modes of variability of monthly rainfall and the climate-oceanographic indices, as well as the persistence of the relationship. The student t-test helped to define the statistical significance with a confidence value of 95% (α = 0.05).
Once the climatic indices with the highest incidence on the NLPCs were determined, we estimated the energy spectrum by wavelet transform to verify the spatio-temporal variability of rainfall, and the cross-wavelet analysis to compare two-time series. The hydro-climate field widely used the Wavelet Transform Methodology for different purposes, such as to: analyze the ENSO influence over the spatio-temporal variability of droughts [12]; study rainfall extremes and its relationship with climate indices [55]; identify teleconnections between monthly streamflow and climate indices [11,27]; and, describe the variability in annual streamflows [56]. In this regard, we used Morlet's Wavelet to obtain the time-frequency representation of NLPCs following the methodology presented by Grinsted et al. [57] and the computational procedure proposed by Torrence and Compo [58]. To avoid spectral leakage and implement the Wavelet Transform, the time-series was filled with zeros to twice the data length. The wavelet function in each scale is standardized to get the energy unit, ensuring the ability to compare the wavelets with others. For more detailed wavelet theory, refer to Torrence and Compo [58]. The Global Wavelet Spectrum (GWS) allows for a better visualization of the distribution of power linked to the frequencies in the climate signals and ensures a better detection of fluctuation/oscillation. GWS is the time mean of all local Wavelet Power Spectrum (WPS) and it is explained with an equation (22) in Torrence and Compo [58].
The Wavelet Coherence (WTC) allows us to compare two-time series (NLPCs vs. climate-indices), due to the fact that there is a correlation coefficient that is estimated in time and space frequency. According to Torrence and Webster [59], X(t) and Y(t) series are given, with Wavelet Transforms, W X (t, s) and W Y (t, s), the WTC can be defined as where is a smoothing operator, W XY (t, s) = W X (t, s)W Y * (t, s) is the cross-wavelet spectrum, (*) indicates the conjugated complex, R 2 n can range from 0 to 1. 0 designates no-correlation between the two-series that are compared, and 1 indicates a perfect correlation. Monte Carlo methods help to estimate the statistical significance of the WTC. More details regarding the estimation of WTC's are available in Grinsted et al. [57]. Figure 3 shows the recovery of missing data from the Remolino (REM) gauge station, which has the most missing data (see Table 1). The results of the estimation of all stations' missing data are available in Canchala et al. [28]. The best topology to recover data was with a network (45-44-45). The RMSE in the estimation of missing data corresponds to 9.8 mm/month. Torrence and Compo [58]. To avoid spectral leakage and implement the Wavelet Transform, the timeseries was filled with zeros to twice the data length. The wavelet function in each scale is standardized to get the energy unit, ensuring the ability to compare the wavelets with others. For more detailed wavelet theory, refer to Torrence and Compo [58]. The Global Wavelet Spectrum (GWS) allows for a better visualization of the distribution of power linked to the frequencies in the climate signals and ensures a better detection of fluctuation/oscillation. GWS is the time mean of all local Wavelet Power Spectrum (WPS) and it is explained with an equation (22) in Torrence and Compo [58]. The Wavelet Coherence (WTC) allows us to compare two-time series (NLPCs vs. climateindices), due to the fact that there is a correlation coefficient that is estimated in time and space frequency. According to Torrence and Webster [59], X(t) and Y(t) series are given, with Wavelet Transforms, ( , ) and ( , ), the WTC can be defined as

Missing Data Estimation
where 〈〉 is a smoothing operator, ( , ) = ( , ) * ( , ) is the cross-wavelet spectrum, (*) indicates the conjugated complex, 2 can range from 0 to 1. 0 designates no-correlation between the two-series that are compared, and 1 indicates a perfect correlation. Monte Carlo methods help to estimate the statistical significance of the WTC. More details regarding the estimation of WTC's are available in Grinsted et al. [57]. Figure 3 shows the recovery of missing data from the Remolino (REM) gauge station, which has the most missing data (see Table 1). The results of the estimation of all stations' missing data are available in Canchala et al. [28]. The best topology to recover data was with a network (45-44-45). The RMSE in the estimation of missing data corresponds to 9.8 mm/month.

Regionalization of Monthly Rainfall
The results of the non-linear approach for rainfall regionalization using NLPCA and SOM showed that Nariño has two clusters: the Andean Region (AR) and the Pacific Region (PR), each one with thirty-three and eleven-gauge stations, respectively (See Figure 4a). The Monopamba (MON) gauge station had to be removed from the analysis due to it is producing noise in the SOM study.
The validation of the homogeneity of climatic regions, using the Fisher's Discriminant Analysis [43], showed a 100% AR and PR classification accuracy. Besides, the K-S test allowed for us to verify the statistical significance of the identified regions, by assuming the null hypothesis that the

Regionalization of Monthly Rainfall
The results of the non-linear approach for rainfall regionalization using NLPCA and SOM showed that Nariño has two clusters: the Andean Region (AR) and the Pacific Region (PR), each one with thirty-three and eleven-gauge stations, respectively (See Figure 4a). The Monopamba (MON) gauge station had to be removed from the analysis due to it is producing noise in the SOM study.
Water 2020, 12, 1863 8 of 20 distribution is the same [45,46]. This test showed a p-value < 0.05, indicating that the precipitation regimes in the regions are different. Figure 4b presents the spatial distribution of the forty-four gauge stations, organized according to the regionalization results. The main feature of rainfall in AR is that it records a bimodal annual cycle with an average monthly rainfall of about 130 mm/month ( Figure  4c), whereas the PR shows a unimodal yearly cycle with an average monthly rainfall of about 350 mm/month (Figure 4d). For a better understanding of the spatial variability of rainfall anomalies in the two regions, refer to Figure 5, where you can find the first spatial pattern that was obtained from NLPCA that clearly shows two different rainfall patterns over the AR and PR.  The validation of the homogeneity of climatic regions, using the Fisher's Discriminant Analysis [43], showed a 100% AR and PR classification accuracy. Besides, the K-S test allowed for us to verify the statistical significance of the identified regions, by assuming the null hypothesis that the distribution is the same [45,46]. This test showed a p-value < 0.05, indicating that the precipitation regimes in the regions are different. Figure 4b presents the spatial distribution of the forty-four gauge stations, organized according to the regionalization results. The main feature of rainfall in AR is that it records a bimodal annual cycle with an average monthly rainfall of about 130 mm/month (Figure 4c), whereas the PR shows a unimodal yearly cycle with an average monthly rainfall of about 350 mm/month (Figure 4d). For a better understanding of the spatial variability of rainfall anomalies in the two regions, refer to Figure 5, where you can find the first spatial pattern that was obtained from NLPCA that clearly shows two different rainfall patterns over the AR and PR. Water 2020, 12, 1863 9 of 20
In SA, the establishment of an LN (EN) episode occurred with decreases (increases) of rainfall over the west coast of Ecuador, southern and southeastern of the continent, and increases (decreases) of them over the north and northeastern of SA [61,62]. Moreover, in the western, central, and northern
progressive flooding [8,63,64]. In contrast, the pattern EN is related to an increase in the average air temperature, a decrease in rainfall, and a reduction in the average flow of the rivers [8,65]. This behavior is a characteristic of large areas of the Caribbean, Andean, and North Pacific Regions. Moreover, in the South Pacific Region, Southwestern Colombia, the Amazon, and some areas of the Piedmont Plains, the establishment of EN (LN) events is associated with increases (decreases) in rainfall [7]. Hence, the results registered in AR (Figure 6a,b) and PR (Figure 6c) are consistent with the above-mentioned, given that the rainfall in these two sub-regions recorded opposite effects when EN and LN events occurred. As an example, the monthly rainfall in AR, decreases when an EN event is registered, while the PR increases. In contrast, when the LN event occurred, the rainfall in AR and PR increases and decreases, respectively. In Colombia, the 1997-1998 EN event is considered to have had the highest intensity and spatial amplitude with devastating effects on different economic sectors. These generated phenomena, such as heatwaves, droughts, forest fires, as well as diminished rainfall and streamflows in some regions; while others recorded more torrential rainfalls, landslides, and floods [66]. In this regard, the registered behavior for AR (PR) NLPCs showed negative (positive) anomalies during the 1997-1998 EN event (See Figure 6). Meanwhile, the effects of this EN event in PR were similar to those that occurred in 1997-1998 in Ecuador's coastal region (a region that borders with the PR of Nariño), which produced significant positive rainfall anomalies, affecting different economic sectors [67]. These ENSO events have been cataloged by the NOAA as a strong event.
According to Quishpe-Vásquez et al. [68], during the EN events, the trade winds weakens along the equator; therefore, the Walker circulation also weakens as atmospheric pressure in western OPT increases, while the atmospheric pressure in the Eastern Pacific decreases. It caused a displacement of the warm waters of the Central Pacific to the east, and the weakening of the Humboldt current In SA, the establishment of an LN (EN) episode occurred with decreases (increases) of rainfall over the west coast of Ecuador, southern and southeastern of the continent, and increases (decreases) of them over the north and northeastern of SA [61,62]. Moreover, in the western, central, and northern territories of Colombia, LN is characterized by high and intense rainfall, increased streamflows, and progressive flooding [8,63,64]. In contrast, the pattern EN is related to an increase in the average air temperature, a decrease in rainfall, and a reduction in the average flow of the rivers [8,65]. This behavior is a characteristic of large areas of the Caribbean, Andean, and North Pacific Regions. Moreover, in the South Pacific Region, Southwestern Colombia, the Amazon, and some areas of the Piedmont Plains, the establishment of EN (LN) events is associated with increases (decreases) in rainfall [7].
Hence, the results registered in AR (Figure 6a,b) and PR (Figure 6c) are consistent with the above-mentioned, given that the rainfall in these two sub-regions recorded opposite effects when EN and LN events occurred. As an example, the monthly rainfall in AR, decreases when an EN event is registered, while the PR increases. In contrast, when the LN event occurred, the rainfall in AR and PR increases and decreases, respectively. In Colombia, the 1997-1998 EN event is considered to have had the highest intensity and spatial amplitude with devastating effects on different economic sectors. These generated phenomena, such as heatwaves, droughts, forest fires, as well as diminished rainfall and streamflows in some regions; while others recorded more torrential rainfalls, landslides, and floods [66]. In this regard, the registered behavior for AR (PR) NLPCs showed negative (positive) anomalies during the 1997-1998 EN event (See Figure 6). Meanwhile, the effects of this EN event in PR were similar to those that occurred in 1997-1998 in Ecuador's coastal region (a region that borders with the PR of Nariño), which produced significant positive rainfall anomalies, affecting different economic sectors [67]. These ENSO events have been cataloged by the NOAA as a strong event.
According to Quishpe-Vásquez et al. [68], during the EN events, the trade winds weakens along the equator; therefore, the Walker circulation also weakens as atmospheric pressure in western OPT increases, while the atmospheric pressure in the Eastern Pacific decreases. It caused a displacement of the warm waters of the Central Pacific to the east, and the weakening of the Humboldt current [69]. Furthermore, the ITCZ movement to the south caused heavy rainfall on the coast of Ecuador [70,71]. The opposite effect can occur with LN events; for example, the 2010-2011 LN event in Nariño showed positive (negative) anomalies for AR (PR). The positive anomalies in AR match anomalies registered in a significant part of Colombia, which recorded extreme precipitation and river discharges that affected approximately 9% of the country's population [65]. In contrast, the negative anomalies in PR match anomalies registered on neighboring regions of Colombia, specifically in the Ecuadorian coast where the drought is associated to the LN event [72].
These results have shown the importance of analyzing the monthly rainfall variability to regional scale and study its relationships with ENSO indices, when considering that, despite there being rainfall regionalization studies throughout the country, such as those developed by Guzman et al. [73], and Estupiñan [74], which coincide with the regionalization performed here, it is the first time that it found significant differences in the relationships between ENSO events and the Southwestern Colombia rainfall, even at very short distances between these sub-regions. Therefore, in the next section, we consider it relevant to study the teleconnections with the different ENSO indices to identify the different time scales required to improve understanding of the findings found so far.
We used the wavelet transform analyses to understand the variability of monthly rainfall in terms of time and frequency with the purpose of studying the teleconnections between AR and PR rainfall with climate-oceanographic indices in their entirety. The Wavelet Analysis applied to NLPC1-AR (Figure 8a Figure 10 depicts the Wavelet Coherence results for AR rainfall and the climate-oceanographic indices that best explain the variability of rainfall. The coherence wavelet between NLPC1-AR and Niño1 + 2 and Niño3 (ONI, MEI, Niño3.4, and Niño4) indices showed similar relationships that were mainly focused at a time scale of 3-8-(3-7-) years from 1989 to 2011 (1989 to 2011), with a phase difference of -145° to 180° (-135° to 180°) (refer to Figure 10a,b, respectively). The -145° and 135° phase relationship for the 3-8-year and 3-7-year interannual scale shows that the EN (LN) event either led to dry (wet) conditions for approximately 4-9 months and 5-11 months, respectively. The 180° phase difference that shows the negative (positive) rainfall anomalies have a synchronous relationship with the mature phase of the EN (LN) events. Besides, there are significant coherences between the SOI index and NLPC1-AR for a 3-7-year inter-annual scale , also present in the same 0° to 45° phase (Figure 10c). A 45° phase relationship indicates that the negative (positive) rainfall anomalies were either led by negative (positive) SOI index by approximately 5-11 months, whereas 0° indicates a synchronous relationship between both series. Meanwhile, the Wavelet Coherence between the NLPC2-AR and Niño1 + 2, Niño3.4, Niño4, and ONI indices showed similar relationships (Figure 10d only shows El Niño1 + 2). The main ones were at a time scale of 1.5-3-years from 1995 to 2003 with a phase difference of 160° to 180°, establishing a phase difference around 1-2 months. Therefore, this implies that EN (LN) event preceded dry (wet) conditions by approximately 1-2 months. Figure 11a,b show the Wavelet Coherence between the NLPC1-PR with Niño3 and SOI indices, respectively. Figure 11a showed significant relationships on the 1.5-3-year interannual scale from 1995 to 2003, with a phase relationship of 10° to 45°, indicating an approximate phase difference of 1-5 months. This implies that there was an increase (decrease) in the SST in the Niño3 region, which led to positive (negative) rainfall anomalies in this region. Furthermore, it shows an 8-12-year decadal-scale from 1995 to 2006, with a −45° phase relationship, indicating that PR rainfalls occur after fluctuations in Niño3 by around 12-18 months. Furthermore, the relationship between SOI index and NLPC1-PR (Figure 11b) did not help to establish a significant coherence wavelet, with the exception of the period of 1992 to 1995, where a phase relationship of 0° to 10° at a 0-1.5-year interannual scale was evidenced, which implied lags of 0-1 month between them.  Figure 10 depicts the Wavelet Coherence results for AR rainfall and the climate-oceanographic indices that best explain the variability of rainfall. The coherence wavelet between NLPC1-AR and Niño1 + 2 and Niño3 (ONI, MEI, Niño3.4, and Niño4) indices showed similar relationships that were mainly focused at a time scale of 3-8-(3-7-) years from 1989 to 2011 (1989 to 2011), with a phase difference of -145 • to 180 • (-135 • to 180 • ) (refer to Figure 10a,b, respectively). The -145 • and 135 • phase relationship for the 3-8-year and 3-7-year interannual scale shows that the EN (LN) event either led to dry (wet) conditions for approximately 4-9 months and 5-11 months, respectively. The 180 • phase difference that shows the negative (positive) rainfall anomalies have a synchronous relationship with the mature phase of the EN (LN) events. Besides, there are significant coherences between the SOI index and NLPC1-AR for a 3-7-year inter-annual scale , also present in the same 0 • to 45 • phase (Figure 10c). A 45 • phase relationship indicates that the negative (positive) rainfall anomalies were either led by negative (positive) SOI index by approximately 5-11 months, whereas 0 • indicates a synchronous relationship between both series. Meanwhile, the Wavelet Coherence between the NLPC2-AR and Niño1 + 2, Niño3.4, Niño4, and ONI indices showed similar relationships (Figure 10d only shows El Niño1 + 2). The main ones were at a time scale of 1.5-3-years from 1995 to 2003 with a phase difference of 160 • to 180 • , establishing a phase difference around 1-2 months. Therefore, this implies that EN (LN) event preceded dry (wet) conditions by approximately 1-2 months.  Overall, the cross-correlation analysis (Figure 7) helped to establish that the rainfall of the AR represented in the primary variability mode (NLPC1-AR), has a meaningful teleconnection with oceanic macroclimatic indices linked to the ENSO phenomenon. In contrast, the PR rainfall of  Overall, the cross-correlation analysis ( Figure 7) helped to establish that the rainfall of the AR represented in the primary variability mode (NLPC1-AR), has a meaningful teleconnection with oceanic macroclimatic indices linked to the ENSO phenomenon. In contrast, the PR rainfall of represented in NLPC1-PR has a weak teleconnection with climate indices that were studied here. The weak direct relationship is mainly associated to the sea surface temperature of the Niño3 region, showing that the ENSO phenomenon does not strongly affect the variability of rainfall in the coastal region of Nariño. These results from the cross-wavelet analysis (Figures 10 and 11) showed that there is a strong (weak) relationship inverse (direct) between AR (PR) rainfall and the SST in Niño3.4, and Niño4 (Niño3) on the 3-7-(1.5-3)-year inter-annual scale. Besides helping to establish that the warm (cold) phase of ENSO precedes dry(wet) conditions by approximately 5-11 months in AR, while the results also helped establish that the warm (cold) phase of ENSO precedes wet(dry) conditions by around 1-5 months in PR.
The results obtained for NLPC's of AR are consistent with those that were obtained by Montealegre [75] and Navarro et al. [76]. They reported a reduction (increase) in the monthly rainfalls linked to positive (negative) anomalies of the oceanic ENSO indices, mainly in the gauge stations of the Andean Zone, the northwest, and north regions of Colombia. These results are also coherent with those that were obtained in the studies of the relationship between the other hydro-climatic variables and large scale climate indices in Colombia [11,21,30,32]. In contrast, the results obtained for PR showed that the monthly rainfall in this region has a direct relationship with the oceanic climate indices, mainly Niño3 and Niño1 + 2. These results are consistent with those that were obtained by Montealegre [75].
They demonstrate a reduction (increase) in the rainfall of the Tumaco gauge station (located in PR) linked to negative (positive) values of these oceanic indices. The relationship between rainfall and oceanic climate indices registered in PR is similar to the connection found on the Ecuadorian coast; here, De Guennie et al. [67] showed a positive correlation between Niño1 + 2, and Niño3 and the rainfall anomalies at lag = 0.  Overall, the cross-correlation analysis (Figure 7) helped to establish that the rainfall of the AR represented in the primary variability mode (NLPC1-AR), has a meaningful teleconnection with oceanic macroclimatic indices linked to the ENSO phenomenon. In contrast, the PR rainfall of With respect to the results of the cross-correlation analysis of the lagged effect to 12 months of the climate-oceanographic indices on rainfall in AR (PR), a strong (weak) persistence of up to six to 12 (five to 12) months was evidenced (See Figure 7). The results were consistent with Navarro et al. [77]. They indicate a lagged effect of up to 6-9 months regarding the relationship between the rainfalls in Western Colombia and the ENSO indices. In the same way, Canchala et al. [11] showed that the teleconnections between streamflows of two large rivers of the Choco Biogeographic Colombian and large scale climate indices are persistent for nine months. Moreover, this relationship was positive (negative) with atmospheric (oceanic) indices. Unlike the NLPC1-PR, the maximum correlations for the NLPC1-AR were obtained with the Niño indices (1 + 2, 3, 3.4, and 4) synchronously. Besides, the Niño3.4 and Niño4 influences on the AR rainfall were more substantial than the Niño1 + 2 and Niño3, which remained for another 10 and eight months, respectively. The result is coherent with the work developed by Poveda and Mesa [63]. They identified that Colombian precipitation has a higher correlation with the SST in the Niño4 region in the Central Pacific than other Niño regions close to SA.
In contrast, for NLPC1-PR, the strongest correlations with ENSO indices were recorded for Niño3.4 and Niño4 in 8, 9, and 10-lags (Figure 7c). Therefore, the ENSO influence occurred for 8, 9, and 10 months after the increase or decrease of SST in TPO. However, the findings indicate that the highest synchronous correlations depend on Niño1 + 2 and Niño3, which were positive, indicating a direct relationship. The influence of Niño1 + 2 is similar to that exerted on the rainfall of the neighboring country (Ecuador). It is directly influenced by SST in the waters near the shore (Niño1 + 2) [70]. On the other hand, AR and PR showed positive relationships with the SOI index, and these were coherent with results that were registered by Ávila et al. [9], Canchala et al. [11], Poveda et al. [16], and Loaiza et al. [78]. The influence of the SOI index is higher in the AR than in the PR. Furthermore, the rainfall influence on the AR is almost synchronous (Lag 1), while the rainfall in PR is substantial in lags 8, 9, and 10.
Other findings show a statistically significant negative relationship between the PDO index and the NLPC1-AR with a maximum value for a nine-month lag (Figure 7a). A similar relationhip was recorded between the streamflows of the two basins of the Choco Biogeographic Colombian by Canchala et al. [11]. According to Labat [79] and Shi et al. [80], ocean-atmospheric oscillations altogether enforce a strong influence on hydro-climatology; hence, climatic variability is subject to the joining and teleconnections of these large-scale processes. The PDO can change the periodicity of ENSO episodes according to the phase, which favoured the occurrence of lower (higher) rainfall, while positive EN/PDO (negative LN/PDO) episodes occur [81].
In summary, the rainfall in two-sub-regions of Southwestern Colombia showed significant differences in their spatio-temporal variability and their teleconnections with climate-oceanographic indices. When considering that climate indices showed high coherence (statistically significant), further study for the adequate management of water resources is suggested, since these can facilitate the prediction of monthly rainfall several months in advance.

Conclusions
Using NLPCA, Pearson's Correlations, and Wavelet Analysis, in this research study, we were able to identify two sub-regions in southwestern Colombia: the Andean Region (AR) and Pacific Region (PR), in order to study their main variability time scales of the monthly rainfalls and their relationships with eight climate-oceanographic indices. Our findings are summarized, as follows: 1.
From the NLPCA, we found two main modes of variability in AR and one in PR that explained around 75% and 48% of the variance of the original datasets, respectively. The positive (negative) periods of NLPC1-AR and NLPC2-AR coincided with the LN (EN) events. In contrast, those positive (negative) periods of NLPC1-PR coincided with EN (LN) events, i.e., the variability of monthly rainfall in the Andean region is contrary to the variability of monthly rainfall in the Nariño coastal region. The variability of monthly rainfall in AR is similar to the variability of rainfall in most of the Colombian territory, and the variability of monthly rainfall in PR is similar to the variability of rainfall that is registered in the coastal region of the border country, Ecuador.

2.
The Pearson's Correlation between NLPCs of AR and PR and eight climate-oceanographic indices showed that the rainfall in AR has a direct (inverse) relationship with the SOI (ONI, MEI, Niño1 + 2, Niño3, Niño3.4, Niño4, and PDO) indices. Furthermore, the magnitude of the correlations was strongest with Niño3.4 and Niño4, and this influence remained for 10 and eight months, respectively. In this regard, the rainfall in PR showed positive synchronous correlations, but non-statistical significance with Niño1 + 2 and Niño3 indices, and an inverse relationship with Niño3.4 and Niño4 in 8, 9, and 10-lags. Overall, indices that were linked to the ENSO phenomenon showed strong teleconnections with the rainfall in AR, and weak ones in PR. The established results allow us to infer that the variability of rainfall in Southwestern Colombia influences the ENSO's inter-annual phenomenon its influence is stronger in east Nariño than on the west. These results are especially important for subsequent prediction analyses, given that the climate-oceanographic-indices with high synchronous and lagging correlation could become potential predictors of monthly rainfall in these sub-regions of Southwestern Colombia.

3.
Wavelet Coherence analysis confirmed the relationships between ENSO indices and rainfall of AR and PR, which showed that the variability of the rain in AR is strongly influenced by the ENSO indices, mainly by the SST in the regions Niño3.4 and Niño4 on the 3-7-year inter-annual scale. Furthermore, this analysis ratified the weak influence of the ENSO phenomenon on PR rainfall. However, we noticed that, in the 1.5-3-year inter-annual scale, there is an influence of the SST in the region Niño3. Overall, the rain in AR is associated to the SST in the east TPO (Niño3.4, Niño4, ONI, and MEI), where the warm (cold) phase of ENSO led to dry (wet) conditions for approximately 5-11 months. Meanwhile, the rainfall of PR is associated to the SST in the center of TPO (Niño3), where the warm (cold) phase of ENSO led to wet (dry) conditions for approximately 1-5 months. 4.
The main findings of this research study are an understanding of the monthly variability rainfall in AR and PR, in Southwestern Colombia concerning the teleconnections with eight climate indices (when considering the multi-scale relations). This study allowed for us to establish for the first time the contrary relationships between ENSO indices and monthly rainfall in the sub-regions studied, where the EN events lead to a decrease (increase) in the monthly rainfall on AR (PR), while the opposite occurred in the LN events. The ENSO influence on the rainfall of the AR is highly correlated, and it has similar behavior to the impact of ENSO over most of the Colombian territory. In contrast, the findings revealed that the ENSO influence on PR is weak, showing similarities with the ENSO influence over the rainfall in the neighboring country of Ecuador. These results show that future studies need to explore the ocean-atmospheric and regional circulation processes that explain the marked difference in the relationships between ENSO and the region's rainfall. Funding: The first author was supported by the Program for Strengthening Regional Capacities in Research, Technological Development and Innovation in the department of Nariño and the CEIBA foundation for doctoral studies. The third author was supported by the Universidad del Valle (Cali-Colombia). The authors thank the Universidad del Valle for financing the research project CI 21010, and Colciencias for funding the research project "Análisis de eventos extremos de precipitación asociados a variabilidad y cambio climático para la implementación de estrategias de adaptación en sistemas productivos agrícolas de Nariño".