Effect of the Non-Stationarity of Rainfall Events on the Design of Hydraulic Structures for Runoff Management and Its Applications to a Case Study at Gordo Creek Watershed in Cartagena de Indias , Colombia

The 24-h maximum rainfall (P24h-max) observations recorded at the synoptic weather station of Rafael Núñez airport (Cartagena de Indias, Colombia) were analyzed, and a linear increasing trend over time was identified. It was also noticed that the occurrence of the rainfall value (over the years of record) for a return period of 10 years under stationary conditions (148.1 mm) increased, which evidences a change in rainfall patterns. In these cases, the typical stationary frequency analysis is unable to capture such a change. So, in order to further evaluate rainfall observations, frequency analyses of P24h-max for stationary and non-stationary conditions were carried out (by using the generalized extreme value distribution). The goodness-of-fit test of Akaike Information Criterion (AIC), with values of 753.3721 and 747.5103 for stationary and non-stationary conditions respectively, showed that the latter best depicts the increasing rainfall pattern. Values of rainfall were later estimated for different return periods (2, 5, 10, 25, 50, and 100 years) to quantify the increase (non-stationary versus stationary condition), which ranged 6% to 12% for return periods from 5 years to 100 years, and 44% for a 2-year return period. The effect of these findings were tested in the Gordo creek watershed by first calculating the resulting direct surface runoff (DSR) for various return periods, and then modeling the hydraulic behavior of the downstream area (composed of a 178.5-m creek’s reach and an existing box-culvert located at the watershed outlet) that undergoes flooding events every year. The resulting DSR increase oscillated between 8% and 19% for return periods from 5 to 100 years, and 77% for a 2-year return period when the non-stationary and stationary scenarios were compared. The results of this study shed light upon to the precautions that designers should take when selecting a design, based upon rainfall observed, as it may result in an underestimation of both the direct surface runoff and the size of the hydraulic structures for runoff and flood management throughout the city.


Introduction
Climate change is gradually coming to affect all living species by giving rise to new sets of meteorological conditions: record breaking high temperatures, melting glaciers, increasing sea water levels, increasing severity of droughts, floods, tornadoes, and hurricanes, to mention just a few.Hydraulic structures have not been the exception; most of these structures were designed was based on rainfall patterns that, in some cases, no longer conform to the realities of today's weather conditions.Generally speaking, sizing a hydraulic structure for an ungauged watershed starts with a hydrological analysis that includes: (a) a rainfall frequency analysis; (b) the watershed delineation and estimation of its morphometric parameters; and (c) the direct surface runoff (DSR) calculation using a rainfall-runoff model, assuming it has the same return period of the generating rainfall.
Usually, a rainfall frequency analysis for an ungauged watershed consists of trying to find the probability distribution that best represents the behavior of the extreme events (of different durations) being analyzed, in order to subsequently obtain the value of the precipitation event associated with a given return period to be later used for the DSR estimation.The traditional methods for frequency analysis of extreme values for hydraulic structures design (streamflow, rainfall or rainfall intensity) assume that: (a) the values of the variable are random and independent; (b) the probability distribution functions are stationary; and (c) the return period (Tr) is the inverse of the probability of exceedance (p).These assumptions imply that the probability of an event producing excessive precipitations does not change over time (the sample is homogeneous).This approach cannot be applied when the variable being analyzed shows noticeable change of pattern (either increase or decline), which may indicate non-stationary conditions [1][2][3][4].In these cases, the typical definition of the return period (Tr = 1/p) has to be redefined to account for non-stationarity.
Cartagena de Indias, a city on the Caribbean coast of Colombia, has been lately undergoing recurring floods that can be mainly attributed to different factors that go from uncontrolled development projects-due to the lack of an updated territorial management plan that fits the city's realities-to poorly designed hydraulic structures for runoff management that do not take future conditions (especially rainfall pattern changes and increase of impervious areas) into account.Most of the runoff of Cartagena de Indias discharges into La Virgen swamp (ciénaga de La Virgen), a waterbody of approximately 502 km 2 , whose southern shoreline has been populated over the years by illegal low-income settlements that suffer the consequences of both a rising sea level (the swamp is connected to the Caribbean Sea) and an obsolete stromwater system.Furthermore, the cities of Cartagena de Indias (downstream) and Turbaco (upstream) share several watersheds that also drain into the La Virgen swamp.The upstream areas of these watersheds (like the one selected in this study) are mostly rural, which have been gradually converted into impervious areas by local developers that offer more competitive prices than those in Cartagena de Indias.This dynamic will most likely (and rapidly) change the landscape, which brings with it more challenging hydrologic, and hydraulic, conditions.
In this context, Colombian legislation has mandated to analyze the possible effects of climate change on the future rainfall pattern if a stormwater system (channels, pipes, culverts, etc.) is to be designed [5].A regional increasing trend of both rainfall and streamflow has been found in some areas of the northern portion of the Colombian Caribbean region [6,7].The studies that carry out assessments for stationary (SC) versus non-stationary (NSC) conditions of streamflow and/or rainfall, typically encompass only performing a statistical analysis (of a trend) without quantifying how the two scenarios may affect real life conditions, both hydrologically and hydraulically speaking.In this study, multiannual 24-h maximum rainfall (P 24h-max ) data (from the synoptic weather station located at the Rafael Núñez airport) was used to: (a) assess the trend of the P 24h-max observations over time, specifically for Cartagena de Indias, (b) quantify the rainfall values obtained for several return periods (2,5,10,25,50, and 100 years) via stationary and non-stationary frequency analyses, (c) carry out a hydrological and hydraulic analysis on the Gordo creek watershed under SC and NSC for different return periods in order to understand the recurring floods reported that affect a commercial and

Study Area
The Gordo creek watershed is located within the jurisdiction of the city of Turbaco (Colombia); the outlet is at the east-southern borderline between the cities of Cartagena de Indias and Turbaco (Figure 1).It has a total area of 258 hectares.The creek runs south to north, with a total length of 2962 m and an average slope of 0.00454.The watershed's centroid and outlet are located at the geographical coordinates 10

Study Area
The Gordo creek watershed is located within the jurisdiction of the city of Turbaco (Colombia); the outlet is at the east-southern borderline between the cities of Cartagena de Indias and Turbaco (Figure 1).It has a total area of 258 hectares.The creek runs south to north, with a total length of 2962 m and an average slope of 0.00454.The watershed's centroid and outlet are located at the geographical coordinates 10°20′04″ N; 75°27′29″ W and 10°22′34.83″N; 75°27′37.51″W respectively.The watershed's soil cover in the upstream area is mostly rural (bush and forest) while the downstream is in a busy commercial, industrial, and residential area (Table 1).The soil texture is predominantly clay-loam with a poor infiltration capacity when saturated.The watershed's downstream area has a main road with a box-culvert that gives access to several industries and businesses.The most critical portion of the watershed's downstream area (circled area in Figure 2) is composed of the last 178.5-mreach (of the creek) and the box-culvert that controls the watershed outlet.Flooding events occur every year during the rainy season (September through November), which impedes regular transit for several hours.Also, people and vehicles have been stranded, and have required rescuing on the flooded road (Figure 3a-d).The watershed's soil cover in the upstream area is mostly rural (bush and forest) while the downstream is in a busy commercial, industrial, and residential area (Table 1).The soil texture is predominantly clay-loam with a poor infiltration capacity when saturated.The watershed's downstream area has a main road with a box-culvert that gives access to several industries and businesses.The most critical portion of the watershed's downstream area (circled area in Figure 2) is composed of the last 178.5-mreach (of the creek) and the box-culvert that controls the watershed outlet.Flooding events occur every year during the rainy season (September through November), which impedes regular transit for several hours.Also, people and vehicles have been stranded, and have required rescuing on the flooded road (Figure 3a-d).

Rainfall Data
Although the Institute of Hydrology, Meteorology, and Environmental Studies of Colombia (IDEAM, in Spanish) operates another two weather stations in Cartagena de Indias (Escuela Naval-CIOH and Santa Ana), the one located at the airport Rafael Núñez has the longest and most complete data of records for several meteorological variables (Table 2); for this reason, it is the most used in hydrological and hydraulic analysis and designs.

Methodology
The methodology used in this research is presented in the flowchart of Figure 4.Each step is described in detail in the following sub-sections.

Rainfall Stationary Frequency Analysis
The stationary frequency analysis is based on the following assumptions: (a) the cumulative distribution functions (CDFs) do not change over time [2] and (b) the variable analyzed must be random and independent [1].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, Log-Pearson 3, among others.The GEV was used for the frequency analysis carried out in this section, as it has shown to work best in Colombia.
The expression for the GEV distribution is given by Equation ( 1), where 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).The GEV distribution may adopt one of the three types of the EV distribution depending upon the value of k [8-10]: (a) when k equals zero, EV is type 1 (Gumbel) [11]; (b) when k is less than zero, EV is type 2 (Frechet) [12]; and (c) when k is greater than zero, EV is type 3 (Weibull) [13].
The estimated GEV's parameters were: α = 31.35mm, β = 78.097mm, and k = 0.06.Table 4 shows the values of the P24h-max for each of the two methods used.Gray cells show maximum values.
Lastly, the CDF goodness-of-fit was carried out via Chi-square [13] and Kolmogorov-Smirnov [14][15][16][17] tests.Values of 0.06001 and 2.2876 were obtained for the Kolmogorov-Smirnov and Chisquare tests respectively.These results indicate that the GEV distribution accurately represents the behavior of the rainfall data analyzed.

Rainfall Stationary Frequency Analysis
The stationary frequency analysis is based on the following assumptions: (a) the cumulative distribution functions (CDFs) do not change over time [2] and (b) the variable analyzed must be random and independent [1].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, Log-Pearson 3, among others.The GEV was used for the frequency analysis carried out in this section, as it has shown to work best in Colombia.
The expression for the GEV distribution is given by Equation ( 1), where 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).The GEV distribution may adopt one of the three types of the EV distribution depending upon the value of k [8-10]: (a) when k equals zero, EV is type 1 (Gumbel) [11]; (b) when k is less than zero, EV is type 2 (Frechet) [12]; and (c) when k is greater than zero, EV is type 3 (Weibull) [13].
The estimated GEV's parameters were: α = 31.35mm, β = 78.097mm, and k = 0.06.Table 4 shows the values of the P 24h-max for each of the two methods used.Gray cells show maximum values.
Lastly, the CDF goodness-of-fit was carried out via Chi-square [13] and Kolmogorov-Smirnov [14][15][16][17] tests.Values of 0.06001 and 2.2876 were obtained for the Kolmogorov-Smirnov and Chi-square tests respectively.These results indicate that the GEV distribution accurately represents the behavior of the rainfall data analyzed.

Rainfall Non-Stationary Frequency Analysis
The non-stationary analysis considers the following assumptions [2][3][4]: (a) the probability of exceeding the random variable over the years (p j ) of the CDFs may increase or decrease, and (b) the extreme values of the analyzed variable are independent in time.Parameters of the CDFs are considered time-varying in order to compute the probability of exceeding occurring, which implies that extreme values are not identically distributed over time.
To define the return period under non-stationary environments (Tr, NSC ), it is necessary to understand the concept of the waiting time.The waiting time is a random variable (X) defined as the occurrence of a value exceeding the design value for the first time.In stationary conditions, the probability of exceeding remains constant over time, which implies that a geometric distribution can be used to compute the expected value (E(x) = Tr, SC = 1/p).Commonly, this is called the return period [8].In contrast, under non-stationary conditions, the return period (Tr) is computed with a non-homogeneous geometric distribution (Equation ( 2)) [18]: where p j is the time-varying probability of exceeding, and the subscript j represents the projecting year.
The GEV distribution with a temporal trend in the location parameter (β t ) can be expressed as (Equation ( 3)) [2,18]: The P 24h-max data at the Rafael Núñez airport weather station (Table 4) exhibits an increasing linear trend over time (dotted black line in Figure 5), which confirms the non-stationary occurrence of the analyzed variable.
To estimate the P 24h-max values, the non-stationary frequency analysis followed these steps [4,18]: (a) selecting the CDF under non-stationary condition (GEV or Gumbel); (b) defining a model (constant, linear, or establishing other models) of each parameter (location, shape, and scale) with its trend (increasing or decreasing); (c) estimating parameters by likelihood method; (d) calculating the goodness-of-fit test of the Akaike Information Criterion (AIC) [19] and then selecting the CFD with the minimum value of it; and (e) estimating the P 24h-max for various return periods.
For this study, the CDF under non-stationary condition (GEV distribution) was selected considering the distribution used in the stationary scenario of Section 3.1.The GEV's parameters were estimated by means of a code [2,18] programmed in R software (version 3.3.1,R Development Core Team, Auckland, New Zealand) with the library nsextremes.The resulting linear trend lines for the location (increasing trend) and scale (remains constant) parameters are shown in Figure 5 (blue and green lines for the location and scale parameters, respectively).The values of the aforementioned parameters were: βt = 79.543mm + [(0.4949 mm/year) (t mean -1978)], α = 31.1 mm, and k = −0.106.The t mean term in the previous expression represents an analyzed year centered over its own mean between 1941 and 2015.A sensitivity analysis showed that variations of both the scale and the shape parameters do not bring an adequate trend for the P 24h-max when the results were compared with the findings of IDEAM.The AIC test was used in order to evaluate the GEV distribution goodness-of-fit by comparing both scenarios (stationary and non-stationary).The obtained values were 747.5103 (non-stationary) and 753.3721 (stationary).The results of the P 24h-max for several return periods are presented in Table 5.To estimate the P24h-max values, the non-stationary frequency analysis followed these steps [4,18]: (a) selecting the CDF under non-stationary condition (GEV or Gumbel); (b) defining a model (constant, linear, or establishing other models) of each parameter (location, shape, and scale) with its trend (increasing or decreasing); (c) estimating parameters by likelihood method; (d) calculating the goodness-of-fit test of the Akaike Information Criterion (AIC) [19] and then selecting the CFD with the minimum value of it; and (e) estimating the P24h-max for various return periods.
For this study, the CDF under non-stationary condition (GEV distribution) was selected considering the distribution used in the stationary scenario of Section 3.1.The GEV's parameters were estimated by means of a code [2,18] programmed in R software (version 3.3.1,R Development Core Team, Auckland, New Zealand) with the library nsextremes.The resulting linear trend lines for the location (increasing trend) and scale (remains constant) parameters are shown in Figure 5 (blue and green lines for the location and scale parameters, respectively).The values of the aforementioned parameters were: βt = 79.543mm + [(0.4949 mm/year) (tmean -1978)], α = 31.1 mm, and k = −0.106.The tmean term in the previous expression represents an analyzed year centered over its own mean between 1941 and 2015.A sensitivity analysis showed that variations of both the scale and the shape parameters do not bring an adequate trend for the P24h-max when the results were compared with the findings of IDEAM.The AIC test was used in order to evaluate the GEV distribution goodness-of-fit by comparing both scenarios (stationary and non-stationary).The obtained values were 747.5103 (non-stationary) and 753.3721 (stationary).The results of the P24h-max for several return periods are presented in Table 5.

Curve Number (CN) Estimation
Given that the Gordo creek watershed is composed of several soil cover types, an area-based composite CN was estimated following the methodology established in the National Engineering Handbook, Section 4 (NEH-4) [20], where a series of lookup tables are given for different combinations of hydrological soil groups (HSGs) and soil cover.The HSG of the watershed is C, as per the description of the soil texture provided in Section 2.1 (clay-loam soil texture).The values shown in those tables represent what is called the average antecedent runoff condition (ARC-II).The ARC is a concept introduced into the CN methodology to try to account for the potential for runoff generation in a watershed due to the soil water content after several rainfall events [21].For that, the total amount of cumulative rainfall of the 5 previous days (5-day antecedent precipitation, P 5 ) is estimated to establish whether the soil is in dry, average/normal or wet (Table 6) [21]; a set of equations are then proposed to adjust the CN accordingly (Equations ( 4) and ( 5)) [8,22].
CN III = CN II 0.427 − 0.00573CN II (5)  In 2010-2011, Colombia underwent the longest and most devastating rainy season ever recorded.Floods were reported throughout the country, especially in the Caribbean coastal region.In this study, a CN adjustment was carried out for ARC III to mimic a saturated soil as the most critical scenario in the Gordo creek watershed in that season, during which more than 53 mm of rain fell in 5 consecutive days.Table 7 shows the values of CN composite obtained for average and wet conditions (gray cells).

Time of Concentration Estimation
The time of concentration may be defined as the time it takes a drop of water to travel from the hydraulically most remote point to the watershed outlet (or any other point of interest) [23].It is one of the most critical parameters for any hydrological analysis and the subsequent design of a hydraulic structure aimed at managing stormwater.Many methodologies (empirical and semi-empirical) have been proposed for its estimation [24,25].For this study, the Kirpich method was used (Equation ( 6)) [26], as the Gordo creek watershed fits some of the morphometric characteristics needed for using this formula [25].
In Equation ( 6), Tc is the time of concentration (min), L is the length of the main watercourse/channel (m) and S is the average slope of the main watercourse/channel (m/m).After plugging the L and S values provided in Section 2.1, the Tc equals 73.3 min.

Direct Surface Runoff (DSR) Estimation
The DSRs produced by the rainfall under both stationary and non-stationary conditions for different return periods (2, 5, 10, 25, 50, and 100 years) were estimated (Table 8) via Hydrologic Engineering Center Hydrologic Modeling System (HEC-HMS), which is a CN-based hydrological model developed by the U.S. Corps Army of Engineers [27].HEC-HMS compiles several methodologies for hydrological analysis integrated into one tool capable of quickly modeling different scenarios at the same time.
Apart from the adjusted CN value estimated in Section 3.3, the HEC-HMS simulations for the various return periods were made under the following conditions: (a) 3-h rainfall as it is the most Fluids 2018, 3, 27 10 of 19 typical duration historically observed in Cartagena de Indias; (b) the 3-h rainfall corresponds to 79% of the P 24h-max , and its distribution was estimated via multiannual probabilistic analysis of 30 pluviographs [8] (Figure 6) (P90% was the one used; it indicates the probability that the rainfall pattern observed falls to the left of the P90% curve); (c) no rainfall area reduction was needed due to the size of the watershed; (d) hydrograph generation via SCS-unit hydrograph method; and (e) a lag time of 0.6 Tc (Tlag = 0.6Tc).Table 8 summarizes the peak flow values for stationary and non-stationary conditions.
different return periods (2, 5, 10, 25, 50, and 100 years) were estimated (Table 8) via Hydrologic Engineering Center Hydrologic Modeling System (HEC-HMS), which is a CN-based hydrological model developed by the U.S. Corps Army of Engineers [27].HEC-HMS compiles several methodologies for hydrological analysis integrated into one tool capable of quickly modeling different scenarios at the same time.
Apart from the adjusted CN value estimated in Section 3.3, the HEC-HMS simulations for the various return periods were made under the following conditions: (a) 3-h rainfall as it is the most typical duration historically observed in Cartagena de Indias; (b) the 3-h rainfall corresponds to 79% of the P24h-max, and its distribution was estimated via multiannual probabilistic analysis of 30 pluviographs [8] (Figure 6) (P90% was the one used; it indicates the probability that the rainfall pattern observed falls to the left of the P90% curve); (c) no rainfall area reduction was needed due to the size of the watershed; (d) hydrograph generation via SCS-unit hydrograph method; and (e) a lag time of 0.6 Tc (Tlag = 0.6Tc).Table 8 summarizes the peak flow values for stationary and nonstationary conditions.53.9 60.9

Channel Hydraulic Modeling
A one-dimension (1-D) hydraulic modeling was carried out using the Hydrologic Engineering Center River Analysis System (HEC-RAS) software, developed by the U.S. Corps Army of Engineers.The simulation was performed under gradually varied flow and subcritical conditions (Froude number less than 1).Also, the following inputs were utilized in the simulation: (a) Gordo creek's bathymetry, that includes ten cross-sections (Figure 7) of the last 178.5 m of the stream (the most critical portion of the watershed) which consist of a natural cross-section with three culverts; (b) design DSR

Channel Hydraulic Modeling
A one-dimension (1-D) hydraulic modeling was carried out using the Hydrologic Engineering Center River Analysis System (HEC-RAS) software, developed by the U.S. Corps Army of Engineers.The simulation was performed under gradually varied flow and subcritical conditions (Froude number less than 1).Also, the following inputs were utilized in the simulation: (a) Gordo creek's bathymetry, that includes ten cross-sections (Figure 7) of the last 178.5 m of the stream (the most critical portion of the watershed) which consist of a natural cross-section with three culverts; (b) design DSR for several return periods estimated in Section 3.5; (c) hydraulic control section governed by a subcritical flow at the downstream box culvert (section K0+000), which produces a backwater profile in the entire channel.The natural channel exhibits different cross-sections, as presented in Figure 8a,b.
Figure 9a,b show the plant view and the longitudinal profile of Gordo creek watershed outlet (gray-shaded areas show the location of culverts).A description of the culverts is presented in Table 9.The entire natural channel presents a Manning's coefficient of 0.050 (natural minor streams with some weeds and stone), while culverts have a value of 0.014 (concrete in culverts with bends, connections, and some debris) [8].for several return periods estimated in Section 3.5; (c) hydraulic control section governed by a subcritical flow at the downstream box culvert (section K0+000), which produces a backwater profile in the entire channel.The natural channel exhibits different cross-sections, as presented in Figure 8a,b.Figure 9a,b show the plant view and the longitudinal profile of Gordo creek watershed outlet (gray-shaded areas show the location of culverts).A description of the culverts is presented in Table 9.The entire natural channel presents a Manning's coefficient of 0.050 (natural minor streams with some weeds and stone), while culverts have a value of 0.014 (concrete in culverts with bends, connections, and some debris) [8].for several return periods estimated in Section 3.5; (c) hydraulic control section governed by a subcritical flow at the downstream box culvert (section K0+000), which produces a backwater profile in the entire channel.The natural channel exhibits different cross-sections, as presented in Figure 8a,b.Figure 9a,b show the plant view and the longitudinal profile of Gordo creek watershed outlet (gray-shaded areas show the location of culverts).A description of the culverts is presented in Table 9.The entire natural channel presents a Manning's coefficient of 0.050 (natural minor streams with some weeds and stone), while culverts have a value of 0.014 (concrete in culverts with bends, connections, and some debris) [8].For the peak flow, the values of the overall increase ratio have the same behavior of that shown by P 24h-max , with the maximum value occurring at a 2-year return period.A 44-percent increase in rainfall (Table 10) caused a 77-percent peak flow rise (Table 11).Once again, this indicates the precautions that the designers must take when using this return period.For the other return periods, it may be observed that the increase in P 24h-max resulted in a peak flow increase of a similar range.Figures 11 and 12 show respectively the non-stationary and stationary difference and ratios between P 24h-max and peak flow for different return periods (DPmax, DQp, Pratio and Qratio).As expected, both show how the behavior of P 24h-max , and the watershed's response to it (the peak flow), have the same trend.For the peak flow, the values of the overall increase ratio have the same behavior of that shown by P24h-max, with the maximum value occurring at a 2-year return period.A 44-percent increase in rainfall (Table 10) caused a 77-percent peak flow rise (Table 11).Once again, this indicates the precautions that the designers must take when using this return period.For the other return periods, it may be observed that the increase in P24h-max resulted in a peak flow increase of a similar range.The water level reached after the simulations carried out under stationary and non-stationary conditions indicates that modeling under non-stationary scenario better represents the current hydrological conditions of the Gordo creek watershed, where rainfall events of low magnitude are generating floods every year during the rainy season (during which the soil is saturated most of the time) at the watershed's outlet, when, in theory, they ought not to.This can be observed, for instance, The water level reached after the simulations carried out under stationary and non-stationary conditions indicates that modeling under non-stationary scenario better represents the current hydrological conditions of the Gordo creek watershed, where rainfall events of low magnitude are generating floods every year during the rainy season (during which the soil is saturated most of the time) at the watershed's outlet, when, in theory, they ought not to.This can be observed, for instance, in the results of the hydraulic simulation at K0+018.5 shown in Table 12 (P 24h-max , peak flow, and water elevation) and Figure 13, where the bankfull level (9.94 m) is reached at a lower return period (floods occur more frequently) under non-stationary conditions (Tr < 2-year).Under stationary conditions, a 2-year return period flow, defined as the mean annual flow (or annual maximum daily flow) [28,29], did not result in a bankfull section, which is not what is currently occurring in the study area (floods every year).The change in the rainfall pattern evidenced in previous sections may be the reason for a shift towards more recurrent and higher-than-usual events that cause floods more frequently than in the past.Likewise, statistically speaking, the ACI results obtained in Section 3.2 (747.5103 for non-stationary and 753.3721 for stationary) demonstrated that the non-stationary condition is more adequate.The AIC test serves to evaluate how good a model is for predicting future values [30].The lower the value of the AIC the better the model.In terms of the flood depth (elevation above bankfull level) at the area of the watershed outlet, Table 13 and Figure 14 show the water depth above the bankfull level for different return periods.This sensitivity analysis reveals the necessity for sizing hydraulic structures under non-stationary conditions to avoid more recurring floods.Notwithstanding that the 1-D hydraulic simulation is incapable of delineating the extension of the flooded area beyond the creek's bank, the difference between NSC and SC elevations in Table 13 is an indication of what locals and workers of the industrial area have been noticing every year: floods in this area have been gradually worsening in terms of both frequency of occurrence and area covered by water (as well as the water depth marks In terms of the flood depth (elevation above bankfull level) at the area of the watershed outlet, Table 13 and Figure 14 show the water depth above the bankfull level for different return periods.This sensitivity analysis reveals the necessity for sizing hydraulic structures under non-stationary Fluids 2018, 3, 27 conditions to avoid more recurring floods.Notwithstanding that the 1-D hydraulic simulation is incapable of delineating the extension of the flooded area beyond the creek's bank, the difference between NSC and SC elevations in Table 13 is an indication of what locals and workers of the industrial area have been noticing every year: floods in this area have been gradually worsening in terms of both frequency of occurrence and area covered by water (as well as the water depth marks left at some points that local people use as reference points).Further research is needed, though.Future work in this area must include a more detailed survey (topographic study) covering more area, so as to assess and more accurately quantify (numerically and spatially) the aforementioned anecdotal evidence.In addition to the already quantified effects of climate change on Cartagena's P24h-max pattern and the Gordo creek watershed's hydrological response, the fact that almost 85% of the watershed is still rural indicates that any increase in the impervious area will both raise the peak flow and reduce the time of concentration.These two variables will worsen the situation at the outlet of Gordo creek watershed unless a series of measures are implemented.For instance, the combination of sustainable urban drainage systems (SUDS) [31,32], stormwater storage vaults/tanks (both online and offline), aquifer storage and recovery (ASR) wells, infiltration ponds, among others, have proved to be effective at managing stormwater by: (a) keeping the peak flow at the same magnitude (or lower) when pre-development and post-development scenarios are compared, (b) avoiding the design of larger hydraulic structures, (c) minimizing the risk of floods, and (d) recharging aquifers, as a plus.

Conclusions
P24h-max values associated to a given return period are a key variable when estimating design DSR for hydraulic structures for stormwater management, especially in ungauged watersheds.The values obtained in the AIC test carried out in this study, indicated that a frequency analysis under nonstationary conditions represents best the behavior of the rainfall patterns of the P24h-max observations at the Rafael Núñez station.A stationary scenario is, in this case, further from the natural reality, when compared to non-stationary conditions.This was also confirmed by both the increase in the occurrence of P24h-max events greater than or equal to 148.1 mm (value for a 10-year return period An option for making decisions when designing some hydraulic structures would be to perform a benefit-cost analysis to compare a design under non-stationary conditions versus under stationary condition, taking into account free board heights, allowing handling of the non-stationary peak flow depending on the design return period. In addition to the already quantified effects of climate change on Cartagena's P 24h-max pattern and the Gordo creek watershed's hydrological response, the fact that almost 85% of the watershed is still rural indicates that any increase in the impervious area will both raise the peak flow and reduce the time of concentration.These two variables will worsen the situation at the outlet of Gordo creek watershed unless a series of measures are implemented.For instance, the combination of sustainable urban drainage systems (SUDS) [31,32], stormwater storage vaults/tanks (both online and offline), aquifer storage and recovery (ASR) wells, infiltration ponds, among others, have proved to be effective at managing stormwater by: (a) keeping the peak flow at the same magnitude (or lower) when pre-development and post-development scenarios are compared, (b) avoiding the design of larger hydraulic structures, (c) minimizing the risk of floods, and (d) recharging aquifers, as a plus.

Conclusions
P 24h-max values associated to a given return period are a key variable when estimating design DSR for hydraulic structures for stormwater management, especially in ungauged watersheds.The values obtained in the AIC test carried out in this study, indicated that a frequency analysis under non-stationary conditions represents best the behavior of the rainfall patterns of the P 24h-max observations at the Rafael Núñez station.A stationary scenario is, in this case, further from the natural reality, when compared to non-stationary conditions.This was also confirmed by both the increase in the occurrence of P 24h-max events greater than or equal to 148.1 mm (value for a 10-year return period under stationary conditions), and the increasing linear trend over time of the overall P 24h-max observations.These findings may be an indication that the typical frequency analysis of rainfall under stationary conditions may no longer be applicable when calculating the design rainfall associated with a given return period for sizing hydraulic structures throughout Cartagena de Indias.Furthermore, designing under stationary conditions may have direct implications in the decisions local authorities will be taking in the coming years, given that millions of dollars will be invested in upgrading the city's stormwater system.The 1-D hydraulic simulation performed herein revealed that the Gordo creek watershed outlet area can be flooded even with an event of 2-year return period (every year)-a situation that had never been observed in the past according to what local people have affirmed.The increase in the P 24h-max , though, is not the sole factor to be taken into account when evaluating the reasons behind the more frequent floods reported.Uncontrolled and unplanned urbanization of vegetated areas may and will exacerbate the recurrent floods registered within the study area.
It is noteworthy to mention that sizing larger hydraulic structures for stormwater management as the only solution for the larger rainfall events estimated under non-stationary conditions might be unsuitable and unsustainable in cases where surface area limitations exist (especially in urban areas).Therefore, implementing best management practices for controlling stormwater sources in old and newly constructed areas shall be one of the city's goals in maintaining the design peak flow values, especially under increasing non-stationary conditions.

Fluids
downstream, and (d) point out the importance of accounting for the effect of climate change in the decision making process in runoff management, especially in flood-prone areas.

Figure 5 .
Figure 5. P 24h-max and GEV parameters (location and scale) behavior over time (Rafael Núñez airport station).

Figure 6 .
Figure 6.3-h rainfall distribution pattern at Rafael Núñez airport weather station.

Figure 6 .
Figure 6.3-h rainfall distribution pattern at Rafael Núñez airport weather station.

Figure 7 .
Figure 7. Plant view of the critical point of flooding.

Figure 7 .
Figure 7. Plant view of the critical point of flooding.

Figure 7 .
Figure 7. Plant view of the critical point of flooding.

Figure 9 .
Figure 9. Hydraulic modeling with HEC-RAS: (a) Plan view; (b) longitudinal profile.WS indicates the water surface level.

Figure 9 .
Figure 9. Hydraulic modeling with HEC-RAS: (a) Plan view; (b) longitudinal profile.WS indicates the water surface level.

Figure 9 .
Figure 9. Hydraulic modeling with HEC-RAS: (a) Plan view; (b) longitudinal profile.WS indicates the water surface level.

Figure 10 .
Figure 10.P24h-max and peak flow increase for various return periods.Figure 10.P 24h-max and peak flow increase for various return periods.

Figure 10 .
Figure 10.P24h-max and peak flow increase for various return periods.Figure 10.P 24h-max and peak flow increase for various return periods.
Figures 11 and 12  show respectively the non-stationary and stationary difference and ratios between P24h-max and peak flow for different return periods (DPmax, DQp, Pratio and Qratio).As expected, both show how the behavior of P24h-max, and the watershed's response to it (the peak flow), have the same trend.

Figure 11 .
Figure 11.P24h-max and peak flow variation for various return periods.DPmax and DQp are the P24h-max and peak flow difference between non-stationary and stationary conditions.Figure 11.P 24h-max and peak flow variation for various return periods.DPmax and DQp are the P 24h-max and peak flow difference between non-stationary and stationary conditions.

Figure 11 .
Figure 11.P24h-max and peak flow variation for various return periods.DPmax and DQp are the P24h-max and peak flow difference between non-stationary and stationary conditions.Figure 11.P 24h-max and peak flow variation for various return periods.DPmax and DQp are the P 24h-max and peak flow difference between non-stationary and stationary conditions.

Fluids 2018, 3 , 19 Figure 12 .
Figure 12.P24h-max and peak flow ratio variation for various return periods.Qratio and Pratio are the ratio of non-stationary to stationary values of peak flow and P24h-max respectively.

Figure 12 .
Figure 12.P 24h-max and peak flow ratio variation for various return periods.Qratio and Pratio are the ratio of non-stationary to stationary values of peak flow and P 24h-max respectively.

Figure 13 .
Figure 13.Water level at cross-section K0+018.5 for different return periods under stationary and nonstationary conditions.

Figure 13 .
Figure 13.Water level at cross-section K0+018.5 for different return periods under stationary and non-stationary conditions.

Table 2 .
Description of Rafael Núñez airport weather station.

Table 3 ;
gray cell is the maximum value registered) were used.Years 2016 and 2017 have not yet been officially reported by IDEAM.

Table 4 .
P 24h-max values for stationary freq.analysis.Bold values indicate the largest value of the two CDF used.

Table 6 .
ARC types based on P 5 .

Table 7 .
Composite CN for Gordo creek watershed.

Table 9 .
Description of culverts.

Table 11 .
Peak flow values comparison.

Table 11 .
Peak flow values comparison.

Is at 9.94 m Max. Water Level Elevation Reached at Peak Flow
Values in bold indicate which year (return period) first exceeded the bankfull level under SC and NSC.
Values in bold indicate which return period first exceeded the bankfull level under SC and NSC.Fluids 2018, 3, x FOR PEER REVIEW 17 of 19