Isohyetal Maps of Daily Maximum Rainfall for Different Return Periods for the Colombian Caribbean Region

: In Colombia, daily maximum multiannual series are one of the main inputs for design streamﬂow calculation, which requires performing a rainfall frequency analysis that involves several prior steps: (a) requesting the datasets, (b) waiting for the information, (c) reviewing the datasets received for missing or data different from the requested variable, and (d) requesting the information once again if it is not correct. To tackle these setbacks, 318 rain gauges located in the Colombian Caribbean region were used to ﬁrst evaluate whether or not the Gumbel distribution was indeed the most suitable by performing frequency analyses using three different distributions (Gumbel, Generalized Extreme Value (GEV), and Log-Pearson 3 (LP3)); secondly, to generate daily maximum isohyetal maps for return periods of 2, 5, 10, 20, 25, 50, and 100 years; and, lastly, to evaluate which interpolation method (IDW, spline, and ordinary kriging) works best in areas with a varying density of data points. GEV was most suitable in 47.2% of the rain gauges, while Gumbel, in spite of being widely used in Colombia, was only suitable in 34.3% of the cases. Regarding the interpolation method, better isohyetals were obtained with the IDW method. In general, the areal maximum daily rainfall estimated showed good agreement when compared to the true values. for the OK the values of 24h-max -IM and


Introduction
Designing hydraulic structures for stormwater management encompasses several tasks, among which are: (a) watershed morphometric analysis, (b) estimation of the time of concentration, (c) calculation of the design rainfall via frequency analysis (typically, under stationary conditions), (d) design flow computation, (e) sizing the hydraulic structure per se, and (f) hydraulic modeling to evaluate the structure's hydraulic performance under various return periods. The design flow may be estimated via either a rainfall-runoff model or regression equations (in ungauged watersheds), or stationary frequency analysis of streamflow data, if available.
Unlike streamflow data, rainfall observations are the most abundant hydrometeorological variable available. Rainfall observations, as a result, are most commonly used when estimating the flow and mitigation, environmental impact assessments, land development, stream restoration, and the design of hydraulic structures for stormwater management. Unfortunately, there is no such map in Colombia. The lack of readily-available and processed information might become an obstacle in the understanding and development of these projects. Furthermore, in the absence of data, detailed isohyetal maps such as the ones developed in this research can be utilized for both regional and global analyses of extreme rainfall event proxies, hydrological regime changes, climate studies, water balance estimation under various scenarios, and the development of water management strategies.
This study intends to contribute to these tasks by: (a) first carrying out rainfall stationary frequency analyses by means of three different Cumulative Distribution Functions (CDF) for extreme values analysis (Gumbel, GEV, and LP3) so as to confirm whether or not the Gumbel distribution is indeed the most suitable for the rainfall data of this study, as it is the most commonly used in Colombia despite the fact that the other two have been also assessed in various studies [18][19][20][21][22][23]; (b) secondly, evaluating which spatial interpolation method is most suitable given the highly-dense spatial distribution of the rain gauges; and (c) lastly, drawing isohyetal maps of P 24h-max -via Geographical Information System (GIS)-for return periods of 2, 5, 10, 20, 25, 50, and 100 years for the Colombian Caribbean region as the first step in a more ambitious project of elaborating maps for the remaining regions (Pacific, Andean, Orinoco, and Amazon). For that, 318 pluviometric stations (rain gauges) operated by IDEAM-with at least 30 years of data-were used. The rain gauges are distributed as follows: 313 throughout the seven departments that compose the Colombian Caribbean region (Guajira, Magdalena, Cesar, Atlántico, Bolívar, Sucre, and Córdoba) and five from neighboring departments (Antioquia, Santander, and Norte de Santander).

Study Area and Data
The Caribbean region of Colombia is comprised of seven political and administrative territorial units (called departments), namely Guajira, Magdalena, Cesar, Atlántico, Bolívar, Sucre, and Córdoba, which sum up a total area of 132,244 km 2 (accounting for approximately 0.012% of the country's total surface area). It has an average annual precipitation that ranges from 0-1500 mm (the region's northern portion) and from 1500 mm up to 5000 mm (in some areas of the southern portion of the departments of Bolívar and Córdoba) [4]. The region is mostly composed of plains, with the exception of the mountain ranges of San Jacinto (located within the departments of Sucre and Bolívar) and Santa Marta (department of Magdalena). The rainiest departments are Bolívar and Sucre, while the driest is Guajira. The average annual temperature is about 30 • C with some areas where it may increase or decrease depending on the climatological season and altitude. The rainfall climate regime of the Colombian Caribbean region has three seasons [24]: dry (December-March), transition-also known as Veranillo de San Juan-(June-July), and rainy (April-June and August-November). The duration of these seasons is affected by the El Niño Southern Oscillation (ENSO) phenomenon [25][26][27][28][29].
The data used in this study correspond to multiannual series of P 24h-max that include a total of 12,828 observations till the year 2015 from 318 pluviometric stations (the oldest station in this study started to operate in 1931). Years 2016 and 2017 were not included because the data from some of the stations have not been officially released by IDEAM. The pluviometric stations selected had to comply with the following criteria: (a) being still operative and (b) having at least 30 years of observations [30]. Additionally, P 24h-max values were eliminated if they came from a year not having a minimum of 150 days of data [31][32][33] and/or missing any of the months of the rainy season. The resulting rainfall data are summarized in Table 1. Figure 1 shows the departments and the geographical location of the rain gauges. Out of the 318 rain gauges, 183 have 30-40 years of data, 113 have 41-50 years of data, and 22 have more than 50 years of rainfall observations. It may also be observed that the southern areas of the departments of Bolívar and Córdoba have no rain gauges that fit the selection criteria. To fill these gaps, it is necessary to improve the isohyetal alignment, so five additional rain gauges were used. The rain gauges are located within the three neighboring departments of Antioquia (three rain gauges), Santander (one rain gauge), and Norte de Santander (one rain gauge). Despite the fact that these three departments do not belong to the Caribbean region, the rain gauges selected are located close to the southern portion of the Caribbean region. Table 2 summarizes the names of the rain gauges analyzed in this study. the rain gauges selected are located close to the southern portion of the Caribbean region. Table 2 summarizes the names of the rain gauges analyzed in this study.    normal, Pearson 3, and log-Pearson 3, among others. When a frequency analysis assumes that the CDF used does not change over time, it is called stationary [34,35]. The CDF represents the probability of no-exceedance (probability that the value analyzed is less than or equal to the rest of the values in the dataset). With respect to the variable (or variate) analyzed, it must be random (uncertainty in its prediction) and independent (its occurrence is not affected by other variables) [36]. Extensive literature exists on this topic [1,2,[37][38][39][40]. Gumbel is a two-parameter unbounded distribution (the shape parameter, k, is zero for this distribution), which uses a double exponential function to estimate the probability of exceedance [F(Z)] [41,42]. In Equations (1), (2), and (3), z is the random variable (rainfall, streamflow, wind, etc.), β is the mode (location parameter), α is the dispersion (scale parameter), and ̅ and are the mean

Stationary Frequency Analysis
In hydrology, there are several CDFs for the analysis of extreme values, namely: gamma, Extreme Value (EV) or Fisher-Tippett (Types 1, 2, and 3), Generalized Extreme Value (GEV), log-normal, Pearson 3, and log-Pearson 3, among others. When a frequency analysis assumes that the CDF used does not change over time, it is called stationary [34,35]. The CDF represents the probability of no-exceedance (probability that the value analyzed is less than or equal to the rest of the values in the dataset). With respect to the variable (or variate) analyzed, it must be random (uncertainty in its prediction) and independent (its occurrence is not affected by other variables) [36]. Extensive literature exists on this topic [1,2,37-40].

Gumbel Distribution (Extreme Value Type 1 or Fisher-Tippett Type 1)
Gumbel is a two-parameter unbounded distribution (the shape parameter, k, is zero for this distribution), which uses a double exponential function to estimate the probability of exceedance [F(Z)] [41,42]. In Equations (1), (2), and (3), z is the random variable (rainfall, streamflow, wind, etc.), β is the mode (location parameter), α is the dispersion (scale parameter), and z and σ z are the mean and the variance of the random variable, respectively. It is one of the most-used CDFs by hydrologist and other water-related professionals in Colombia. It is also used in New Zealand for rainfall frequency analysis.

Generalized Extreme Value Distribution
The GEV is a three-parameter distribution that compiles all three types of the EV distributions into one formula (Equation (4)). In Equation (4), F(Z) is the CDF, z is the random variable, k is the shape (or shift) parameter, β is the mode (location parameter), and α is the dispersion (scale parameter, always assumed to be greater than zero). The GEV distribution may adopt one of the three EV distributions depending on the value of k [1,2,40]: (a) when k equals zero, EV is Type 1 (Gumbel) [41,42]; (b) when k is less than zero, EV is Type 2 (Fréchet) [43,44]; and (c) when k greater than zero, EV is Type 3 (Weibull) [45].
This distribution is one of the most widely used. Countries like the United States (especially in the eastern portion), the United Kingdom, Australia, Bangladesh, Canada, South Africa, and New Zealand have adopted it in some areas [46][47][48][49][50][51].

Goodness-of-Fit Test
The goodness-of-fit of the three CDFs used in the frequency analysis was evaluated via the chi-squared test (X 2 ) with a significance level of 0.05.
In Equation (9), O i and P i are the observed and predicted distributions, respectively. The theory behind this test may be found throughout the literature on statistical hydrology [1,2,38,40]. In general, when testing several CDFs for frequency analysis, the lower the value of X 2 , the better the CDF fits the dataset.

Estimation of P 24h-max for Different Return Periods
All datasets for each of the 318 rain gauges were subjected to a stationary frequency analysis by using the above-described distributions. The return period (Tr) or recurrence interval is a concept commonly misinterpreted, as some describe it as the time (in years) it takes an event (rainfall, streamflow, etc.) to occur again. Instead, Tr must be understood as the occurrence of a given rainfall event, in any particular year, that may be equal to or exceeded by some percentage. The return period and the probability of exceedance (Pe) are inversely proportional (Tr = 1/Pe). The Tr is used in many fields: hydrology, hydraulic structures design, protection of water bodies receiving wastewater discharges (low flow indices' estimation), and ecology, just to mention some. The parameters of the Gumbel and GEV distributions were estimated via the maximum likelihood method [1,2,40,[58][59][60], and the Sundry Averages Method (SAM) method was used for LP3 [1,2,40,[52][53][54].
The values of P 24h-max selected (to be later used in the drawing of the isohyetals) for the return periods of 2, 5, 10, 20, 25, 50, and 100 years were those that came from the CDF showing the best result in the goodness-of-fit test. Additionally, a test for outliers was performed ([1], pp. 403-405). Two rain gauges were found to have outliers, namely San Cayetano (department of Bolívar) and Irán (department of Magdalena). The outliers were eliminated, and a new frequency analysis was made for these two rain gauges.

Drawing of Isohyetals for Different Return Periods
A spreadsheet with the geographical location (latitude and longitude) of the 318 rain gauges along with their corresponding estimated P 24h-max values for the different return periods was exported into ArcGIS (Version 10.4.1, ESRI Inc., Redlands, CA, USA) to create a layer for further processing. Given the large density of rain gauges in some areas, three interpolation methods were used to generate the isohyetals: Inverse Distance Weighting (IDW), Ordinary Kriging (OK), and Spline (SP), in order to evaluate their performance prior to selecting one of the methods [61][62][63]. Default ArcGIS inputs were used in all methods since a sensitivity analysis done for each of the methods by changing the various inputs showed no major changes or improvements in the areas with noticeable alignment discrepancies. The ArcGIS inputs for the three methods are shown in Table 3.

Assessing the Interpolation Methods
The accuracy of each interpolation method in predicting a P 24h-max value for a given return period was evaluated by means of the Root Mean Squared Error (RMSE) (Equation (10)), the Relative Error (REr) (Equation (11)), the Bias or Mean Deviation (MDv) (Equation (12)), and the Nash-Sutcliffe Coefficient (NSC) (Equation (13)). In Equations (10)- (13), P sim is the simulated areal precipitation from each of the Interpolation Methods (P 24h-max -IM), P sim is the average simulated areal precipitation, and P true is P 24h-max -RG, which has been considered the true value due to the proximity to the rain gauges. RMSE and REr measure the accuracy, whereas MDv tests the bias. The lower the values of RMSE, REr, and MDv, the better the interpolation method. The NSC, in particular, is widely used in hydrology to assess the prediction power of a model, and it ranges from −∞-one, where negative values indicates that it is better to use the mean of the measured data (true value) than the predicted/simulated value; a value of zero (or close to zero) indicates that either the mean of the measured values or the predicted/simulated value can be used indistinctively; and a value equal to one indicates good agreement between the predicted/simulated and the measured data (true value) [9,64,65].

Multiannual Time Series of P 24h-max Values
Scatter plots of the multiannual series of P 24h-max of each rain gauge revealed that there might be regionalization of the daily maximum rainfall trends within the departments that need to be further explored and analyzed, as it was observed that rainfall observations showed a noticeable increasing or decreasing trend line over time, which may indicate (a) a change in the rainfall pattern due to, among others, anthropogenic factors and (b) that a non-stationary frequency analysis is more suitable for the rainfall data of those rain gauges at a local level [66,67]. Figure 3 shows the scatter plot and the trend lines of the rain gauges Puerto Giraldo and Los Campanos located in the department of Atlántico. values indicates that it is better to use the mean of the measured data (true value) than the predicted/simulated value; a value of zero (or close to zero) indicates that either the mean of the measured values or the predicted/simulated value can be used indistinctively; and a value equal to one indicates good agreement between the predicted/simulated and the measured data (true value) [9,64,65].

Multiannual Time Series of P24h-max Values
Scatter plots of the multiannual series of P24h-max of each rain gauge revealed that there might be regionalization of the daily maximum rainfall trends within the departments that need to be further explored and analyzed, as it was observed that rainfall observations showed a noticeable increasing or decreasing trend line over time, which may indicate (a) a change in the rainfall pattern due to, among others, anthropogenic factors and (b) that a non-stationary frequency analysis is more suitable for the rainfall data of those rain gauges at a local level [66,67]. Figure 3 shows the scatter plot and the trend lines of the rain gauges Puerto Giraldo and Los Campanos located in the department of Atlántico.

CDFs and Frequency Analysis
As previously mentioned, in Colombia, it is a common practice among hydrologist to use the Gumbel distribution solely for rainfall frequency analysis. However, results presented in Table 4 show that GEV is the CDF that best fit the majority of the rain gauges analyzed in this study with 47.2%, followed by Gumbel with 34.3% and LP3 with 18.5%.

CDFs and Frequency Analysis
As previously mentioned, in Colombia, it is a common practice among hydrologist to use the Gumbel distribution solely for rainfall frequency analysis. However, results presented in Table 4 show that GEV is the CDF that best fit the majority of the rain gauges analyzed in this study with 47.2%, followed by Gumbel with 34.3% and LP3 with 18.5%.
Furthermore, based on the value obtained for the shape parameter (k), the 150 rain gauges where GEV was found to be the best fit were further analyzed to determine whether the GEV corresponded to a Type 2 (Fréchet) or Type 3 (Weibull) Extreme Value (EV) distribution. Table 5 summarizes the EV Type 2 and 3 totals for each department. Taking for granted that all rainfall time series fit the Gumbel distribution could lead to either under-or over-estimation of the design rainfall. This is clearly shown in Table 6, where, according to the chi-squared test, the Gumbel distribution was not the best fit in any of the cases. The situation became more critical at higher return periods such as 50 years and 100 years (two of the most-used return periods in hydrological analysis for stormwater management and in consideration for flood studies). For the 100-year return period, the differences observed between GEV and Gumbel distributions at rain gauges Santa Ana and Palo Alto were 54.1 mm and 15.2 mm, respectively. Not selecting the most appropriate distribution can negatively impact, for instance, the design of hydraulic structures for stormwater management and/or the delineation of flood-prone areas. In general, the GEV and Gumbel distributions were demonstrated to be the most suitable for rainfall frequency analysis in the Caribbean region of Colombia. The LP3 distribution showed both poor performance (in some cases, the datasets did not even fit) in most of the rain gauges analyzed and a tendency to underestimate when compared to the values obtained by either GEV or Gumbel distributions, irrespective of being (or not) the best fit of the three.

Interpolation Methods Assessment
After visually inspecting the isohyetal alignments, it was observed that the IDW method had less inconsistencies, among all. Nonetheless, in some areas, it was necessary to adjust manually the isohyetals' alignment generated by the IDW method (the manually-adjusted IDW method will be referred to as IDW adjusted). Figure 4 depicts the differences among the interpolation methods. Isohyetals by the OK method ( Figure 4c) evidenced: (a) the presence of small oval-shaped isohyetals with the same rainfall value of the larger oval they were within and (b) large areas between neighboring rain gauges with no isohyetals ("dead or no-variation zones"), which could affect the estimation of the areal precipitation for a given watershed (the problem was most evident as the return period decreased). The spline method (Figure 4d) generated, in some cases, isohyetals with negative values. Isohyetals drawn by the IDW method (Figure 4b), despite the fact that a few minor adjustments were manually made in some areas, did not show the setbacks of the other two methods. The IDW method, in spite of its simplicity compared to other interpolation methods like the OK method, is recommended when the data are irregularly distributed like the rain gauges used in this study [68,69].
In addition to the visual inspection, the accuracy of all generated maps (by each of the interpolation methods) to estimate P 24h-max values for a given return period was evaluated in eight watersheds with various area sizes and located at different distances from neighboring rain gauges. Five watersheds are located in the northern part of the department of Bolívar. Watersheds 1, 2, and 3 (W-1, W-2, and W-3) are close to the rain gauges Bayunca, Cañaveral, and Escuela Naval-CIOH, respectively. The rain gauges Bayunca and Cañaveral are located within Watersheds 4 and 5 (W-4 and W-5). Furthermore, three watersheds (W-6, W-7, and W-8) located in the vicinity of the rain gauge Loma Grande (department of Atlántico) were selected to validate the performance of the isohyetal methods in areas with nearby rain gauges not included in this study (Loma Grande had less than thirty years of observations). Table 7 summarizes the watersheds' area and distance from the centroid to the nearest rain gauge. Figure 5 depicts the location and the isohyetals at each of the watersheds for a 100-year return period.
Watersheds areal P 24h-max for return periods of 2, 5, 10, 20, 25, 50, and 100 years were estimated from each of the isohyetal maps (P 24h-max -IM) and compared to the resulting P 24h-max values calculated via frequency analysis of the time series of each of the nearby rain gauges (P 24h-max -RG).
The estimated values of P 24h-max -RG and P 24h-max -IM are presented in Table 8. P 24h-max -RG is the value obtained via frequency analysis of the precipitation data for all return periods at each of the rain gauges used for validation. P 24h-max -IM is the estimated areal precipitation of each watershed by means of the isohyetal method. The observed two-year P 24h-max -RG values ranged from 75.4 mm (Loma Grande) to 100.5 mm (Bayunca). Among all rain gauges, Loma Grande had the lowest P 24h-max -RG for almost all return periods, except for 100 years. The highest P 24h-max -RG values calculated were for the rain gauges of Bayunca (for 2 years, 5 years, and 10 years) and Cañaveral (for 20 years, 25 years, 50 years, and 100 years). The minimum two-year P 24h-max -IM values estimated were 91.     The estimated values of P24h-max-RG and P24h-max-IM are presented in Table 8. P24h-max-RG is the value obtained via frequency analysis of the precipitation data for all return periods at each of the rain gauges used for validation. P24h-max-IM is the estimated areal precipitation of each watershed by means of the isohyetal method. The observed two-year P24h-max-RG values ranged from 75.4 mm (Loma Grande) to 100.5 mm (Bayunca). Among all rain gauges, Loma Grande had the lowest P24h-max-RG for almost all return periods, except for 100 years. The highest P24h-max-RG values calculated were for the rain gauges of Bayunca (for 2 years, 5 years, and 10 years) and Cañaveral (for 20 years, 25   With the values estimated in Table 8, the performance of the interpolation methods at each of the watershed was tested via REr, and the overall performance of each interpolation method in predicting P 24h-max for a given return period in all watersheds (n = 8) was tested via RMSE, MDv, and NSC (Table 9). The REr values in Table 9 show how the IDW adjusted outperformed the other two methods in almost all watersheds except in 18 cases (gray cells) out of 168 (nine of those 18 cases occurred in watersheds where the rain gauge was within them). With respect to the watersheds having a rain gauge within them (W-4, W-5, and W-6), better REr results were obtained in the IDW adjusted method, except for W-6, where the OK method predicted values closer to the ones calculated via frequency analysis. The IDW adjusted method showed REr estimates ranging from 0.8-4% for W-4, 0.6-6.2% for W-5, and 1.3-14.6% for W-6. The spline method had REr estimates that oscillated from 0.4-3.2% for W-4, 0.4-8.2% for W-5, and 3.4-17.3%. The REr values obtained through the OK method ranged from 0.8-12.5% for W-4, from 0.6-26.4% for W-5, and from 0.0-7.5% for W-6. The IDW adjusted method reported fifteen REr values above 5%, of which only two cases were greater than 10% (the highest was 14.6%). The spline method reported 36 cases where REr values were above 5%, with 17 values above 10% (the highest was 20.8%). Finally, the OK method had 27 cases where REr estimates were above 5%, with 11 values above 10% (the highest was 26.4%). When the performance of the interpolation methods was evaluated based on the distance between the watershed and the rain gauge, once again, the IDW adjusted method exhibited better results, with a few exceptions observed (gray cells in W-2, W-7, and W-8). For W-1 and W-3 (respectively, the closest and the furthest), lower REr values were obtained with the IDW adjusted method. For W-2, located at 9.2 km of its corresponding closest rain gauge, the IDW adjusted method had lower estimates of REr for all return periods, except at two years, where the OK method was best. However, the IDW adjusted method reported an error of only 3.7%. For W-8, also located at 9.2 km, for return periods of 10 years, 25 years, and 50 years, the OK and spline methods were best. Nonetheless, all REr values observed for the IDW adjusted method were lower than 5%, except for the return period of 100 years with a value of 11%. For W-7, located at 6.0 km from the rain gauge Loma Grande, the OK method reported lower REr values, which ranged from 0.0-7.5%. The IDW adjusted and spline methods reported similar REr values (less than 5%) in most of the return periods, except for the return period of 100 years, where values of 11% and 20.8% were obtained for the IDW adjusted and spline methods, respectively. As for the return periods of 50 years and 100 years, two of the most-used return periods when evaluating the performance of hydraulic structures (e.g., culverts, bridges, and open channels) and in flood mitigation and stream restoration projects, the IDW adjusted method reported better results, with only six cases (four for the spline method and two for the OK method) out of 42, where the other two methods outperformed. The IDW adjusted method reported four cases above 10% (the highest was 14%), while spline had six cases (the highest was 20.8%). Although the OK method also reported four cases with REr values greater than 10%, the maximum value was 26.4%, which is the largest of the three methods. Irrespective of the watershed size and/or distance from a nearby rain gauge, these results are not only indicative that the IDW adjusted method performed best in the majority of the watersheds where the rain gauge was within them (W-4, W-5, and W-6), but also-and most importantly-in areas distant from a rain gauge where the design rainfall estimation was most needed, since in watersheds where a rain gauge is within them, it is always recommended to calculate the design rainfall directly from the rainfall observations.
As for the overall performance of the interpolation methods, the IDW adjusted method had the best results (RMSE, MDv, and NSC) in the majority of the cases. The OK method showed lower bias (MDv) values for two occasions (−1.07 mm and 0.75 mm for return periods of two and five years; negative values indicate that average simulated value was higher than the true value). This notwithstanding, the IDW adjusted method estimates were low as well (2.05 mm and 1.53 mm for return periods of two and five years). A closer look at the isohyetal alignment of both methods within the assessed area revealed that the OK method alignment for isohyetals of 70 mm, 80 mm, 90 mm, and 100 mm did not pass through several rain gauges they were supposed to have (dead or no-variation zones). This would have resulted in high error values for the OK method had other watersheds been selected. In general, the results indicate that the values from the IDW adjusted method showed less bias (2.5 mm-9.5 mm) than the other two methods (5.74-16.52 mm and −1.07-16.78 mm for the spline and OK methods, respectively), where the bias was more noticeable as the return period increased. This can lead to serious implications, especially at return periods of 50 years and 100 years, which are two of the most used in the design of hydraulic structures for both stormwater management and flood mitigation projects. As for the prediction power, the IDW adjusted method outperformed the other two methods with NSC values greater than or equal to 0.39, indicating good agreement between the true and simulated variables. For the spline method, the NSC was above zero only in return periods of two and five years. The remaining return periods showed negative values, implying that the average value was a better predictor. Similar results were observed for the OK method, with values closer to zero as the return period increased, which signifies that either the average or the true values were better predictors. The NSC values for the spline and OK methods are consistent with the bias results: poor performance in predicting P 24h-max when compared with the IDW adjusted method.
In sum, the results above confirm that the performance of an interpolation method should always be assessed prior to its selection [68][69][70][71][72][73][74][75]. This is of great importance, especially, because the OK method has become one of the most widely-preferred interpolation methods to the point that some do not even question its adequacy due to it typically showing good results.

Isohyetal Maps for Different Return Periods
The seven isohyetals maps of daily maximum precipitation for return periods of 2, 5, 10, 20, 25, 50, and 100 years drawn by means of the IDW adjusted method for the Colombian Caribbean region are shown in Figures 6-12. The statistical analysis results demonstrated that areal P 24h-max can be estimated by means of the isohyetal maps proposed in this study. However, the final decision of using the maps shall be at the user's discretion based on his/her experience and knowledge of the area where the areal P 24h-max is intended to be estimated. It is then recommended that a frequency analysis of the multiannual series be performed, and the results shall be compared to those estimated through the isohyetal maps in order to rule out any major discrepancies that could potentially affect the calculation of the design streamflow value and subsequent sizing of hydraulic structures.
For the department of Guajira, P 24h-max isohyetals ranged from 50-120 mm for two years, 70-150 mm for 5 years, 90-160 mm for 10 years, from 100-180 mm for 20 years, 120-200 mm for 25 years, from 120-240 mm for 50 years, and from 120-280 mm for 100 years. The lowest and highest values were observed in the northeastern and northwestern areas, respectively.
For the department of Magdalena, P 24h-max isohyetals ranged from 80-130 mm for two years, 90-160 mm for five years, 110-170 mm for 10 years, from 120-180 mm for 20 years, 120-200 mm for 25 years, from 120-240 mm for 50 years, and from 120-260 mm for 100 years. The lowest and highest values were observed in the eastern and northern areas, respectively.
For the department of Cesar, P 24h-max isohyetals ranged from 80-120 mm for two years, 90-150 mm for five years, 100-170 mm for 10 years, from 120-180 mm for 20 years, 120-200 mm for 25 years, from 120-220 mm for 50 years, and from 120-260 mm for 100 years. The lowest and highest values were observed in the northeastern and southern areas, respectively.
For the department of Atlántico, P 24h-max isohyetals ranged from 70-90 mm for two years, 100-110 mm for five years, 100-120 mm for 10 years, from 120-140 mm for 20 years, 120-140 mm for 25 years, from 140-160 mm for 50 years, and from 120-180 mm for 100 years. The lowest and highest values were observed in the central-eastern and southern areas, respectively. This was the department with the smallest ranges of P 24h-max values mainly due to its total area (3385.1 km 2 ) and topography.
For the department of Bolívar, P 24h-max isohyetals ranged from 90-140 mm for two years, 90-170 mm for five years, 100-190 mm for 10 years, from 120-200 mm for 20 years, 120-200 mm for 25 years, from 120-240 mm for 50 years, and from 140-270 mm for 100 years. The lowest and highest values were observed in the northern and southwestern areas, respectively.
For the department of Sucre, P 24h-max isohyetals ranged from 90-130 mm for two years, 90-160 mm for five years, 100-190 mm for 10 years, from 120-200 mm for 20 years, 120-200 mm for 25 years, from 120-240 mm for 50 years, and from 140-250 mm for 100 years. The lowest and highest values were observed in the northeastern and southeastern areas, respectively.
For the department of Córdoba, P 24h-max isohyetals ranged from 80-130 mm for two years, 100-160 mm for five years, 110-190 mm for 10 years, from 120-200 mm for 20 years, 120-200 mm for 25 years, from 120-220 mm for 50 years, and from 140-240 mm for 100 years. The lowest and highest values were observed in the northwestern and southeastern areas, respectively. Water 2019, 11, x FOR PEER REVIEW 18 of 29

Conclusions
Unlike other countries, Colombia, currently, does not have any official document that compiles a series of recommended methodologies for frequency analysis for a particular region. In this context and contrary to what is the common practice, the Gumbel distribution was not the best fit for most of the time series analyzed. Instead, and according to the chi-squared test, the GEV distribution was

Conclusions
Unlike other countries, Colombia, currently, does not have any official document that compiles a series of recommended methodologies for frequency analysis for a particular region. In this context and contrary to what is the common practice, the Gumbel distribution was not the best fit for most of the time series analyzed. Instead, and according to the chi-squared test, the GEV distribution was shown to be the best fit among the three CDFs used in majority of the datasets. Only 34.3% of the rain gauges fit the Gumbel distribution, while 47.2% of them fit the GEV distribution. LP3, on the other hand, did not work well among the rain gauges analyzed. Based on the results of this study, GEV and Gumbel are the most recommended distributions for the Caribbean region of Colombia.
With respect to the best interpolation method for generating isohyetals in the study area, the IDW method outperformed both the spline and the ordinary kriging methods. These results demonstrate that geostatistically-based interpolation methods (e.g., ordinary kriging) are not always the best selection as many typically take for granted.
According to the results of the REr, RMSE, MDv, and NSC, the areal P 24h-max estimated with the resulting isohyetal maps (by the IDW adjusted method) of this study were close to the true value and can be used to select a design rainfall for a given return period. This is of particular significance given the large amount of ungauged areas within the Colombian Caribbean region where many hydrological studies (made by either international or local consulting companies) are based on the utilization of neighboring rain gauges that sometimes are located at remote distances that do not guarantee the reliability and/or accuracy of the derived design rainfall. The isohyetal maps generated in this study are the first step in developing similar ones for the remaining four regions of Colombia. The 318 multiannual P 24h-max time series analyzed were assumed to be stationary. However, increasing and decreasing trends were observed in some of the time series, suggesting the presence of non-stationarity, which, if confirmed, could result in higher or lower values of point P 24h-max (at the rain gauge). The possible modifications of the isohyetal alignments in certain areas may or may not alter the total areal P 24h-max of a given watershed, since, unlike point rainfall, the areal rainfall is calculated differently. It depends on both the isohyetal alignment (which, in this case, will vary based on the stationary versus the non-stationary P 24h-max value obtained) and the watershed size and distance from the nearest rain gauges, which define how much area of the watershed is covered by each isohyetal (the weighted area). In other words, a change in the value of point rainfall does not necessarily imply a major change in the areal rainfall. Because of that, the authors' future work will address: (a) P 24h-max regionalization at each of the seven departments of the Colombian Caribbean region, (b) the monotonic trend of the multiannual P 24h-max time series via Mann-Kendall and Spearman's rho tests, and (c) the stationarity and non-stationarity of the time series of the rain gauges analyzed in this study so as to compare point versus areal rainfall values and their impact on the design rainfall selection for different return periods at the local and regional level.
Finally, rather than being considered as the sole source for P 24h-max estimation, these maps are intended to be used as a reference in the hydrological and hydraulic analysis, mainly for stormwater management and flood mitigation projects.