Non-stationary analysis for road drainage design under land-use and climate change scenarios.

Frequency analysis has been the most widely used tool worldwide to dimension water-related infrastructures and evaluate flood risks. The concept of stationarity has been a common and practical hypothesis in hydrology for many years. However, in recent decades due to climate change pressure and changes in land use, it has been related to the presence of time-series trends that in hydrology indicate non-stationary effects. In this sense, the need to comprehensively address non-stationary frequency analysis has been identified. This study proposes to incorporate the non-stationary flood frequency analysis into the dimensioning process of road structures with the following objectives: i) evaluate the effect of land use on peak flow in a simulated period of 129 years, ii) evaluate covariates related to land use, and iii) evaluate covariates related to climate change. To this end, road drainage simulation exercises were carried out in three sections of the Ibagué-Cajamarca road located in Colombia. Likewise, the Generalized Additive Models for Location, Scale and Shape was implemented for the non-stationary analysis, and covariates related to climate variability were included, such as El Niño-Southern Oscillation indices (ONI12, ONI3.4, MEI, and SOI), and the Pacific Decadal Oscillation (PDO) index, as well as some related to the evolution of land use such as hydraulic conductivity, soil water storage in the root zone, and infiltration capacity represented in the curve number. The results indicate that the non-stationary analysis improves the prediction of maximum flows, and it is possible to obtain road drainage dimensioning that adjusts to climate and land-use variations.

Introduction

In recent decades, flood frequency analysis with a non-stationary approach has been the subject of extensive debate in the hydrological scientific community. It has positioned itself as one of the main topics of analysis since the study of Milly et al. (2008). Accordingly, some researchers state that assuming a stationarity premise should not be discarded, and it is a rational and useful framework to project hydrological risk in the future (Serinaldi and Kilsby, 2015; Montanari and Koutsoyiannis, 2014). On the other hand, others consider it appropriate to introduce the non-stationarity hypothesis as a strategy to improve the projection of hydrological risk (e.g., Ilić et al., 2021; Matalas, 2012; Milly et al., 2008), given the potential climate change and the land-use effects on the maximum flow regime (Sadeghi et al., 2020; Mwangi et al., 2016; Zhang et al., 2015).

The traditional design process for road drainage structures involves estimating direct runoff for a projected design period from observed surface stream flows or using records of maximum daily precipitation in the study area. To this end, the stationary flood frequency analysis has broadly been used to project the hydrological risk of structures (Ul et al., 2019) and drainage structures design (e.g., Department of Public Works, 2017; INVIAS, 2009). In this sense, it has been recognized that climate change introduces variations in the frequency and intensity of rainfall (e.g., Li and Fang, 2017; Narsimlu et al., 2013), as well as effects related to the evolution of land use on the magnitude of surface runoff (e.g., Omer et al., 2020; Siswanto and Francés, 2019). Therefore, the need to address the flood frequency analysis under the assumption of non-stationarity has been identified (Nasr et al., 2019; Salas et al., 2018), given the alterations that climate change and the evolution of land use produce on the original design conditions (Nasr et al., 2019; Ehsani et al., 2017) and, consequently, generate scenarios of probable failures on the capacity of hydraulic structures (Hui et al., 2018; Mondal and Mujumdar, 2015; Salas and Obeysekera, 2014).

In this context, the implementation of statistical models that allow the adjustment of the parameters of the distribution functions through covariates that explain the trends in the hydrological series can improve the results of the flood frequency analysis and the estimation of the design flow (e.g., Agilan and Umamahesh, 2017; Nasri et al., 2017; Šraj et al., 2016). In the literature, various techniques have been developed to carry out non-stationary flood frequency analysis (Cannon, 2010; Villarini et al., 2010) and have mainly been applied to the design of hydraulic structures (e.g., Nasri et al., 2017; Mondal and Mujumdar, 2015). However, the academic literature has not reported NSFFA applications for the dimensioning process of road drainage structures.

Accordingly, the aim of this work is to evaluate the non-stationary effect of land use and climate changes on the capacity of road drainage structures through a case study analysis located in the Ibagué-Cajamarca road in Colombia. Likewise, the results of the stationary analysis of maximum flows are evaluated with respect to the non-stationary analysis, incorporating climate change scenarios in the hydrological time series projected until the year 2100 (IPCC, 2014; IDEAM et al., 2015). Thus, covariates related to land use were considered, such as soil water storage capacity in the root zone (Hu), vertical saturated hydraulic conductivity (Ks), and infiltration capacity represented in the curve number (CN) (USDA-SCS, 1972). Also, the climatic indices Multivariate ENSO Index (MEI), Southern Oscillation Index (ONI12, ONI3.4), Southern Oscillation Index (SOI), and Pacific Decadal Oscillation (PDO) were evaluated as covariates related to climate variability. Therefore, the evaluation of the capacity of hydraulic structures in a variable land-use and climate change environment was carried out by applying a non-stationary flood frequency analysis using generalized additive models for location, scale, and shape (GAMLSS). The results indicate that it is possible to obtain better maximum flow predictions when a non-stationary analysis is incorporated, so such analysis can contribute to achieving hydraulic designs more adjusted to the changing conditions of the hydrographic basins.

Methodology

Study area

The capacity of three road drainage structures located at the abscissas K4+940, K8+838, and K11 + 010 of the Ibagué-Armenia road in Colombia was evaluated and identified in this study as E1, E2, and E3, respectively (Figure 1). An average of 2,157,352 vehicles per year transit through this road and 20,298.05 tons/year of cargo is transported. Daily rainfall records for the last 40 years from the meteorological stations 21210160 and 21245040 managed by the Institute of Hydrology, Meteorology and Environmental Studies (IDEAM) indicate a bimodal precipitation behavior in the study area. The periods with the highest amount of rainfall are from November to January (NDJ) and from March to May (MAM), while the dry periods, with a reduction in rainfall events, are from June to August (JJA) and from September to November (SON).

Figure 1

Location map of hydraulic structures. E1, E2, and E3 are the drainage structures with catchment area; MQD: soil with high infiltration capacity and coarse texture; MQO: soil with fine to medium texture and moderate infiltration capacity; Hu: soil water storage capacity in the root zone, Ks: saturated hydraulic conductivity, and CN: curve number.

The afferent areas that contribute with runoff to structures E1, E2, and E3 were delineated using GIS techniques and a digital elevation model. The curve number (CN) assignment requires determining the soil types and land uses in the watersheds (USDA-SCS 1972). Figure 1 shows the soil types in the study watersheds, type A (MQD) and B (MQO). The percentages of vegetation cover in the study area (e.g., Márquez et al., 2021), terrain slope, and soil type for each watershed are shown in Table 1. The curve number generation process was carried out based on the USDA-SCS manual (1972). The classification of plant cover types was carried out based on the available land-use maps of the study zone corresponding to the years 2006, 2013, and 2018, which includes forest, crops, grassland, and impervious surfaces. Table 1 presents the main characteristics of the structures assessed and their corresponding afferent areas.

Table 1

Main characteristics of the road drainage structures analyzed.

Structure	Shape	Geometry	Material	Roughness coefficient simulated	Structure slope (%)	Afferent area (ha)	Soil typeSCS classification	Terrain slop (%)	Land cover	Area (ha)
E1	Circle	900 mm Diameter	Concrete	0.0013	0.05	2.25	B	36.22	Impervious	0.08
							A	82.12	Forest	0.97
							B	47.22	Forest	1.19
E2	Rectangular	1.0 m height	Concrete	0.0013	0.05	29.18	B	8.96	Impervious	0.003
		1.0 m width					A	58.66	Pasture	22.73
		110.3 m length					B	47.7	Pasture	6.45
E3	Rectangular	2.0 m height3.0 m width90.6 m length	Concrete	0.0013	0.02	105.38	A	49.81	Impervious	0.01
							B	30.38	Impervious	0.03
							A	57.02	Crops	0.228
							vA	56.71	Forest	105.05
							B	28.1	Forest	0.058

The catchments contributing to the hydraulic structures of the study are characterized by small areas and steeply sloping ground (35–50%). The times of concentration were estimated using the Kirpich equation. The length of the main current and its slope are used in the equation, and the results show that the basins respond suddenly, i.e., E1 (2.3 min), E2 (7.4 min), and E3 (9.1 min).

Simulation of land-use change scenarios

The vegetation cover maps corresponding to the years 2006, 2013, and 2018 were obtained to simulate the rainfall-runoff process behavior considering the temporal variation of land use in the study area, useful to project land-use evolution for the years 2040, 2070, and 2100 through linear regression, correlating areas and land uses in each period. For this purpose, the following restrictions were considered: i) no use may overlap the motorway/highway polygon (Kalkhajeh, 2019); ii) polygons with urban planning may not evolve in areas with slopes greater than 10% (Gáfaro, 2015); iii) urban planning may not evolve in areas less than 30 m from surface currents. Likewise, elements of the land use plan of the city of Ibagué (Colombia) were incorporated relating concepts of urban growth, soil type, coverage, topography, channels or watercourses, floods, and landslides, infrastructure, and social and productive conditions in the study area, aspects that other authors have incorporated in similar studies (e.g., Romero et al., 2020; Brown et al., 2012).

Simulation of precipitation and climate change scenarios

The analysis of the effect of climate change on the drainage capacity of hydraulic structures considered in this study includes the generation of precipitation scenarios for the 2011–2040, 2041–2070, and 2071–2100 periods. These scenarios are generated by simulating precipitation on stochastic autoregressive moving average (ARMA) models to generate synthetic series through the conceptualization Escalante and Reyes (2002) proposed. In this sense, IDEAM et al. (2015) carried out the downscaling analysis of the RCP scenarios (2.6, 4.5, 6.0, and 8.5) contained in the IPCC AR5 report (IPCC, 2014) to generate more representative scenarios for the regions in Colombia. Hence, it was possible to incorporate a higher spatial resolution of precipitation anomalies in the study area.

Thus, 15 hydrological simulation scenarios that included five land-use change scenarios were analyzed by implementing 12 stationary and non-stationary models with three covariates associated with land-use evolution and 15 covariates related to climate change and weather. These analyses are summarized in Table 2.

Table 2

Combination of non-stationary models.

Run off simulation		NSA
Land use change mapLUC	Climate change scenarioCC	Non-stationary modelNSM	Code	Probability density functionPDF	Code	Soil covariate	Climate covariateCCc
2006	2011–2040	Mu (K)_ Sigma (Cov)	1	LogNormal	LN	CN	Multivariate ENSO Index -MEI-
2013	2041–2070	Mu (K)_ Sigma (Cov NL1)	2	Gumbel	GB	HU	NIÑO zone 1 + 2
2018	2071–2100	Mu (Cov)_Sigma(K)	3	Gamma	GM	KS	NIÑO zone 3.4
2040		Mu (Cov)_Sigma(Cov)	4			t_CN	Southern Oscillation Index –SOI–
2070		Mu (Cov)_Sigma(CovNL,1)	5			t_Hu	MEI_NIÑO 1 + 2
		Mu (CovNL)_Sigma(K)	6			t_Ks	MEI_NIÑO3.4
		Mu (CovNL)_Sigma(Cov)	7			CN_Ks	MEI_SOI
		Mu (CovNL,1)_Sigma(CovNL,2)	8			CN_Hu	NIÑO1+2_SOI
		Mu (CovNL,2)_Sigma(CovNL,2)	9			Hu_Ks	NIÑO3.4_SOI
		K (Stationary)	10			t_CN_Hu	MEI_NIÑO1+2_SOI
		Mu (CovNL,2)_Sigma(K)	11			CN_Hu_Ks	NIÑO1+2_NIÑO3.4
		Mu (K)_ Sigma (CovNL,2)	12			t_Hu_Ks	MEI_NIÑO1+2_NIÑO3.4
						t_CN_Ks	MEI_NIÑO3.4_SOI
						t_CN_Hu_Ks	SOI_NIÑO1+2_NIÑO3.4
						t(time)	MEI_NIÑO1+2_NIÑO3.4_SOI

Mu and Sigma represent the location and shape parameters in the PDF; K indicates that the parameter remains constant; Cov indicates the relationship of the parameter with soil or climate covariates; NL1 and NL2 represent the spline cubic non-linearity function in degree 1 and NL2 in grade 2, respectively; LN, GB and GM are the distribution functions applied.

Modeling of the rainfall-runoff process

In the study area, the meteorological station 21215180 has records for the 2010–2017 period at a 10-min scale. A typical storm was identified to have duration of 4 h after analyzing the historical records of storms. The analysis of the histograms allowed defining the duration and shape of the rainfall design used in this study.

The surface runoff produced in the area afferent to each drainage structure was estimated by applying the conceptual model Storm Water Management Model (SWMM), widely used in drainage system modeling (e.g., Babaei et al., 2018, Birgani and Yazdandoost, 2014), land-use evolution (e.g., Hongxiang and Findlay, 2013; Bhaduri et al., 2000), and climate change scenarios (Hung et al., 2020). Fifteen simulation scenarios were established to evaluate the evolution of land use and increases in precipitation in climate change scenarios in the rainfall-runoff process (Table 2). In this sense, temperature variation was not included in the hydrological modeling because the effect of evapotranspiration on storm-scale runoff production is minimal (Williams and Allman, 1969).

In this sense, the hydrological modeling was carried out by events considering the maximum rainfall of each series per year, land-use changes by periods, and climate change scenarios.

Stationary flood frequency analysis

Trends in precipitation series were evaluated with the Mann-Kendall test that has been widely used in hydrology (Zhang et al., 2008; Mwangi et al., 2016). This part of the study aims is to assess the presence of linear trends on time series data related to maximum annual daily rainfall. The significance level used in the study is alpha = 0.05, and the method selected to remove the possible presence of serial correlation is pre-whitening. This procedure assumes that the underlying mechanism generating the observation conforms to a first-order autoregressive process.

On the other hand, the maximum flow series was obtained as a result of the hydrological simulation, and the peak flow frequency analysis considered the application of the LogNormal and Gamma functions of two parameters that have been applied in flood studies (e.g., Swetapadma and Ojha, 2020; Nagy et al., 2017).

Non-stationary flood frequency analysis (NSFFA)

The non-stationary analysis was carried out by implementing the Generalized Additive Models for Location, Scale and Shape (GAMLSS) (Stasinopoulos and Rigby, 2007), applied to evaluate the effect of climate change on hydrological systems (López and Francés, 2013). Based on the results of the simulated hydrological scenarios (Table 2), the non-stationary flood frequency analysis was performed considering 12 models in combination with the LogNormal, Gumbel, and Gamma functions, covariates such as time (t), hydraulic soil properties such as curve number (CN), saturated hydraulic conductivity (Ks), soil water storage capacity in the root zone (Hu) used to evaluate the effect of land-use change on surface flows (e.g., Márquez et al., 2021; Siswanto and Francés, 2019), and covariates related to climate variability such as the Multivariate ENSO Index (MEI), NIÑO1+2, NIÑO3.4, and the Southern Oscillation Index (SOI), considered to assess effects on surface runoff in basins (Pasquini and Depetris, 2007). In this case, 16,200 simulations were carried out in the combinations presented in Table 2.

The soil water storage capacity in the root zone (Hu) represents the useful water content in the soil plus the surface storage (mm) (Eq. (1)).where (ρb) is the apparent density of the soil profile (g.cm−ɜ), (ρw) is the water density (g.cm−ɜ), (p) indicates the minimum between soil thickness and effective root depth (m), (AW) is the water available to the plants in the soil profile in percentage, (As) represents the storage of surface water retained by the effect of soil roughness (mm), and (Imax) is the maximum interception capacity in the vegetation cover canopy.

The overall methodological process is shown in Figure 2.

Figure 2

Methodological process followed.

Results and discussion

Land-use evolution

Land-use evolution was simulated for the periods 2018–2040, 2041–2070, and 2071–2100, based on land-use maps for 2006, 2013, and 2018, available for the study area (Figure 3). Figure 1 shows the values obtained for hydraulic soil properties. In the 1971–2006 period, the area related to E1 was predominantly forest, and later in the 2006–2018 period, conversion to impervious zones was observed in 4.00% of the area. Therefore, the land-use evolution simulation increased the impervious area by up to 21.77% in the year 2100. Therefore, land-use dynamics produce changes in hydraulic soil properties, in this case, manifested in a decrease of the infiltration capacity (Figures 1 and 3).

Figure 3

Land-use evolution in the catchment basins of the drainage structures under assessment for the 2006–2100 period.

On the other hand, in the 1971–2006 period, the area related to E2 corresponds exclusively to grassland, and in the 2006–2018 period, conversion to forests and impervious areas was observed in 50.96% and 0.10%, respectively. Hence, the land-use evolution for the year 2100 indicates a growth of impervious zones of 0.75% and a decrease of grasslands of 48.29% (Figures 1 and 3).

The area afferent to E3 shows a decrease in forest cover of 82.73% between 2006-2018 so that by 2100, it is expected that 0.23% of the total area will have forest cover given the growth observed in grassland and crops areas between 2006-2018, which was 5.84% and 11.52%, respectively. Thus, by 2100, a reduction in the forest area and an increase in cultivated areas, grassland, and impervious zones are expected (Figure 3), related to a decrease in infiltration capacity (Figure 1). Similar land-use evolution rates have been reported in simulations based on 10, 15, and 29-year-periods (Liu et al., 2017; Romero et al., 2020).

Stochastic simulation of precipitation

Stochastic ARMA-type models were implemented from the historical precipitation series to generate synthetic series. The model type selection was based on the Akaike information criterion (AIC), resulting in the best ARMA model (1,1). The models were applied in the weather station 21245040 for the 1971–2100 period, where no positive trends were detected in the historical time series, and in the weather station 2121060 for the 1984–2100 period that exhibited trends. The evaluation of trends in precipitation series was carried out by applying the Mann-Kendall test.

Hydrological response in land-use change scenarios

Land-use evolution has been reported to produce effects on surface runoff (Hu et al., 2021; Peña et al., 2016). This study analyzed the effect of land-use evolution on the peak flow and capacity of structures E1, E2, and E3.

In the case of E1, there is an increase of 8.64% in maximum flow during the 2006–2100 period, related to a decrease in hydraulic conductivity (Ks) of 18.93%; the capacity of E1 experienced a decrease from 98.00% to 87.71%. Likewise, structure E2 showed a decrease in peak flow of 4.62% during the 2006–2100 period, explained as a consequence of an increase in infiltration capacity from 8.05 to 49.54 mm/h. This increase is due to an expansion in forest cover in the runoff contribution area to E2, which increased the capacity by 12.53% in the structure (Table 3 and Figure 4).

Table 3

Structure-weighted soil characteristics and maximum flow per structure for land-use change scenarios.

Structure	Parameter	2006	2013	2018	2040	2070	2100
E1	CN	57.58	57.58	57.6	58.71	60.59	62.56
	Ks (mm h-1)	60.79	60.79	60.77	58.14	53.82	49.28
	Hu (mm)	127.56	127.56	127.41	122.18	113.1	103.56
	Q (m3s-1) (Station 2121060)	1.458	1.455	1.453	1.487	1.535	1.584
	Q (m3s-1) (Station 21245040)	1.173	1.171	1.169	1.201	1.247	1.326
E2	CN	70.43	60.09	60.09	60.29	60.31	60.33
	Ks (mm h-1)	8.05	49.68	49.68	49.56	49.55	49.54
	Hu (mm)	13.65	76.01	76.01	76.09	76.07	76.04
	Q (m3s-1) (Station 2121060)	21.171	20.371	20.348	20.208	20.211	20.193
	Q (m3s-1) (Station 21245040)	16.407	16.310	16.210	19.160	16.165	19.060
E3	CN	45.06	48.22	48.22	52.06	57.57	63.35
	Ks (mm h-1)	114.4	97.81	97.81	77.71	48.71	18.26
	Hu (mm)	139.79	119.68	119.68	93.77	60.12	23.19
	Q (m3s-1) (Station 2121060)	60.008	65.761	65.761	71.656	77.492	80.276
	Q (m3s-1) (Station 21245040)	45.868	51.172	51.172	56.498	61.698	63.900

Figure 4

Effect of land-use scenario in weather station 2121060 on a) Discharge in E1, b) Discharge in E2, c) Discharge in E3, d) Capacity in E1, e) Capacity in E2, and f) Capacity in E3.

On the other hand, the contribution area to E3 showed a forest cover decrease of 99.77% and an increase in areas with grassland by 33.72%. Therefore, the infiltration capacity (Ks) decreased by 15.96%, related to the increase in peak flow of 33.77% during the 2006–2100 period (Table 3 and Figure 4).

Similar behavior has been reported in observations of paired basins in transition from grassland to forests, where there has been a decrease in the runoff volume in a range of 48–78% (Marshall et al., 2014), as well as an increase of 5–30 times the volume of surface runoff in transition from forested to impervious areas (Hurni et al., 2005). It has also been pointed out in the literature that the evaluation of the effect of land use on surface runoff through hydrological modeling indicates that when there is an increase in areas with forest, a decrease in peak flow is experienced in ranges from 6 to 14% (Kabeja et al., 2020; Ruman et al., 2021). Accordingly, Kalantari et al. (2014) found that the change from forests to grassland in a period of 50 years generated an increase of 5% in the water layer height that produced a decrease in the hydraulic capacity in road drainage structures.

Non-stationary modeling of the maximum flow regime in land-use change and climate change scenarios

This section presents the results obtained from modeling the maximum flow regime in stationary and non-stationary contexts in land-use change scenarios considering covariates such as soil water storage capacity in the root zone (Hu), hydraulic conductivity (Ks), and the infiltration capacity represented as the curve number (CN). In the current research, the maximum flows simulated in the 1976–2100 period were evaluated in the scenarios described in Table 2 based on the precipitation series of the weather station 2121060, which shows an increasing trend, and the weather station 21245040 with no trend detected.

The results indicate that the time series of maximum flows generated in the areas afferent to structures E1, E2, and E3 in land-use change scenarios show non-stationarity; the Gamma distribution function is the one that presented the best fit in the analysis related to station 2121060 and LogNormal concerning station 21245040. Another aspect to highlight is that the applied models that best describe the variation in the magnitude and frequency of peak flows for E1 correspond to GM4 and LN4. These involve t and t_Ks, respectively, as a combination of covariates. For E2, the models are GM7 and LN8 that include Hu_Ks and CN_Hu as covariates. For E3, the models are GM4 and LN7 that adopt the covariates t and CN_Hu (Table 4 and Figure 5).

Table 4

Summary for the best fit models indicating the selected distribution, significant covariates, model code, and the Akaike Information Criterion (AIC).

Weather Station	2121060		21245040		2121060		21245040		2121060		21245040
Land-use change covariates	AIC	Code	AIC	Code	AIC	Code	AIC	Code	AIC	Code	AIC	Code
Stationary	46.180	GM	-60.139	LN	649.38	GM	592.102	LN	1015.337	GM	996.915	LN
t	23.053	GM4	-116.329	LN7	637.094	GM4	566.106	LN9	920.279	GM4	896.048	LN9
CN	24.460	GM4	-112.054	LN4	634.628	GM7	565.92	LN4	932.648	GM5	900.915	LN9
HU	24.437	GM4	-108.388	LN9	634.763	GM7	565.676	LN4	932.968	GM5	899.132	LN9
KS	24.493	GM4	-112.08	LN4	634.739	GM7	565.712	LN4	932.628	GM5	900.914	LN9
t_CN	26.609	GM4	-116.795	LN4	635.888	GM4	566.585	LN4	920.602	GM5	895.649	LN9
t_Hu	26.601	GM4	-114.588	LN7	635.976	GM4	566.585	LN4	922.714	GM4	907.785	LN4
t_Ks	26.622	GM4	-116.804	LN4	635.957	GM4	566.587	LN4	922.634	GM4	908.604	LN4
CN_Ks	24.383	GM5	-112.263	LN8	634.568	GM7	568.003	LN4	932.688	GM5	899.851	LN4
CN_Hu	26.864	GM9	-109.223	LN4	638.926	GM8	565.63	LN8	930.097	GM5	892.05	LN7
Hu_Ks	26.894	GM9	-109.25	LN4	634.392	GM7	568.39	LN4	930.129	GM5	899.887	LN9
t_CN_Hu	27.532	GM5	-115.506	LN7	635.505	GM4	567.394	LN9	922.137	GM4	905.529	LN4
CN_Hu_Ks	26.430	GM5	-108.344	LN9	636.390	GM7	571.634	LN9	932.096	GM5	899.789	LN9
t_Hu_Ks	27.406	GM5	-113.321	LN7	635.131	GM4	567.888	LN9	926.284	GM4	905.998	LN4
t_CN_Ks	25.298	GM5	-116.108	LN9	635.627	GM4	566.964	LN9	920.855	GM5	899.193	LN4
t_CN_Hu_Ks	24.607	GM5	-113.210	LN9	638.027	GM9	567.835	LN9	923.340	GM5	903.613	LN4
Climate change covariates	AIC	Code	AIC	Code	AIC	Code	AIC	Code	AIC	Code	AIC	Code
Stationary	10.936	GM	-32.4	LN	155.615	GM	220.28	LN	213.032	GM	323.208	LN
t	10.936	GM10	-35.187	LN4	155.615	GM10	215.272	LN4	213.032	GM10	315.072	LN4
MEI	-15.929	GM8	-31.155	LN7	155.219	GM1	221.397	LN1	212.851	GM1	324.404	LN1
NIÑO1+2	4.764	GB8	-39.328	LN8	156.595	GM3	217.801	LN5	214.036	GM11	319.059	LN5
NIÑO3.4	5.121	GM8	-30.611	LN1	155.368	GM1	221.928	LN1	213.035	GM1	324.895	LN3
SOI 5	-6.055	GM8	-29.481	LN12	150.298	GM2	222.203	LN1	213.545	GM3	325.009	LN1
MEI_NIÑO 1+2	7.582	GM8	-34.38	LN5	156.258	GB6	218.879	LN5	214.086	GM1	322.964	LN3
MEI_NIÑO3.4	-2.140	GM8	-28.336	LN12	157.29	GM12	224.47	LN12	214.767	GM1	326.005	LN11
MEI_SOI	-6.501	GB8	-30.392	LN6	156.284	GM1	223.259	LN1	214.849	GM1	324.615	LN6
NIÑO1+2_SOI	8.067	LN6	-31.443	LN11	154.56	GM12	222.634	LN6	211.024	LN6	325.237	LN4
NIÑO3.4_SOI	-2.561	GB8	-35.38	LN8	155.242	GM5	218.478	LN6	214.512	GM12	319.141	LN8
MEI_NIÑO1+2_SOI	10.936	GM10	-31.97	LN3	157.175	GM1	221.401	LN11	214.144	GM7	325.124	LN3
NIÑO1+2_NIÑO3.4	-2.174	GB8	-34.535	LN5	156.178	GM1	218.141	LN5	213.701	GM7	318.097	LN5
MEI_NIÑO1+2_NIÑO3.4	6.657	GM9	-32.706	LN8	157.961	GM1	220.739	LN8	215.572	GM11	321.282	GM2
MEI_NIÑO3.4_SOI	10.936	GM10	-30.515	LN11	157.712	GM1	223.979	LN2	215.456	GM12	3,25,119	LN11
SOI_NIÑO1+2_NIÑO3.4	9.389	GM6	-31.842	LN6	156.068	GM1	220.414	LN11	212.848	GM6	324.782	LN3
MEI_NIÑO1+2_NIÑO3.4_SOI	6.161	LN8	-30.446	LN11	154.908	GM8	224.814	LN12	210.895	GM8	325.763	LN12

Underlined and bold numbers and codes are the best values obtained.

Figure 5

Summary of the results of modeling annual maximum peak discharge in the three representative structures used in the study with models under stationary and non-stationary conditions. The results show the median estimates (dark line) and the 5th, 25th, 75th, 87.5th, 95th, and 99th percentiles.

On the other hand, the Kendall, Spearman, and Pearson correlation tests were applied to define the covariates best related to the precipitation behavior in the study area and observe the degree of dependence. The climatic indices used as covariates were the Multivariate ENSO Index (MEI), NIÑO1+2, NIÑO3.4, and the Southern Oscillation Index (SOI). The implementation of the different statistical models proposed to model the maximum flows regime generated in the areas related to structures E1, E2, and E3 showed a strong dependence of the parameters of the distributions relating to the climatic indices. The results indicate the need to address the analysis of frequencies in the maximum flow regime in a non-stationary context, incorporating the effects of changes in climate. The models that best describe the variation of peak flow magnitude and frequency for stations 2121060 and 21245040 in the case of E1 correspond to GM8 and LN8, which incorporate the combination of covariates MEI and NIÑO1+2, respectively. For E2, models GM2 and LN4 that include SOI and t as covariates, respectively; and for the E3 structure, models GM8 and LN4, which adopt MEI,_NIÑO1+2,_NIÑO3.4, and SOI, and t as covariates (Table 4 and Figure 5). In this sense, incorporating climatic covariates can improve the magnitude and frequency variability description of the maximum flow regime, as reported by Nasri et al. (2017) and López and Francés (2013). Likewise, the dependence of the distribution function parameters with respect to time was identified when working with the precipitation series that experienced a trend (station 2121060). This result shows the need to incorporate models that assume gradual changes in the hydrological series, as previously reported by Obeysekera and Salas (2016).

The models that involve two covariates showed the best adjustments in 66.67% of the cases when performed in analyses with covariates associated with land-use change. In contrast, the use of covariates related to climate change indicated that 83.33% of the models applied presented the best adjustments when a single covariate was used (Table 4 and Figure 5).

The goodness of fit evaluation of the selected models was carried out by verifying the normality and independence of residuals. The process consisted in verifying the first four statistical moments of the residuals (mean, variance, skewness, and kurtosis) and the Filliben correlation coefficient (Table 5). A visual inspection of quantile graphs without a “worm plot” trend was also carried out where the shape of the chart allows inspecting visually how the data differs from the arranged distribution, a normal one in this case (Figure 6). This analysis is approached to verify that the selected models can adequately describe the systematic part and compare the models using land-use and climate change covariates. In this context, the analysis of residuals in the current study indicates that when covariates related to climate change are incorporated, the description of the variation of the maximum flow regime generated in the area afferent to each structure analyzed will improve in comparison with the results obtained when using covariates associated with land-use change (Figure 6 and Table 5).

Table 5

Normality of residuals for annual maximum land-use change and climate change scenarios.

Weather station	Structure	Non-stationary Model	Mean	Variance	Skewness	Kurtosis	Filliben correlation
2121060	E1	t_GM4	0.001	1.009	-0.376	2.978	0.994
	E2	Hu_Ks_GM7	0.001	1.009	-0.350	2.578	0.989
	E3	t_GM4	-0.004	1.005	-0.223	2.767	0.996
	E1	MEI_GM8	2.774 e-13	1.039	-0.086	2.591	0.991
	E2	SOI_GM2	-0.117	0.965	-0.128	2.038	0.984
	E3	MEI_NIÑO1+2_NIÑO3.4_SOI_GM8	-0.002	1.035	-0.173	2.298	0.991
21245040	E1	t_Ks_LN4	-0.006	1.008	0.166	2.819	0.998
	E2	CN_Hu_LN8	1.84E-05	1.008	0.198	2.824	0.994
	E3	CN_Hu_LN7	-0.022	1.007	0.135	3.485	0.991
	E1	NIÑO 1+2 _LN8	-0.024	1.021	0.371	2.148	0.983
	E2	t_LN4	0.029	1.021	0.220	2.122	0.988
	E3	t_LN4	0.021	1.021	-0.089	2.280	0.988

Figure 6

QQ-plot residuals for non-stationary models with land-use change and climate change covariates.

The NSFFA of maximum flows in each road drainage structure indicates that when considering the hydrological series of station 2121060, the best adjustments are associated with covariates related to climate change. Particularly for the E1 structure, the model that best describes the maximum flows in climate change scenarios with the precipitation series for station 2121060 corresponds to GM8, which incorporates the behavior of MEI as a covariate. In contrast, in land-use change scenarios, the best fit is obtained with the LN4 model that considers infiltration capacity (Ks) and time (t) as covariates with the precipitation series for station 21245040. Likewise, the structure E2 with station 2121060 shows better results using covariate SOI, and the GM2 and LN8 models present the best adjustments with the CN_Hu covariates for station 21245040. Finally, for the E3 structure with station 2121060, the results improve with the GM8 model using MEI_NIÑO1+2_NIÑO3.4_SOI covariates; and with respect to land use, the LN7 model with covariates CN_Hu for station 21245040.

In this case, the results of the NSFFA indicate that the best adjustments are found when covariates associated with climate change are applied in the precipitation series where trends are detected, as in the case of station 2121060 (Table 4 and Figure 7). Likewise, when the NSFFA was carried out based on the precipitation series of station 21245040 that does not show trends, the best results are related to incorporating covariates associated with land-use change (Table 4 and Figures 5, 6, and 7).

Figure 7

Comparison of different models with soil and climate covariates for E1; a) Covariates of land-use change station 2121060 model GM4; b) Covariates of climate change station 2121060 model GM8; c) Covariates of land-use change station 21245040 model LN4; d) Covariates of climate change station 21245040 model LN8.

On the other hand, the dependence of soil water storage capacity in the root zone (Hu) and the infiltration rate of the vertical flow (Ks) with respect to the magnitude of the maximum flows generated in the areas afferent to the three drainage structures are identified. These covariates improve the NSFFA related to the effect of land-use change on runoff production and, in turn, to the design of hydraulic structures (Chen et al., 2017; Wijesekara et al., 2010). Likewise, variations in climate and soil generate non-stationarity effects on the maximum flow regime in the three afferent areas analyzed, indicating that the flows obtained through stationary models are not adequate for the design of the hydraulic infrastructure agreeing with the work of Cheng et al. (2014). This could indicate the need to incorporate a non-stationary analysis through the use of covariates related to land-use change (Siswanto and Francés, 2019) or climatic indices (López and Francés, 2013), or a combination of the two, to obtain better adjustments in the design of hydraulic structures (Kalantari et al., 2014).

The best results of the NSFFA are obtained using climate change covariates such as NIÑO 1 + 2, and for land-use change, the behavior of Hu. All cases are due to particular conditions, such as the structure type, basin shape, and hydroclimatology, explaining that there is no specific covariate or a model that optimally fits all the design conditions of hydraulic structures (Serago and Vogel, 2018; Bertoni, 2010).

Conclusions

To better understand the frequency analysis of peak flows in land-use change and climate change scenarios and try to establish a framework to address the non-stationary flood frequency analysis for the dimensioning of road drainage structures, this study performs a comparison of the performance of statistical models in the following contexts: i) stationary, ii) non-stationary (land-use changes) and iii) non-stationary (climate change), through statistical simulation.

Comparing the different typologies of proposed statistical models showed that non-stationary models could be more efficient in representing the variation of the maximum flow regime when non-stationarity is detected in time series. In this sense, the presence of non-stationarity in time series can increase accuracy deterioration when the stationary frequency analysis is approached.

The results indicate that the change in land use produces a non-stationary effect on the production of runoff in the areas afferent to each drainage structure according to the magnitude and evolution of the type of soil cover. Therefore, incorporating covariates related to hydraulic soil properties contributes to improving the descriptive capacity of non-stationary models in the frequency analysis of maximum flows. Mainly, the combination of two covariates shows the best results, at least in the context of this study. Similarly, combining two covariates associated with climate change can obtain excellent approximations to the behavior and variation of the maximum flows in the analyzed drainage areas.

In this study, the use of covariates related to land-use change describes the best adjustments in the non-stationary frequency analysis when the hydrological response of the basin is related to precipitation series in which no trends are detected. Likewise, the use of covariates related to climate change contributed to obtaining the best descriptive capacity in the analysis of the frequency of maximum non-stationary flows when the hydrological response of the basin is related to a precipitation series in which a growing trend was detected. Therefore, other studies should include basin-scale tests using hydrological modeling or observations in paired basins, which can help confirm or refute these findings.

The results of this study contribute to generating a framework that allows establishing the selection between a stationary model and a non-stationary model, which is urgent for its implementation in hydrological praxis. However, the understanding of the problem is incomplete since more in-depth research is necessary to select the best models in the context of accuracy and the uncertainty associated with the models. The non-stationary frequency analysis should be promoted when diagnosing the presence of non-stationarity in the hydrological time series, but this study demonstrates the need to identify a suitable model structure for its practical application.

Declarations

Author contribution statement

Mónica Jiménez-U: Performed the experiments; Contributed reagents, materials, analysis tools or data; Analyzed and interpreted the data; Wrote the paper.

Luis E. Peña & Jesús López: Conceived and designed the experiments; Analyzed and interpreted the data; Wrote the paper.

Funding statement

This work was supported by the (12-262-COL00) and the Universidad de Colima (PRODEP-UCOL-PTC-210).

Data availability statement

Data included in article/supplementary material/referenced in article.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.