Analysis of the Behavior of Daily Maximum Rainfall within the Department of Atlántico, Colombia

: In the Colombian Caribbean region, there are few studies that evaluated the behavior of one of the most commonly used variables in hydrological analyses: the maximum daily rainfall (P max-24h ). In this study, multiannual P max-24h time series from 19 rain gauges, located within the department of Atlántico, were analyzed to (a) determine possible increasing/decreasing trends over time, (b) identify regions with homogeneous behavior of P max-24h , (c) assess whether the time series are better suited under either a stationary or non-stationary frequency analysis, (d) generate isohyetal maps under stationary, non-stationary, and mixed conditions, and (e) evaluate the isohyetal maps by means of the calculation of areal rainfall (P areal ) in nine watersheds. In spite of the presence of both increasing and decreasing trends, only the Puerto Giraldo rain gauge showed a significant decreasing trend. Also, three regions (east, central, and west) with similar P max-24h behavior were identified. According to the Akaike information criterion test, 79% of the rain gauges showed better fit under stationary conditions. Finally, statistical analysis revealed that, under stationary conditions, the errors in the calculation of P areal were more frequent, while the magnitude of the errors was larger under non-stationary conditions, especially in the central–south region. simulated or true value better represents the areal rainfall value. For the central and east regions, negative values prevailed in most of the cases for either stationary or non-stationary conditions (only one value was above 0.5). Within the central region for stationary conditions, all return periods had negative values, while, for the non-stationary conditions, this behavior was seen in 80% of the cases. In the east region, four out of five return periods showed negative values for stationary conditions, and three out of five return periods showed negative values for non-stationary conditions. These results indicate that the average of the true value (mixed isohyetal maps) is a better predictor for these two regions.


Introduction
Globally, changes in the pattern of behavior of hydrometeorological variables (e.g., precipitation, temperature, runoff, relative humidity, etc.) are influenced, among others, by population growth, watershed land-use/land-cover (LULC) changes, and the increase in greenhouse gases emissions [1,2]. In Colombia, the Institute of Hydrology, Meteorology, and Environmental Studies (IDEAM, in Spanish) conducted several studies focused on evaluating the changes in the behavioral patterns of some hydrometeorological variables [3][4][5][6][7]. IDEAM [3] analyzed the annual average rainfall trend over the periods 2011-2040, 2041-2070, and 2071-2100 for the different departments (political and administrative territorial units) that compose the five regions of Colombia (Caribbean, Pacific, Andean, Orinoco, and Amazon). According to the study, the department of Atlántico (located within the Caribbean region) will experience an annual average rainfall decrease ranging from 7.39% through 11.26% during the 2011-2100 period. In addition, it was predicted that some of the municipalities located in the southeast of the department will be the most affected.
Such changes in the hydrological cycle can lead to (a) decreases in the water supply (both for human consumption and for the different sectors of the economy), (b) possible water supply cost increase, and (c) under-or oversized hydraulic structures for stormwater management [8]. Several studies analyzed the rainfall behavior within the department of Atlántico [9][10][11][12][13]. However, none The department of Atlántico has a warm and dry climate, with an average annual temperature of approximately 28 °C and maximum temperatures that can reach up to 40 °C. The average annual rainfall ranges from 500 mm to 1500 mm [15]. The rainfall regime has three seasons: dry (December to March), transition (late April to June), and rainy (August to early December) [22].
In this research, multi-annual series of maximum daily rainfall were used, from 19 rain gauges operated by IDEAM (Table 2), totaling 728 observations from 1940 to 2015 (records from 2016-2019 were not included as some rain gauges do not yet have the information available for those years). The rain gauges used in this study were selected under the following criteria: (a) time series with a minimum of 20 years of data, (b) exclusion of rain gauges with less than 25% of the rainfall information in any given year, (c) exclusion of rain gauges that did not have information from the months corresponding to the rainy season [23], and (d) elimination of outliers by means of the Water Resources Council method [24]. The selected rain gauges are shown in Table 2.

Methodology
After selecting the rain gauges that met the abovementioned criteria, the rainfall time series underwent further analysis. Figure 2 shows the steps that make up the methodology proposed in this research.

Trend Analysis
In this study, the monotonic trend detection was performed through the nonparametric tests of Mann-Kendall and Spearman's rho [25], with a significance level of 5%. Nonparametric testing has the advantage of being able to detect trends, independently of whether the data has a normal distribution or not, as with hydrometeorological variables [26]. In addition, the analysis was complemented by determining the magnitude of the trends' slopes identified in the Pmax-24h series using the Theil-Sen slope [27].

Mann-Kendall (MK) Test
The test considers a null hypothesis (H0) when no trend (to increase or decrease) exists and an alternative hypothesis (H1) that there is a trend. The calculation of Mann-Kendall's statistics S and Pre-processing of the P max-24h multiannual series standardized Z uses the following set of formulas (Equation (1)): where Xi and Xj represent the time series data in chronological order, n is the number of data points in the time series, tp is the number of links for the p-th value, and q is the number of links. Positive Z values (MK statistic) indicate the presence of increasing trends, and negative values indicate decreasing trends. If |Z| > Z1 − α/2, the null hypothesis is rejected, indicating a statistically significant trend. The critical value of Z1 − α/2 for a significance level of 5% is 1.96 [26]. That is, a trend is considered increasing or decreasing, at a significance level of 5%, only if ZMK is greater than |1.96|. Otherwise, it is considered trendless (constant).

Spearman's Rho (SR) Test
This test assumes that the data are independent and identically distributed. Null and alternative hypotheses are defined the same way as in the Mann-Kendall test [28]. The RSR and ZSR statistical variables are calculated using Equations (2) and (3).
where Di represents the i-th observation, i is the chronological order of the number, n is the number of observations, ZSR is the value of the t-student distribution with n − 2 degrees of freedom. Positive and negative ZSR values represent increasing and decreasing trends, respectively. If |ZSR|> t(n−2,1−α/2), the null hypothesis is rejected, demonstrating statistically significant trends [10,29].

Estimator of Theil-Sen Slope
This estimator allows determining the actual slope of the trends of a given time series. The principle of this estimator is based on the assumption that, when a series of data shows a linear trend, the median of the slope of the linear trend can be calculated using the slopes of several data points, and this value represents the slope of its trend (Equation (4)) [30,31].
where βTS represents the slope between the xj and xk points in the time series (which corresponds to the time points j and k, with j > k).

Delimiting Homogeneous Regions
The delimitation of regions with similar hydrometeorological conditions involved a statistical analysis of the data observed in the study area [32,33] by means of the regionalization method suggested by Hosking and Wallis [34]. This method defines statistical parameters similarly to the traditional L-moments, which, in turn, involves the calculation of the β estimators (Equations (5)-(8)).
where Xi represents the value of the Pmax-24h series, i is the rank of each data point arranged from highest to lowest, and n is the number of data points in the series of each rain gauge j. Subsequently, the L-moments (represented with λ) are obtained using Equations (9)-(12).

= ,
Finally, the dimensionless L-moments were calculated using Equations (13)- (15), where τ2 represents the coefficient of variation, τ3 is the asymmetry coefficient, and τ4 is the kurtosis coefficient.
Homogeneous regions were then formed via cluster analysis by means of the K-means method [35], which allows the dimensional L-moments to be related to the elevation and location parameters of each of the rain gauges. This way, the clusters sharing similar characteristics are detected so that the homogeneous regions can be defined. Additionally, varying the number of clusters in the Kmeans method helps with finding a geographically consistent configuration. Finally, the selected cluster configuration was reassessed using the methodology used by Hosking et al. [36] in order to corroborate the region's homogeneity.

Stationary and Non-Stationary Rainfall Frequency Analysis
The point rainfall values used for the Pmax-24h isohyetals were estimated via frequency analysis for both stationary (SC) and non-stationary (NSC) conditions for return periods of five, 10, 25, 50, and 100 years.
Based on the results obtained by González-Alvarez et al. [14] for the Colombian Caribbean region, generalized extreme value (GEV) [37] (Equation (16)) and Gumbel [38] (Equations (17)-(19)) were used for the stationary frequency analysis. The theoretical basis for these two cumulative distribution functions (CDF) can be widely found in the literature [39][40][41]. In Equations (17)-(19), k is the shape parameter, β is the mode (or location), and α is the scale (always greater than zero). GEV can take one of three extreme value (EV) distributions depending on the value of k: (a) if k equals zero, it takes the form of Type 1 EV (Gumbel); (b) if less than zero, it takes the form of Type 2 (Fréchet); (c) it takes the form of Type 3 (Weibull), if greater than zero. The Pmax-24h value selected for the stationary frequency analysis was the one that came from the function showing the best chi-square test [42] result.
The non-stationary frequency analysis was performed according to the methodology proposed by Obeysekera and Salas [1,43,44], which uses (a) the GEV function by varying the location parameter over time and maintaining the constant parameters of scale and shape (called GEVmu), and (b) a definition of the return period (Tr) according to the geometric distribution given by Equation (20), where Pj is the is the time-varying exceedance probability, and j represents the year to be projected [1,45]. The GEVmu function was already tested by Gonzalez-Alvarez et al. [45] in the Colombian Caribbean region, where a sensitivity analysis showed that varying the shape and/or scale parameters did not bring any improvement in the performance of either GEV or Gumbel distributions.
Subsequently, a linear trend model of each parameter (location, shape, and scale) was defined to estimate its value using a code programmed in the R software (Version 3.3.1, R Development Core Team, Auckland, New Zealand) with the library nsextremes [46].
After the Pmax-24h values were calculated for stationary and non-stationary conditions, the Akaike information criterion (AIC) goodness-of-fit test [47] was used to determine which of the two conditions (stationary and non-stationary) better represented the multiannual series for each of the rain gauges analyzed. The Pmax-24h values for SC, NSC, and the better of the two conditions were later used for the generation of the stationary, non-stationary, and mixed isohyetal maps, respectively (Section 3.4).

Generation of Isohyetals Maps
After obtaining the Pmax-24h values (for stationary and non-stationary conditions) for each of the rain gauges, three different types of Pmax-24h isohyetals (for return periods of five, 10, 25, 50, and 100 years) were generated (Table 3), using ArcGIS (Version 16.1, ESRI Inc., Redlands, CA, USA). For this, the inverse distance weighting (IDW) interpolation method with manual adjustment was utilized, based on the findings of González-Álvarez et al. [14] for the Colombian Caribbean region (inputs: Zvalue = 2; cell size = 0.021; search radius variable, and number of points = 12). Table 3. Types of isohyetal maps.

Condition Description Stationary
Isohyetal maps generated from Pmax-24h values under stationary conditions. Non-stationary Isohyetal maps generated from Pmax-24h values under non-stationary conditions.

Mixed
Isohyetal maps generated from the Pmax-24h value corresponding to the best fit according to the Akaike information criterion (AIC) test.

Evaluation of the Different Isohyetal Maps
The performance of the isohyetal maps was tested in nine watersheds (three watersheds per homogeneous region) by estimating their corresponding areal precipitation under SC, NSC, and mixed conditions for return periods of five, 10, 20, 50, and 100 years. The selected watersheds had various sizes and were located at different distances from the nearest rain gauge.
REr measures the error percentage between the true and simulated values; its optimal value is zero. RSR also evaluates the error, with an optimal value of zero; however, it does so in a standardized manner by dividing the root-mean-square error (RMSE) of the true and simulated values by the standard deviation. PBIAS estimates the bias as a percentage, with an optimal value of zero; negative values indicate overestimation, while positive values indicate underestimation. NSE is an indicator of the predictive power of a model (range of values from −∞ to one); it measures how the simulated values resemble true values (dispersion around the 1:1 line). The optimal value of NSE is one (perfect fit). Negative values indicate that it is better to use the average of true values than simulated values. Values of zero (or close to zero) indicate that either the average of true values or the simulated values could be used.
For this study, areal Pmax-24h values from the mixed isohyetals were assumed to be the true value, given that these (the isohyetals) were derived from the point Pmax-24h data of the distribution functions that performed best according to the AIC test. In Equations (21)-(24), Ptrue represents the true areal Pmax-24h, Psim corresponds to the areal Pmax-24h estimated from both the stationary and the non-stationary isohyetals, n is the number of watersheds analyzed, and i is the watershed analyzed.

Trend Detection
The results of the MK and SR tests (Table 4) showed that only the Pmax-24h time series of the Puerto Giraldo rain gauge (gray cell in Table 4) had a significant trend at a 5% level of confidence (gray cell in Table 4). The Theil-Sen slope value of −0.89 corroborated the results obtained from the MK and SR tests.
0.13 C = constant or no significant trend; DC = significant decreasing trend; IC = significant increasing trend.
Despite the fact that Puerto Giraldo was the only rain gauge with a significant trend, there were also other rain gauges with either increasing or decreasing trends. Out of the 19 rain gauges, 10 showed increasing trends, five of which had ZMK values greater than one. Likewise, eight rain gauges had decreasing trends, three of them with values below one. The trends of these time series, although currently considered as not significant (values less than |1.96|), should be evaluated in the coming years to determine any increment of the estimated values. With respect to the Theil-Sen slope results, San José and Hda El Rabón had values less than one. These two rain gauges are both located in the southern part of the department, where IDEAM [3] predicted a rainfall decrease. These findings help with (a) a better understanding of the rainfall regime (both annual and daily maximum), and (b) confirming the hypothesis raised by González-Álvarez et al. [37] as to the existence of Pmax-24h trends within the Colombian Caribbean coast. Table 5 presents the results of the dimensionless L-moments τ2, τ3, and τ4 for each of the rain gauges analyzed, which were used for identifying the homogeneous regions. Subsequently and with the purpose of defining the best homogeneous region, the K-means method was performed using clusters with different set-ups of rain gauge groups. Groups of three, four, and five clusters were defined with respect to the homogeneity presented in the variables shown in Table 5. Finally, the best configuration was selected using (a) the Hosking et al. [36] methodology, and (b) a geographical comparison among the homogeneous group distributions.

Identification and Delimitation of Homogeneous Regions
Figures 3a-c depict the rain gauge spatial distribution grouped into three, four, and five clusters, respectively. In Figures 3b, c, the green ovals show how a rain gauge belonging to another group (yellow diamond) is within a different cluster (blue circles). This indicates that grouping rain gauges into either a four-or five-rain-gauge cluster introduces geographically inconsistent distributions. In fact, Hosking et al. [37] affirmed that these types of plots (those that relate τ2, τ3, and τ4) can sometimes contain overlapping groups, which makes it difficult to select the number of clusters that adequately represent rain gauges with similar characteristics. The three-rain-gauge cluster group (Figure 3a) did not show that problem. Also, the adequacy of the three-rain-gauge cluster was further verified by analyzing the behavior of the geographic location variables (latitude and longitude) with respect to the coefficient asymmetry (τ3). Figure 4 evidences three well-defined clusters, which represent those rain gauges with similar characteristics.  Once the clusters were defined, three homogeneous regions were delineated through the IDW interpolation method (via ArcGIS) ( Figure 5), namely, east (Cluster 1), central (Cluster 2), and west (Cluster 3). Table 6 summarizes the rain gauges that make up each of the homogeneous regions.   Table 7 and Figure 6 show the rain gauges (in each of the regions) where CDFs Gumbel and GEV better represented the time series. Overall, the Gumbel distribution was the best fit in 63.2% of the 19 rain gauges analyzed, while GEV was the best fit in only 36.8%. These results represent a shift from the findings by González-Álvarez et al. [45], where GEV was best in 53.8% of the 13 rain gauges assessed. This is due to the fact that the Gumbel distribution was the best CDF among the new additional six rain gauges analyzed in this study. Based on the findings and despite the fact that the Gumbel distribution was best in the majority of the cases, there was not a unique CDF that better represented all the time series within a particular region. Pmax-24h values under stationary conditions for each of the 19 rain gauges are presented in Table A1 (Appendix A).

Non-Stationary Frequency Analysis
Pmax-24h estimates under non-stationary conditions (by means of the GEVmu distribution) for the 19 rain gauges are compiled in Table A1 (Appendix A). Differences of up to 58.9 mm were observed at the El Porvenir rain gauge between the Pmax-24h values under stationary and non-stationary conditions for a 100-year return period.

Selecting the Best Pmax-24h Value
After the estimation of the Pmax-24h under both stationary and non-stationary conditions, it was necessary to determine-via AIC test-which of the two conditions better represented the time series (Table 8). The values obtained in this section were later used for the generation of isohyetal maps called mixed (derived from the best rainfall value of the two conditions) explained in the next section. Table 8 shows the best condition (stationary or non-stationary) and CDF, as well as the AIC test values obtained. Figure 7 depicts the time series of rain gauges at Casa de Bombas, Hda El Rabón, Puerto Giraldo, and Los Campanos. For the first two rain gauges, a stationary condition frequency analysis is best, according to the AIC test, while the last two suit a non-stationary one.

Isohyetal Maps
With the Pmax-24h values of the 19 rain gauges obtained in Sections 4.3.1-4.3.3, maps of stationary (Figure 8), non-stationary (Figure 9), and mixed (the best value of the two conditions according to the AIC test, Figure 10) isohyetals were drawn for return periods of five, 10, 25, 50, and 100 years.    Figure 8 shows that isohyetal values ranged from 100-110 mm for five years, from 110-120 mm for 10 years, from 120-150 mm for 25 years, from 130-160 mm for 50 years, and from 140-180 mm for 100 years. Figure 9 shows values of 90-120 mm for a return period of five years, 100-130 mm for 10 years, 110-160 mm for 25 years, 110-180 mm for 50 years, and 110-210 mm for 100 years.
Finally, in Figure 10, it can be observed that the isohyetals ranged from 100-120 mm for five years, from 110-130 mm for 10 years, from 120-150 mm for 25 years, from 130-180 mm for 50 years, and from 140-190 mm for 100 years. The highest values were observed in the south zone and the lowest values were observed in the north.

Isohyetal Map Assessment
In order to assess how the use of stationary and/or non-stationary isohyetal maps could affect the calculation of the areal Pmax-24h (Pareal) in a given watershed (W), nine watersheds were selected, three in each of the three homogeneous regions (light-green areas in Figure 11). W1, W2, and W3 are located in the north, central, and southern areas of the west homogeneous region, respectively. W4, W5, and W6 are within the central region, located in the northern, central, and southern areas, respectively. Finally, W7 (north), W8 (center), and W9 (south) are located within the east region. Table  9 summarizes the watershed area, the nearest rain gauge, and its distance to each of the watersheds.  Pareal values for all watersheds (return periods of five, 10, 25, 50, and 100 years) for stationary and non-stationary conditions, as well as mixed ones, are shown in Table 10. Likewise, stationary and non-stationary Pareal values were compared, through their differences (mixed minus the SC and NSC values), with Pareal values of mixed isohyetal maps (Table 10).  The frequency analysis under stationary conditions is most commonly used by hydrologists to estimate the design rainfall for hydraulic structures for stormwater management. Nonetheless, the differences observed in Table 10 show that the stationary frequency analysis underestimated the values of areal Pmax-24h for the study area. This behavior occurred in 60% of the cases evaluated (27 out of 45) throughout the department of Atlántico. This implies that, if a designer decides to use a stationary design rainfall, the subsequent estimation of the design flow for a given hydraulic structure could end up as an underestimated value. It can also be observed that these underestimations reach their most critical values within the central region, where the highest areal Pmax-24h differences range from 11.3 to 24.2 mm (both in W5) for return periods of 25 (which are the most used in the design of drainage hydraulic structures). The west region exhibited only two cases where the difference in Pareal had values greater than 5 mm (7.7 mm in W1 and 8.9 mm in W2). In the eastern region, underestimations occurred in 60% of cases (six cases out of 10) with values ranging from 6.8 to 11.2 mm (W7 and W8). Like the central region, these values came from return periods of 25, 50, and 100 years. These results also show that, for stationary conditions, the probability of underestimating the values of Pareal is higher in these three return periods, with values of 66%, 88%, and 66%, respectively. Rainfall differences with values less than 5 mm occurred in 57.8% of the cases (26 cases distributed as follows: nine for the west region, nine for the east region, and eight for the central region, where two values of 5.1 mm were given that could be considered within the rank).
For the non-stationary scenario, the tendency to underestimate Pareal occurred in 46.7% of the cases (21 out of 45): seven cases for the west region, nine for the central region, and five for the east region. It was also noted that, for this scenario, more Pareal values with differences less than or equal to 5 mm are estimated, which was observed in 66.7% of the cases (11 for the west region, nine for the central region, and 10 for the eastern region). On the other hand, when each of the regions was analyzed individually, the western region showed a slight tendency to overestimate the Pareal values (eight out of 15, negative values in red in Table 10). This behavior was more noticeable for the 100year return period, where the three watersheds evaluated presented values ranging from 2.2 mm (W2) to 13.3 mm (W3). In the central region, underestimated values of Pareal (60% of the cases or nine out of 15) were observed, particularly in the southern part of this region (W6), with values ranging from 11.1 to 26.5 mm for return periods 10, 25, 50, and 100 years. For the eastern region, there was also a tendency to overestimate the Pareal (66.7% of the cases, 10 of 15), with values ranging from 2.1 mm (W9) to 19.5 mm (W8).
Regarding the real values of Pareal (mixed isohyetals), the southern area of the department of Atlántico exhibited the lowest values, specifically in W3 (located in the south of the west region) and W6 (located in the south of the eastern region). This is clearly evident, for example, when comparing the Pareal values for the 100-year return period between W6 and W7 (located northeast). W6, despite having an area three times larger than W7, has a lower Pareal value (153.3 mm versus the 155.0 mm for W7). These results coincide with the findings of IDEAM [3], who determined that municipalities located in the southeast of the department will likely be the most affected by rainfall decrease.
The isohyetal map (stationary and non-stationary) performance assessment within each of the watersheds was carried out through the relative error percentage (REr). Additionally, the performance of each of the regions was assessed by RSR, PBIAS, and NSE. The results of this statistical analysis are presented in Table 11. For the stationary isohyetal maps, REr values ranged from 0.00% (W2) to 9.52% (W3) for the fiveyear return period, 0.77% (W6) to 6.17% (W4) for the 10-yearperiod, 0.13% (W3) to 8.84% (W5) for the 25-year period, 0.17% (W6) to 9.59% (W5) for the 50-year period, and 0.23% (W9) to 15.64% (W5) for the 100-year period. The highest REr values were observed in two of the return periods (50 and 100 years) most commonly used in the design of hydraulic structures for runoff management (50 and 100 years). No relationship was observed between the watershed area and REr.
In general, for stationary conditions, there were 14 cases where REr was greater than or equal to 5.00% (gray cells). The maximum value was 15.64%, with only one case where REr was greater than or equal to 10%. For the non-stationary conditions, 10 cases were observed where the REr was greater than or equal to 5.00%. The maximum value was 20.85%, with five cases where REr was greater than or equal to 10%. Furthermore, when the stationary and non-stationary conditions were compared one-to-one, it was observed that, in 57.8% of the cases (26 out of 45), the REr values for the nonstationary conditions were lower than their stationary counterparts. These results suggest that (a) the error might be more frequent when using the stationary condition isohyetal maps, and (b) additional attention should be paid during the design of hydraulic structures under stationary frequency analysis, especially as it was also found that this scenario tends to underestimate the Pareal (Table 10).
With respect to the overall performance of all regions, the stationary conditions resulted in lower values of RSR (10 in total) than those for the non-stationary conditions (five in total). At first glance, this may indicate less error under stationary conditions (which contradicts the results previously obtained when the REr was analyzed for each watershed). Nonetheless, a closer look at Table 10 revealed that, despite the fact that each condition had five REr values greater than or equal to 5.00%, the non-stationary condition had REr values of up to 20.85%, which contributed to having an overall larger RSR value. Such large REr values were due to the fact that W6 happened to have a rain gauge (Hda El Rabón) with a time series better suited to a stationary frequency analysis ( Table 8). As for the individual performance of each region, the west region showed 90% (nine out of 10) of the RSR values below one, followed by the east region with three values. The central region had more stationary condition values (in four out of the five return periods) that outperformed the non-stationary ones.
Regarding the Pareal tendency to under-or overestimate, PBIAS values indicate that isohyetal maps under stationary conditions tend to underestimate (black positive values in Table 11) in the majority of the cases (66.7% or 10 out of 15), which is more evident in the central region. These results corroborate what was previously found in Table 10. The underestimated results also observed in the central region for the non-stationary conditions (which are opposite to the results of both Table 10 and REr in Table 11) were mainly caused by the large Pareal differences found in W6. The nonstationary isohyetal maps tend to underestimate (60% in all regions, or nine out of 15). In the central region, the underestimation occurred in all return periods for both stationary and non-stationary conditions, especially for 50-and 100-year periods, which are two of the return periods most used in the design of hydraulic structures. In the east region, a tendency to overestimate (red values in Table  11) was detected for non-stationary conditions in four out of five return periods. A different behavior was observed for the stationary conditions within the same region (east) where underestimation prevailed. Overall, the west region showed less bias when compared with the other two regions, with values ranging from 0.61% (underestimation for the 25-year return period for stationary conditions) to −4.21% (overestimation for the five-year return period for non-stationary conditions). Central and east regions showed values oscillating from 0.72% (underestimation for the five-year return period for non-stationary conditions) to 5.80% (underestimation for the 100-year return period for stationary conditions).
With regard to the prediction power of the isohyetal maps under stationary and non-stationary conditions, better NSE results were obtained within the west region in the majority of cases. Both conditions (stationary and non-stationary) had three return periods with NSE values above 0.5 (blue cells in Table 11). Among the return periods most used for the design of hydraulic structures for stormwater management (25,50, and 100 years), 25-and 50-year periods showed values close to one (indicator of a good performance), with values above 0.70 for both conditions. For the 100-year period, a value of 0.39 was observed for stationary conditions, denoting good performance as well. However, for non-stationary conditions, a value of 0.09 denotes both that the simulated value is far from the 1:1 line and that the average value of either the simulated or true value better represents the areal rainfall value. For the central and east regions, negative values prevailed in most of the cases for either stationary or non-stationary conditions (only one value was above 0.5). Within the central region for stationary conditions, all return periods had negative values, while, for the non-stationary conditions, this behavior was seen in 80% of the cases. In the east region, four out of five return periods showed negative values for stationary conditions, and three out of five return periods showed negative values for non-stationary conditions. These results indicate that the average of the true value (mixed isohyetal maps) is a better predictor for these two regions.
In general, lower values of REr, RSR, and PBIAS were observed within the west region, especially for stationary conditions, which suggests that a stationary frequency analysis might be used in watersheds within this region. This was also confirmed by the NSE results obtained in four out of five of the return periods. For the central and east regions, the use of a stationary frequency analysis (typically and widely used in hydrology), according to the results obtained, might introduce errors in the calculation of Pareal, which could affect, for instance, the magnitude of the estimated runoff for water balances (for agriculture, livestock, and energy water demand, among other uses), hydraulic structures for stormwater management, flash flood guidance, and flood risk assessment.

Conclusions
With respect to the Pmax-24h behavior, three regions were determined, namely, east, central, and west. The regionalization will be of great help for the Pmax-24h analysis in ungauged areas given the fact that the department of Atlántico is, among the remaining six departments of the Colombian Caribbean region, the one with the lowest rain gauge density (only 19 with statistically representative time series).
Increasing and decreasing trends were identified among the 19 Pmax-24h time series analyzed within the department of Atlántico. Furthermore, only one rain gauge showed a significant decreasing trend with values of ZSR, ZMK, and βTS of 1.36, −2.06, and −0.89, respectively. However, other rain gauges also showed increasing and decreasing trends, for which, despite not being significant, five of them showed ZMK greater than one and three had values less than negative one. This suggests the need for future trend analysis in the coming five-year periods to determine any further trend increase/decrease. Overall, the southern area of the central and west regions showed the most noticeable decreasing trend. These results are in agreement with IDEAM [3] findings.
As to which frequency analysis-stationary or non-stationary-better represented the 19 time series analyzed, the AIC test revealed that 79% of them suited a stationary one. In terms of the performance of the isohyetal maps under stationary and non-stationary conditions when compared with the mixed (stationary along with non-stationary), REr values indicate that, while the error under stationary conditions can be observed more frequently, that under non-stationary conditions could be more significant in terms of magnitude, especially in the southern central region. This was also confirmed by the RSR and PBIAS results, where the non-stationary condition, despite having less cases with REr greater than 10% among the nine watersheds evaluated, the results showed how the magnitude of the error impacts the overall results within a given region. In sum, the west region had fewer cases (watersheds) with REr values above 10% under both stationary and non-stationary conditions. Likewise, RSR, PBIAS, and NSE also indicated that either a stationary or a non-stationary frequency analysis might be performed in the estimation of the areal Pmax-24h, which represents a contribution to the hydrological analysis given that, according to the results of this study, a stationary frequency analysis (the most commonly used) might be safely performed within the west region. On the other hand, the other two regions presented a tendency for underestimation, especially under stationary conditions, which indicates, for example, that hydraulic structures for stormwater management should be designed with precaution.
The findings of this study shed some light on the need of a better understanding of both the regional hydrological behavior and the impact of climate change on future water-related projects.