Risk Zone Prediction in Meandering Rivers by Using a Multivariate Approach

River meandering prediction is one of the most complex physical phenomena regarding river hydraulics. The constant cause—effect interaction between the propagation of a flood wave and the geologic conditions of a cross section of a river are, with no doubt, the genesis of river meandering. As a consequence of this interaction, certain processes like lateral migration (owing to margins erosion), river widening, aggradations and degradation of the river bed, land forms of the bottom, evolution and variation on suspended sediment concentration, are fluvial processes that occur along and across the functional components of a river, transmitting with this, a unique characteristic to the flood plains and the mouth of the river. All of these fluvial processes are known, generically, as river migration. At the present, based on this concept, a number of studies are being carried out to predict river behavior. This is important, but also the hydrodynamic characteristics of the river impose a greater degree of complexity, especially when such models are used to predict meandering river migration (Abad and Garcia 2006).

General studies on the patterns that a river can follow along its meandering course can be found in the works of Hickin (1974), Hooke (1984, 1995) and Hooke and Redmond (1992). Several methodologies to characterize a meandering river have also been proposed (Hooke 1984; Oneill and Abrahams 1986; Howard and Hemberger 1991; MacDonald et al. 1992), however, only a few computer software packages have been developed specifically to describe and predict river migration (Lagasse et al. 2004; Abad and Garcia 2006). In Latin America, RVR Meander model (Abad and Garcia 2004) constitutes the most precise tool for meandering prediction. Three modules, one of preprocessing, a statistical one and the river migration tool itself, constitutes this software. It is a model that focuses on the morphologic processes at a medium temporal scale. Its application has been successful in describing the evolution of Bermejo River in Argentina because of the influence of the construction of the Lavalle Bridge. The dominium of the study zone is around 100 km, and the results were validated with satellite images for different time steps.

A procedure to characterize river meandering must consider a complete description of the fluvial cycle. That is, a characterization of the youth, mature and full development stages of the river (Mueller 1968). This procedure is based on the identification of the two types of sinuosity: topographic and hydraulic. In the latter, the form of the meander is associated to the delimitation of the fluvial corridor and to the configuration of the geologic platform of the stream.

For the case of topographic sinuosity, the so-called Kinoshita curves are used. These curves are based on an empirical formulation developed by Beck (1988). The formulation of Kinoshita curves takes into account a number of physiographic characteristics the configuration of the river; however, it is the curve amplitude, which represents the main characteristic of this formulation (Parker et al. 1983; Parker and Andrews 1986). This main characteristic is known as the angular sinuosity coefficient ($\theta$). In this paper, this formulation is used, and an alternative expression for meander prediction, based on a stochastic multivariate analysis of the geomorphologic and physiographic characteristics of a river, is proposed. In this sense, the main objective of this paper is to show a methodology for the establishment of a risk zone based in the new meander prediction methodology. The method has been validated with stochastic models by using 3480 simulated occurrences of the sixteen characteristics proposed to characterize the meandering of Cahuacán River in the Mexican state of Chiapas.

The Mexican state of Chiapas is located at the southeast of the country (Fig. 1). Its limits are to the north, the state of Tabasco; to the south, the Pacific Ocean; to the east, Guatemala and to the west, the states of Oaxaca and Veracruz. The state has a total surface of $\mathrm{73,724}{\mathrm{km}}^2$, which represents 3.8% of the national territory. It is constituted by three hydrologic regions (Table 1): Coast of Chiapas, Grijalva—Usumacinta and Coatzacoalcos; although the Coatzacoalcos region is merely symbolic because it reaches only 0.03% of the state surface.

The particular study zone is located in the Coast of Chiapas. This region covers the south coast of the state running along the Pacific Ocean. It is composed of a number of rather small basins such as Suchiate, Coatán, Huixtla, Cacaluta and Novillero among others. The first two has their beginning in Guatemalan territory; Huixtla basin comprises the Cintalapa River and Los Cerritos water body. The Pijijiapan basin is also located in the region, with its main stream, Pijijiapan River and La Joya and Buenavista water bodies. There also exist several small rivers and at the end, a coastal lagoon called Mar Muerto. The region has $\mathrm{11,748}{\mathrm{km}}^2$, 10,566 of which belong to Chiapas State and 1,182 to Oaxaca State. In this study, only de Chiapas part was considered.

The study zone is located on the Cahuacan River, specifically between the municipalities of Tapachula and Tuxtla Chico (Fig. 2), covering almost 2.5 km along its course over the coastal floodplain near to Guatemala border, on a zone of sedimentary deposits and basic extrusive igneous rocks.

The conformation of the river morphology in space and time is a phenomenon that involves a number of physiographic and dynamic variables that in most of the cases, are difficult to determine. The group of variables and morphologic components of a river are unique for each case, and the number of variable that influence its behavior is so high that their determination is almost impossible. This is because many of them are unknown or have not been estimated as dominant relationships in the meandering process of the river. A procedure, proposed by Abad and Garcia (2006), consists in identifying three main stages. A preprocessing, which consist in obtaining the exact location of the $n$ points, that all together, conform the trace of the river along its axis. This stage includes for example, the digitalization of the river, the use of aerial photography and satellite imagery that can be incorporated to Geographic Information Systems. In several cases, it can be necessary to screen and correct those points that have an infield provable deviation. This stage of preprocessing is used in the software RVR Meander, to represent the horizontal coordinates of the river through parameterized third-order polynomials, which allow the calculation of the first and second derivative of the function in a continuous form. Once the first stage is fulfilled, the configuration of the river is completed. Fig. 3 shows the graphical representation of a 230-km reach of the Cahuacan River.

Second stage involves a statistical processing. In this step, the evolution of the river is appropriately estimated. Parameter such as the lateral and longitudinal displacements of the river, the amplitude, the sinuosity, and the mean curvature (MacDonald et al. 1992) are estimated. Additionally, some other parameters, corresponding to dynamical aspects such as the washed load and the dominant discharge, are used. Table 2 and Fig. 4, show the morphologic and fluvial characteristics for the analysis of the river migration by means of a stochastic sinusoidal scheme. Tables 3 and 4 show the values of such characteristics for the Cahuacan River. Fig. 4 shows each one of such parameters.

The third stage is properly the prediction of the river migration. In this step it is important to define the formulation of the curves configuration: symmetrical (Langbein and Leopold 1966) or asymmetrical trends (Kinoshita 1961; Kinoshita and Miwa 1974; Parker et al. 1983; Yamaoka and Hasegawa 1984; Parker and Andrews 1986). In the first case, it is recommended to adjust to a sinusoidal function. For example, in the RVR Meander model, in this stage, lateral migration is predicted. In this model, the hydraulic hypothesis of quasi-steady flow is assumed, which is valid only for erodible rivers with a constant bed width at all times of the simulation. To model erosion of vertical walls, the concept of excess velocity is employed. In these models, some simplifications are often used, such as the no inclusion of the continuity equation for sediments or the consideration of a lineal profile in the bed of the river (Abad and Garcia 2006). The same considerations are assumed in this paper, however, in difference of the RVR Meander model, which uses the Kinoshita curve as the main formulation, a new trigonometric expression to generate curves is proposed. The new procedure does use the form of the Kinoshita curve, but it affects a succession of sines and cosines directly with the morphologic characteristics of the river.

With this basic idea, it is usual to define the mean variation downstream the valley ($z$) and the width of the cross section of the same valley ($y$) (MacDonald 1991; MacDonald et al. 1992) aswhere $n$ = component, perpendicular to the river axis, along the river; $s$ = component, perpendicular to $n$, that forms an angle $\theta$ with the reference horizontal plane; $\mathrm{\Delta }t$ = time step of occurrence.

In this paper, a new equation with a succession of sines and cosines is proposed. The amplitude of each trigonometric function is formed by the values of each one of the morphologic and characteristics previously introduced ($X_i$). Furthermore, the frequency of sines and cosines is the product of each one of the principal components $|{\mathrm{\Omega }}_i|$ multiplied by the development longitudinal distance of the meander. This means that the width of the cross section of the valley ($y$) is a function of the meandering variation along the valley ($z$) and the morphologic and fluvial characteristicswhere $X_i^{″}$ are the proposed central standardized values of the morphologic and fluvial characteristics $X_i^{″}=(X_i-{\mu }_x)/{\sigma }_x$; $|{\mathrm{\Omega }}_i|$ $i$-th principal component.

A similar expression is the one proposed by Kinoshita to generate curves in meanders (Beck 1988)where $J_s$ = skewness coefficient; $J_f$ = kurtosis coefficient; ${\theta }_0$ = maximum angle of amplitude; $\lambda$ = arc of curvature of the channel; $s$ = coordinate of rotation.

Beck (1988) introduced a qualitative comparison between the mean sinuosity and the riverbed elevation with the symmetric curvature. Hooke (1974) presented these relationships by means of three equations known as the Beck relationships (Abad and Garcia 2006). The same comparison was made in rivers with high sinuosity (Hills 1987). The results suggest that there exists a component, along the river, that relates the maximum cross slope of the river to the maximum ratio of curvature. However, the other topographic components are well described. Similarly, a comparison between the Kinoshita formulation [Eq. (3)] and the proposed expression in this paper [Eq. (2)] is also presented. Fig. 5 shows the sinuosity obtained by the use of the values presented by Beck formulation and by using the same parameters, the sinuosity, through Eq. (2) for $i=1$ (a single pair of sines and cosines).

The similarity between the two formulations is quite evident, even though a lag is noticed. In the proposed method, this can be fixed by adding extra pairs of sines and cosines, while they are needed until the desired precision is reached. It can be mentioned that in both cases, the used values were those presented by Beck (1988) and reproduced with detail by (Abad and Garcia 2004, 2006), ${\theta }_0=20°\mathrm{Js}=1/32\mathrm{Jf}=1/192$ $y$ ${\theta }_0=100°\mathrm{Js}=1/32\mathrm{Jf}=1/192$.

The fundamental idea of the proposed equation is the use of a number of pairs of sines and cosines that allow the best description of the path and migration of the meander. Additionally, these trigonometric functions will be affected by the maxima of morphologic and fluvial characteristics given in Table 2, once they have been selected and prioritized through a multivariate analysis. For example, in the case under study, with the estimation of the sinuosity of the Cahuacan River, the formulation can be extended up to 15 morphologic and fluvial characteristics, which gives the result of having seven pairs of both sines and cosines ($n=7$) plus an independent parameter. This latter parameter will be the most significant variable.

The proposed formulation is based additionally on the idea that a group of variables is able to specify by themselves, a certain resemblance, that when graphed in a diagram, they will reveal similarities or grouping characteristics among them. This consideration is often employed in a systematic way in the delineation of hydrologic homogeneous regions (Everitt 1978; Wiltshire and Beran 1987; Donald 1988; Nathan and McMahon 1990).

This kind of multivariate development foresees a previous discrimination of the variables, which allow emphasizing the importance of each one of the characteristics. A classic example of this approach is the technique known as multidimensional trace (Andrews 1972) which is given bywhere $X_i$ represent each one of the analyzed and prioritized variables.

A characteristic of this expression is that its results depend on the magnitude of the selected variables. This fact is usually not noticed at once. For this reason, it is important to standardize the values of the characteristics before any multidimensional trace is done. However, the additive property of sines and cosines gives cyclic characteristics. In this form, the first characteristics are associated to cyclic components of the low frequencies and the last, to components of the high frequencies. This certainly can help to associate the variable to its importance. Low frequencies are much more difficult to observe and so, $X_1$ represents the most important variable, $X_2$ the second in importance and so on. This formulation is similar to Eq. (3) and with the help of a succession of sines and cosines, affected by the principal components [Eq. (2)], the expression to estimate river meandering is obtainedwhere $X_i^{″}$ are the central standardized morphologic and fluvial characteristics; $|{\mathrm{\Omega }}_i|$ is the $i$-th principal component; $n$ is the number of pairs of sines and cosines.

As was already mentioned, before using the proposed scheme, it is necessary to prioritize the variables. With the values of the morphologic and fluvial characteristics shown in Table 2, an EOF analysis was carried out. In this analysis, it was found that the first two principal components explain 70.5% of the variance of the morphologic and fluvial characteristics of the Chauacan River. Fig. 6 shows the circle of correlation coefficients among the 15 proposed variables.

EOF results show the existence of three groups of variables that define the behavior of the Cahuacan River (Fig. 6). The first group is formed by the morphologic characteristics of the river like the curvature ratio, the sinuosity and the lengths of curve and wave ($r/b$, $r$, ML, Lc and b2). The second group is composed by the merely hydrologic features which include the discharge and the basin surface (Qd, FA and Ac). The third group represents both morphologic and geometric features of the river (a, B, L, Gc, Ang $\theta$, $b$ and $b1$). This analysis can offer additionally a prioritization of the characteristics based on the projection of each variable over the axis of the principal components $|{\mathrm{\Omega }}_i|$. This procedure corresponds to the traditional interpretation of an EOF.

Likewise, the values of the characteristics shown in Table 2 are standardized, obtaining with their maximum and minimum.

These results are shown as highlighted in bold font of Table 5, in which it can be observed that the variable with highest importance for the first principal component $|{\mathrm{\Omega }}_1|$ is the sinuosity ($r/b$), and that of highest relevance for the second component $|{\mathrm{\Omega }}_2|$ is the meander width (B).

On the other hand, the maximum and minimum values of each one of the characteristics of the Cahuacan River are estimated. These values allow obtaining the limit conditions up to date, of the possible morphologic configuration of its meanders (Table 5). Later, the prioritized characteristics and the principal components are used to represent the river meandering in the study reach. For the case of Cahuacan River, nine features $X_{i,i=1,\dots ,9}$ ($n=4$ plus an independent term) and the first four principal components $|{\mathrm{\Omega }}_{i,i=1,\dots ,4}|$ were selected. With these values, Eq. (5) can be developed in the following form:where $\dot{z}$ is the variation of the meander conditions along the river axis (2.3 km).

By substituting the maximum and minimum morphologic characteristics of the river in the last expression, a meandering prediction as a graphical representation of the sinuosity of the river is obtained. Figs. 7 and 8 show the results of this procedure.

It is worthy of mention that the use of four pairs of sines and cosines cause that signal to distort and then white noise arises. To smooth, the trace it is necessary to take the moving average of each condition. The sinuosity thus obtained has the advantage that it can be displaced along the longitudinal axis of the river, which allows the reconfiguration of the meander, according to the analysis when comparing the proposed expression in this paper versus Beck formulation (1988).

Fig. 9 shows the comparison of the maximum and minimum conditions and the natural configuration of the Cahuacan River meandering. It is interesting to observe how, once smoothed, the trace of the behavior of these two conditions result close to the natural condition. Even though the crests and valleys are not totally reproduced, the representation is a very close approximation to the actual conditions of the river meandering.

The maximum and minimum conditions used in this sinuosity formulation give by themselves a limit condition in meandering. However, to ensure these results have infield practical utility it is particularly interesting to compute certain confidence limits involving both conditions. So, the use of the classic formulation to define such limits is proposed. By selecting a tolerance of 95%, the limits for each condition will be given by the following expressions:where ${\sigma }_{\dot{y}}{|}_{\dot{z}}$ = standard deviation of the values of $y$ along the river meander axis.

Confidence limits are built using Eq. (8) for each one of the conditions already mentioned. Figs. 10 and 11 show the results.

It is accepted that the trace of the limits comes from a formulation, which involves the morphologic and fluvial characteristics or the Cahuacan River, then the border obtained in this way could delimit the extreme conditions of the river movement. However, it is worth mentioning that in a rigorous way, the construction of these limits is based in a purely statistical condition and do not involve the physical conditions of the river. This can be solved for example, by adding to the proposed methodology, an extra stage of stochastic simulation to define directly the borders or meandering zones under risk.

The conception of idealizing and representing the evolution of the possible tendency of river meandering is very complex. Various authors have introduced several schemes of this process (Duan 1998; Darby and Delbono 2002; Jang and Shimizu 2005). Lots of them draw qualitative comments, whereas some others are completely quantitative, but in all cases, the consideration of the curvature of the meander is the key element for the correct representation of the phenomenon. In this order of ideas, this common consideration has been proved through the principal component analysis, which put the curvature ratio ($r$) and the bed width ($B$) as the two most important and representative features of the meandering process. An idealized trace, considering these morphologic characteristics, has the advantage of the exhaustive proof through statistical testing to demonstrate its validity or even to systematize the expressions in diverse numerical models (Bradshaw 1969, 1973; Chebbi et al. 1998; Blanckaert 2002; Merwade et al. 2005; Letrung et al. 2011). However, it is worth mentioning that when idealizing the trace, several considerations are omitted, affecting with this the river dynamics. For examples, the equilibrium condition on the channel width and the sediment load (Ackers and Charlon 1970; Whiting and Dietrich 1993a, b; Garcia and Niño 1993; Shields et al. 2003; Abad and Garcia 2004; Mizumura 2011).

Results obtained through the proposed formulation show an acceptable coherence with the traditional formulation (Beck 1988) which used the Kinoshita approach. Trace smoothing using a moving average process is adequate to describe a meander completely, and can be directly used in the field to define risky zones. On the other hand, the consideration of using a multivariate analysis to determine the relative importance of each morphologic and fluvial characteristic is highly recommendable and gives the methodology an extra degree of confidence in the numeric sense. Furthermore, an improvement on the proposed methodology can be reached by adding a stochastic component. There are two ways to consider this component. The first one consists of the adjustment of an autoregressive model to each one of the historical observation series (morphologic variables) of the river $y$ with these equations to generate samples or predict the future conditions of meandering. The second form is to directly add an extra term to the formulation as it is following shown:where $\epsilon (\dot{z})$ is a random number $\epsilon ={\widehat{\sigma }}_{\epsilon }{\mathrm{\Omega }}_{1,2}$; ${\mathrm{\Omega }}_1$, ${\mathrm{\Omega }}_2$ are random numbers (standard normal); ${\alpha }_1$, ${\alpha }_2$ are random numbers with uniform distribution [0,1].

By considering as example, the equation proposed by Box and Muller (1958), random number can be generated by the following expressions:

Fig. 12 shows the results of the stochastic simulation (ARMA models) of some years in the close future by taking into account three conditions (Anderson et al. 2007). The first is based on the historical conditions of the river, the second is based on the maximum conditions, and the third one is based on the minimum meandering condition of the Cahuacan River. This means that results of stochastic simulation represent what could happen by 2016.

In the same order of ideas, an analysis of moving averages to smooth the confidence limits was carried out. This representation constitutes the main result of this work because this image when put on top of an aerial photograph of the zone under study (Fig. 13), allows the exact location of regions or lots under risk because of meandering dynamics. Finally, it is worthwhile to mention that this analysis is based on 15 morphologic and fluvial features, however, the inclusion of other dynamic characteristics or sediment load can be of advantage because the EOF will give them the associate importance.

To sum up the proposed methodology and some recommendations in its use, here a list of the necessary steps:

In this study, the sinuosity of the river was considered as a variant to the Kinoshita formulation, and the results obtained by using four pairs of sine and cosine plus an independent term were more than acceptable. Although the proposed formulation does not include cinematic components for the meander behavior; it is a good approximation to the evolution of sinuosity in the river. Parameters regarding sediment transport have not been included either, but the method allows the incorporation of variables related to this phenomenon. A principal component analysis determines the position and importance of these variables.

The prioritization of the variables obtained through EOF showed clearly the existence of three groups of parameters, which altogether explain the behavior of the meandering of Cahuacan River. The first group is formed by the morphologic characteristics of the river. The second group corresponds to the hydrologic features of the basin, and the third one to the morphologic and geometric characteristics of the river.

The use of the maximum and minimum conditions that the meandering in the river has reached historically, gives us an adequate panorama to start with the estimation of the potentially damaged zone. In this way, the computation of the confidence limits, although in a statistical manner, constitutes a good tool to consolidate the arguments that define the zones at potential risk because of Cahuacan River meandering. Likewise, the stochastic simulation of the future conditions of the river allows the precise definition of the zones directly in the field. Finally, the proposed formulation is simple and can be implemented easily with a computer program. Besides, it allows the inclusions of other variables or a combination of them, that can describe for example, the effects that have not been possible to model with the traditional empirical formulations, like sediment distribution along the cross section of a river.