2.1 Flow field and sediment transport

Figure 2.

Sketch of a meandering reach with width variations and notation used in the mathematical model. In the cross-section the super-elevation of the free surface is represented (blue continuous line) with respect to the reference flow (blue dashed line). Notation as in Luchi et al. (Reference Luchi, Zolezzi and Tubino2011).

Table 1.

Perturbation solutions at different orders of approximation used in the present work. The perturbation parameters $(\unicode[STIX]{x1D708},\unicode[STIX]{x1D6FF})$ represent the dimensionless amplitudes of curvature and width oscillations, as defined in § 3.

We consider a periodic, sine-generated meandering river with variable width, a constant slope and uniform grain size  $d_{s}^{\ast }$ . Figure 2 depicts the situation. Hydrodynamics and sediment transport are described by the two-dimensional depth-averaged model derived from the quasi-steady three-dimensional Reynolds-averaged Navier–Stokes equations coupled with the Exner equation. We adapt the above model to an intrinsic curvilinear coordinate reference system $(s^{\ast },n^{\ast })$ , where $s^{\ast }$ represents the longitudinal coordinate measured along the channel axis and $n^{\ast }$ is the transverse coordinate orthogonal to $s^{\ast }$ (hereafter an asterisk denotes dimensional quantities).

Let us define

the intrinsic curvature of the river reach measured along the channel axis, where $\unicode[STIX]{x1D717}$ is the inflection angle of channel centreline with respect to a Cartesian axis $x^{\ast }$ ; $B^{\ast }(s^{\ast })$ the half-width of the channel, measured as orthogonal to the channel axis, and $\unicode[STIX]{x1D711}$ the local angle of the bank profile with respect to the channel axis, such that

We also conveniently define the following dimensionless quantities

where $(U^{\ast },V^{\ast })$ are the components of the depth-averaged velocity vector, $H^{\ast }$ denotes the water surface elevation with respect to a reference horizontal plane, $D^{\ast }$ is the water depth, $\boldsymbol{q}^{\ast }=(q_{s}^{\ast },q_{n}^{\ast })$ is the bedload discharge vector per unit width. All these quantities are spatially variable, i.e. they depend on the spatial coordinates $(s,n)$ . Furthermore, $\unicode[STIX]{x1D6E5}$ is the relative submerged density, and $g^{\ast }$ is the gravity. The transverse coordinate $n$ has been stretched to match the local channel width, therefore lateral boundary conditions are to be imposed at $n=\pm 1$ . Furthermore, $U_{0}^{\ast },D_{0}^{\ast }$ and $F_{0}$ denote the mean velocity, depth and Froude number of the uniform flow that would occur, for the given discharge, grain size and slope, in an infinitely long rectangular channel of constant width $2B_{0}^{\ast }$ .

Such reference conditions can be unambiguously described by fixing the values of the following dimensionless parameters:

namely the aspect ratio, the Shields parameter and the relative roughness, where $C_{0}$ is the friction coefficient of the reference flow.

We then split the half-channel width in the form

and conveniently define the longitudinal distributions of centreline curvature and half-width in the following dimensionless form:

where $R_{0}^{\ast }$ is a representative value of the channel radius of curvature (in the following twice the radius of curvature at the bend apex of a sine-generated channel centreline) and $B_{max}^{\ast }$ is the maximum value of the half-width within a meander wavelength. Substituting (2.5) and (2.6) into (2.1) and (2.2) we obtain

where two dimensionless parameters arise, namely the curvature ratio $\unicode[STIX]{x1D708}$ and the amplitude of the channel width oscillations $\unicode[STIX]{x1D6FF}$ ,

Hereafter, the formulation for the flow field and sediment transport follows closely the one reported in Zolezzi et al. (Reference Zolezzi, Luchi and Tubino2012b ), which the interested reader is referred to. Here we briefly outline the main assumptions and closures adopted.

1. (i) Secondary flows are accounted for through a procedure originally introduced by Zolezzi & Seminara (Reference Zolezzi and Seminara2001), in which both dispersive terms and bed shear stresses associated with spiral flows due to streamline curvature are parameterized by coefficients $k_{i},i=0,\ldots ,5$ (see Zolezzi & Seminara Reference Zolezzi and Seminara2001, pp. 189–191). The total streamline curvature ${\mathcal{C}}_{s}(s,n)$ is computed as the sum of the curvature of the channel centreline and the curvature associated with local streamline deviation from channel centreline in the form

   (2.9) $$\begin{eqnarray}{\mathcal{C}}_{s}=\unicode[STIX]{x1D708}{\mathcal{C}}-\frac{UV,_{s}-VU,_{s}}{U^{2}+V^{2}}.\end{eqnarray}$$
   By employing streamline curvature we account for centrifugal effects (Luchi et al. Reference Luchi, Zolezzi and Tubino2011), while normal and tangential Reynolds stresses are neglected under the assumption of slowly varying flow field (Zolezzi & Seminara Reference Zolezzi and Seminara2001).

2. (ii) The local bed shear stress vector $\unicode[STIX]{x1D749}(s,n)=(\unicode[STIX]{x1D70F}_{s}(s,n),\unicode[STIX]{x1D70F}_{n}(s,n))$ is defined as

   (2.10) $$\begin{eqnarray}(\unicode[STIX]{x1D70F}_{s},\unicode[STIX]{x1D70F}_{n})=C\sqrt{U^{2}+V^{2}}(U,V+\unicode[STIX]{x1D70F}_{k}),\end{eqnarray}$$
   where
   (2.11) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{k}=\left[\frac{U{\mathcal{C}}_{s}D}{\unicode[STIX]{x1D6FD}\sqrt{C_{0}}}\left(k_{3}+\frac{k_{5}D,_{s}}{\unicode[STIX]{x1D6FD}\sqrt{C_{0}}}\right)+\frac{k_{4}D^{2}(U{\mathcal{C}}_{s}),_{s}}{\unicode[STIX]{x1D6FD}^{2}C_{0}}\right],\end{eqnarray}$$
   represents the secondary flow contribution and $C$ is the local value of the friction coefficient. A commonly used logarithmic formula for rough open channel flow is employed to relate the spatially variable friction coefficient $C(s,n)$ to the local water depth values
   (2.12) $$\begin{eqnarray}C(s,n)=\left[6+2.5\ln \left(\frac{D(s,n)}{2.5d_{s}}\right)\right]^{-2}.\end{eqnarray}$$

3. (iii) The components of the local sediment discharge $[q_{s}(s,n),q_{n}(s,n)]$ , are expressed by means of the following relations which account for both shear stress and the downslope pull of gravity:

   (2.13a ) $$\begin{eqnarray}\displaystyle & \displaystyle q_{s}=\unicode[STIX]{x1D6F7}\left\{\cos \unicode[STIX]{x1D6FF}_{q}\left[1-\frac{1}{2}\left(\frac{\unicode[STIX]{x1D702},_{s}}{\unicode[STIX]{x1D6FD}(1+\unicode[STIX]{x1D708}\hat{n}{\mathcal{C}})}\right)^{2}\right]-\sin \unicode[STIX]{x1D6FF}_{q}\left(\frac{\unicode[STIX]{x1D702},_{s}\unicode[STIX]{x1D702},_{\hat{n}}}{\unicode[STIX]{x1D6FD}^{2}(1+\unicode[STIX]{x1D708}\hat{n}{\mathcal{C}})}\right)\right\},\qquad & \displaystyle\end{eqnarray}$$

   (2.13b ) $$\begin{eqnarray}\displaystyle & \displaystyle q_{n}=\unicode[STIX]{x1D6F7}\sin \unicode[STIX]{x1D6FF}_{q}\left[1-\frac{1}{2}\left(\frac{\unicode[STIX]{x1D702},_{n}}{\unicode[STIX]{x1D6FD}}\right)^{2}\right], & \displaystyle\end{eqnarray}$$
   where
   (2.14) $$\begin{eqnarray}\sin \unicode[STIX]{x1D6FF}_{q}(s,n)=\left(\frac{\unicode[STIX]{x1D70F}_{n}}{|\unicode[STIX]{x1D749}|}-\frac{r}{\unicode[STIX]{x1D6FD}\sqrt{\unicode[STIX]{x1D703}}}\unicode[STIX]{x1D702},_{\hat{n}}\right),\end{eqnarray}$$
   $\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D703})$ is the bed load function, $\unicode[STIX]{x1D702}=F_{0}^{2}H-D$ is the bed elevation, $\hat{n}=nB$ , $r$ is an empirical coefficient approximately equal to 0.5 (Talmon, Struiksma & Van Mierlo Reference Talmon, Struiksma and Van Mierlo1995), $\unicode[STIX]{x1D703}(s,n)$ is the local value of the Shields parameter and $\unicode[STIX]{x1D6FF}_{q}(s,n)$ denotes the local (spatially variable) angle formed by the local sediment transport rate vector with the streamwise coordinate $s$ .

   The formula of Parker (Reference Parker1990) is employed to compute sediment transport, as representative of a class of transport formulae in which the exponent of the local Shields parameter in the dimensionless bedload function tends to the classical $3/2$ value when $\unicode[STIX]{x1D703}$ is far enough from the critical conditions for sediment mobility, the same exponent increasing for $\unicode[STIX]{x1D703}$ approaching the critical value. Changing the bedload formula within the same class (e.g. Meyer-Peter & Müller Reference Meyer-Peter and Müller1948, Bagnold Reference Bagnold1966) has been shown not to qualitatively affect the outcomes of meander morphodynamic models (Zolezzi, Luchi & Tubino Reference Zolezzi, Luchi and Tubino2009). A local application of Parker’s formula, i.e. based on the local value of the Shields parameter $\unicode[STIX]{x1D703}(s,n)$ has been chosen because it yields analogous results to the other ones and it is more suitable for an analytical implementation, thanks to the formal absence of a critical threshold.

The resulting two-dimensional model for the flow field and sediment transport reads

where

arise as a consequence of the metric coefficients of the adopted curvilinear system and of the stretching of the transverse variable needed to accommodate to width variations. In the system (2.15), a comma has been used to denote the partial derivative of an unknown with respect to an independent variable $(s,n)$ . System (2.15) has already been employed in several studies (see e.g. Zolezzi et al. Reference Zolezzi, Luchi and Tubino2012b ) as the mathematical basis for deriving a spatially variable analytical solution for meandering rivers with variable channel width and curvature, and it contains all relevant physical effects relevant to such morphodynamic problem, which integrate inertial, redistribution and topographic effects due to the spatial variations of the forcing terms ( ${\mathcal{C}}(s),b(s)$ ) to the classical mechanisms controlling flow and bed topography in a constant curvature bend (transverse slope, centrifugal acceleration, superelevation of the free surface).

System (2.15) is associated with four boundary conditions in the lateral direction $n$ , imposing lateral impermeability to both flow and sediment discharge,

where $\boldsymbol{n}_{b}$ denotes the unit vector locally orthogonal to the banklines. In the longitudinal direction we impose periodic boundary conditions.

2.2 Planform evolution

To model the planform development of meandering rivers with spatially and temporally variable width, we integrate the planform evolution approach of Seminara et al. (Reference Seminara, Zolezzi, Tubino and Zardi2001) through an additional equation describing the temporal evolution of the spatial width oscillations.

The analysis is based on three key simplifying hypotheses. First, we assume the time scale of planform development to be sufficiently long to allow for the instantaneous adaptation of the flow field and bed topography to planform changes (curvature, width and length). This implies that as $\unicode[STIX]{x1D708}$ and $\unicode[STIX]{x1D6FF}$ coevolve in time, the flow and bed topography can be obtained from the steady solution of (2.15) and therefore they change in time only parametrically. Furthermore, on such a long time scale the bank migration process can be viewed as a continuous, rather than intermittent, process. Last, we focus our attention on width changes due to perturbations of flow field and bed topography with respect to the reference state. Therefore we assume that such a state, which changes slowly in time, does not produce any variation of the average channel width $2B_{0}^{\ast }$ .

Figure 3.

(a) Sketch of the model for planform development of channel axis and banks. (b) The model for bank development projected onto the curvilinear reference frame (the versor $\hat{\boldsymbol{n}}_{b}$ is locally orthogonal to the channel bank).

We refer to the geometrical construction reported in figure 3 to develop an evolutionary equation for the bank profile. Let us consider a local rate of migration of the channel bank $\unicode[STIX]{x1D701}_{b}$ , orthogonal to the channel axis, scaled with the reference velocity $U_{0}^{\ast }$ . The migration rate $\unicode[STIX]{x1D701}_{b}$ can be decomposed into a transversally antisymmetric contribution $\unicode[STIX]{x1D701}$ and a transversally symmetric contribution $\unicode[STIX]{x1D709}$ , such that

where $n=\pm 1$ denote the outer and inner banks, respectively. Therefore, $\unicode[STIX]{x1D701}$ represents the lateral migration rate of the channel axis, while $\unicode[STIX]{x1D709}$ represents the local migration rate of the bank with respect to the channel axis. The additional equation for $\unicode[STIX]{x1D709}$ is obtained in the following form (see figure 3 b):

where $t$ is the dimensionless time, scaled with the flow time scale $B_{0}^{\ast }/U_{0}^{\ast }$ .

Hence, the planform evolution model for a meandering river with freely variable width is described by the following system of equations:

Equation (2.20a ), which was first derived by Seminara et al. (Reference Seminara, Zolezzi, Tubino and Zardi2001), implies a cubic geometric nonlinearity in the evolution of the channel centreline. The integro-differential term in (2.20a ) accounts for the history of the deformation process of the channel centreline which affects the evolution of the channel axis. The cosine term in (2.20b ) is due to the different elongation of the bank with respect to the centreline and to the adopted definition of the bank migration rate $\unicode[STIX]{x1D709}$ as orthogonal to the channel axis. With $\unicode[STIX]{x1D701}$ the odd part of $\unicode[STIX]{x1D701}_{b}$ , we expect $\unicode[STIX]{x1D717}$ to be an odd function. This fact is enforced by the structure of (2.20a ), implying the cubic nonlinearity. On the other hand, we expect $\unicode[STIX]{x1D711}$ to be an even function, with $\unicode[STIX]{x1D709}$ being the even part of $\unicode[STIX]{x1D701}_{b}$ . The latter is enforced by the presence of the cosine term in the bank inflection equation (2.20b ).

When periodic, sine-generated meanders are of interest, the structure of equations (2.20) suggests a Fourier expansion of the variables $\unicode[STIX]{x1D717}$ and $\unicode[STIX]{x1D711}$ of the form

where the complex notation is conveniently employed, $\text{c.c.}$ represents the conjugate of a complex number and $\unicode[STIX]{x1D706}$ is the dimensionless meander wavenumber, scaled with the reference channel half-width $B_{0}^{\ast }$ . We note that $\unicode[STIX]{x1D706}$ is time dependent because meander sinuosity changes as curvature develops in time. We also note that only odd harmonics $(m=1,3,\ldots )$ need to be considered in (2.21a ) as a consequence of the cubic nonlinearity of (2.20a ).

Recalling (2.7b ), one can readily see that expansion (2.21b ) implies that the fundamental wavenumber of width oscillations is chosen as twice the wavenumber of curvature variations. This hypothesis relies on the assumption that each meander bend corresponds to an oscillation of channel width, which therefore displays two maxima over one meander wavelength. This is consistent with the few existing observations on real meanders (e.g. Seminara Reference Seminara2006, Luchi et al. Reference Luchi, Zolezzi and Tubino2011, Eke et al. Reference Eke, Czapiga, Viparelli, Shimizu, Imran, Sun and Parker2014a ) and is also implicit in Brice’s (Reference Brice1975) classification scheme for meandering rivers. Furthermore, the spatial structure set by (2.21b ) arises as a natural consequence of the autogenic nonlinear mechanism that produces width oscillations, as discussed in the next § 3.1.

If we keep only the first term of the expansion (2.21a ) we recover the classical sine-generated distribution. Using (2.7a ) we then obtain

Including also the second term of (2.21a ) we recover the Kinoshita (Reference Kinoshita1961) distribution of curvature, in dimensionless form

where

are the fattening and skewing coefficient, respectively.

Finally, a closure relation for the bank migration rate $\unicode[STIX]{x1D701}_{b}$ is required in order to solve (2.20). During the past years a number of models able to reproduce bank erosion and accretion mechanics leading to variations in channel width were developed. Mechanistic approaches for bank erosion processes were proposed by Duan & Julien (Reference Duan and Julien2005) and Chen & Duan (Reference Chen and Duan2006), who accounted for the sediment transport close to the banklines, and by Parker et al. (Reference Parker, Shimizu, Wilkerson, Eke, Abad, Lauer, Paola, Dietrich and Voller2011) and Motta, Abad & Garcia (Reference Motta, Abad and Garcia2012), who provided a physical description of the mechanics of bank processes such as the failure of slump blocks. Simplified models for bank erosion were proposed by Mosselman (Reference Mosselman1992), Darby et al. (Reference Darby, Alabyan and Van de Wiel2002), Jang & Shimizu (Reference Jang and Shimizu2005), Rüther & Olsen (Reference Rüther and Olsen2007), Darby, Rinaldi & Dapporto (Reference Darby, Rinaldi and Dapporto2007), while a rational, process-based description of bank accretion processes was provided by Zen et al. (Reference Zen, Zolezzi, Toffolon and Gurnell2016). However, including complex bank submodels into a perturbation scheme such as the one adopted in this paper is not straightforward. Simplified models able to approximately describe the key processes on longer time scales were widely used in the literature (Partheniades & Paaswell Reference Partheniades and Paaswell1970; Ikeda et al. Reference Ikeda, Parker and Sawai1981; Seminara et al. Reference Seminara, Zolezzi, Tubino and Zardi2001; Zolezzi et al. Reference Zolezzi, Luchi and Tubino2012b ). They generally neglect bank boundary layers and the corresponding distribution of stresses at the bank and simply relate the lateral bank shift to the excess (defect) of flow velocity or stress at the bank toe. The latter assumption can be reasonable to estimate bank retreat even at the time scale of single flood events driving morphological change, while it strongly oversimplifies the complex depositional dynamics of bank advance that involves vegetation encroachment. However, besides the complexity of including a bank erosion/accretion model with a high level of detail into an analytical framework, we consider that on the long time scale of planform evolution of meandering rivers the differences in the processes that locally rule the advance and retreat of channel banks tend to be mitigated.

In Seminara et al. (Reference Seminara, Zolezzi, Tubino and Zardi2001) the centreline’s migration rate is assumed to be proportional to the difference between the velocity at the outer bank and that at the inner bank. A novel element in our approach, which is needed to allow width to locally evolve, is the possibility for each bank to migrate independently. The concept underpinning such novel formulation is graphically illustrated in figure 4. Namely, the migration rate of each bank is set as proportional to the difference between the velocity at the bank and a threshold value, here assumed to coincide with the reference flow velocity $U_{0}^{\ast }$ . In dimensionless form, such a bank migration model reads

where $E$ is a lateral migration coefficient, which typically takes values in the range $10^{-7}{-}10^{-8}$ (Camporeale et al. Reference Camporeale, Perona, Porporato and Ridolfi2005), and the factor 2 has been added for consistency with the traditional method. It must be noted that the time scale of simulated meander planform development in a real case study crucially depends on the choice of the lateral migration coefficient $E$ , which has been shown to actually represent a calibration parameter (Crosato Reference Crosato2009) for model testing against multi-temporal observations on channel centreline evolution.

We also note that the ratio between the planform time scale $T_{E}\sim B_{0}^{\ast }/EU_{0}^{\ast }$ and the bed time scale $T_{B}\sim B_{0}^{\ast }D_{0}^{\ast }/\sqrt{g^{\ast }\unicode[STIX]{x0394}{D_{s}^{\ast }}^{3}}\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D703}_{0})$ takes the order of magnitude of $10^{1}{-}10^{3}$ in natural meandering rivers, which justifies the assumption of instantaneous adaptation of flow and bed topography to planform changes.

Figure 4.

Sketch of decoupled bank migration model and notation.

Equations (2.25) are then combined is such a way that half their difference and their average give rise to the centreline migration rate $\unicode[STIX]{x1D701}$ and to the local width change rate $\unicode[STIX]{x1D709}$ , respectively,

Equation (2.26a ) recovers the classical lateral migration law for the channel centreline originally proposed by Ikeda et al. (Reference Ikeda, Parker and Sawai1981), while (2.26b ) describes the temporal evolution of the spatial channel width oscillations. We note that setting the threshold velocity at the fixed value $U_{0}^{\ast }$ implies that we neglect the adaptation of the average channel width to the slow time change of the reference flow due to meander elongation. Although the latter process could be easily incorporated into the model (see for instance Eke et al. (Reference Eke, Czapiga, Viparelli, Shimizu, Imran, Sun and Parker2014a )), we adopt here such approximation in order to single out the role of flow and bed perturbations.

3.1 Perturbation approach

The two-dimensional hydrodynamics and sediment transport problem (2.15) is solved analytically by means of a second-order two parameter ( $\unicode[STIX]{x1D708}$ , $\unicode[STIX]{x1D6FF}$ ) expansion. The analysis is based on the recognition that the values of $\unicode[STIX]{x1D708}$ and $\unicode[STIX]{x1D6FF}$ are typically much smaller than unity in real meandering rivers. For instance, Zolezzi et al. (Reference Zolezzi, Luchi and Tubino2012b ) examined the spatial series of curvature and width for the class-C bends reported in the large dataset of Lagasse et al. (Reference Lagasse, Zevenbergen, Spitz, Thorne, Associates and Collins2004) showing that average $\unicode[STIX]{x1D6FF}$ values are typically slightly less than $0.1$ , with values up to $\unicode[STIX]{x1D6FF}\simeq 0.2$ for very few bends, while $\unicode[STIX]{x1D708}$ falls in the range 0.05–0.2 and rarely exceeds $0.25$ .

Hence, we employ the following expansion of the unknown variables:

where $\boldsymbol{V}_{ij}^{(k)}=(U_{ij}^{(k)},\boldsymbol{V}_{ij}^{(k)},H_{ij}^{(k)},D_{ij}^{(k)})^{\text{T}}$ , $\boldsymbol{V}_{0}=(1,0,H_{0}(s),1)^{\text{T}}$ and the overbar ( $\,\bar{\;}\,$ ) denotes the conjugate of a complex number. Note that analogous perturbation expansions are performed for the Shields parameter, friction coefficient, bed shear stresses, sediment fluxes, which are therefore variable along the meandering channel as they are functions of ( $s,n$ ) through their nonlinear dependences on the local values of the four fundamental unknowns ( $U$ , $V$ , $H$ , $D$ ), as expressed by the closure relations (2.9)–(2.13) and on the local curvature values, as expressed by (3.1). In order to keep the analysis as simple as possible, yet retaining the fundamental interactions between curvature and width variations, we have included the first two terms of the curvature distribution arising from (2.21a ) and the fundamental component of width oscillations

Furthermore, in (3.1) we have only considered the nonlinear terms ( $\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FF}$ , $\unicode[STIX]{x1D708}^{2}$ ) that arise from the interaction between the fundamental components ( $e_{1}$ , $e_{2}$ ) of curvature and width distributions and we have neglected the $O(\unicode[STIX]{x1D708}^{2},\unicode[STIX]{x1D6FF}^{2})$ nonlinear distortion of the basic reference flow.

Finally, let us briefly explain why the second-order contribution $O(\unicode[STIX]{x1D6FF}^{2})$ has not been incorporated. In a dominant-mode periodic analysis, we only consider the first harmonic of each longitudinal signal. The $O(\unicode[STIX]{x1D708})$ contribution has a $\unicode[STIX]{x1D706}$ periodicity, while the $O(\unicode[STIX]{x1D6FF})$ contribution has a $2\unicode[STIX]{x1D706}$ periodicity; $N$ th-order terms have $N\unicode[STIX]{x1D706}$ periodicity, i.e. $O(\unicode[STIX]{x1D708}^{2})$ has $2\unicode[STIX]{x1D706}$ periodicity (hence it interacts with the width oscillation), while $O(\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FF})$ has both a $\unicode[STIX]{x1D706}$ and a $3\unicode[STIX]{x1D706}$ periodicity; $O(\unicode[STIX]{x1D6FF}^{2})$ has a $4\unicode[STIX]{x1D706}$ periodicity, hence it does not affect the fundamental and can be left out from the analysis.

The solution (3.1) is used to feed the planform evolution model (2.20) through the bank and centreline migration relations (2.26).

The $O(\unicode[STIX]{x1D708})$ linear solution represents the classical flow field solution forced by the presence of linearly forced bars in a regular sequence of sine-generated meander bends; its analytical solution has been employed for a long time in the analysis of meander stability (Blondeaux & Seminara Reference Blondeaux and Seminara1985) and evolution (Seminara et al. Reference Seminara, Zolezzi, Tubino and Zardi2001). Its transversal structure is antisymmetric and oscillates longitudinally with the same frequency as the curvature.

The $O(\unicode[STIX]{x1D6FF})$ linear solution was first derived by Repetto et al. (Reference Repetto, Tubino and Paola2002) and represents the flow-bed topography response to a periodically varying width along a straight reach. Its transverse structure is symmetrical, giving rise to a spatial sequence of mid-channel bars.

The $O(\unicode[STIX]{x1D708}^{2})$ nonlinear self-interaction has a transversal symmetric structure and it may be the cause of the presence of mid-channel bars in meandering rivers, as shown by Luchi et al. (Reference Luchi, Hooke, Zolezzi and Bertoldi2010a ). Being transversally symmetrical, it is the leading curvature perturbation that can interact with width variation in the processes of local widening and narrowing because of the same transverse structure. Moreover, its wavenumber intrinsically emerges as twice that of the curvature, consistently with the assumption made on channel width in (2.21b ).

The $O(\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FF})$ nonlinear mixed interaction is instead antisymmetrical and embodies the effect of width variations on bend dynamics (Luchi et al. Reference Luchi, Zolezzi and Tubino2011), because its spatial structure exactly corresponds to that of the $O(\unicode[STIX]{x1D708})$ first-order solution. Expansion (3.1) suggests that width oscillations also affect the development of the third longitudinal harmonic $e_{3}$ of curvature distribution, which is responsible for fattening and skewing of meander shape, through the term $\boldsymbol{V}_{11}^{(3)}$ . The role of such contribution, which was not previously investigated, may turn out to be relevant in the final stage of meander evolution as the role of $e_{3}$ component becomes significant at relatively high sinuosity.

Table 2.

Symmetry properties and wavenumbers of the perturbation solutions for the longitudinal velocity $U_{ij}$ .

Symmetry properties of the perturbation solution at different orders are summarized in table 2. We note that, from the mathematical point of view, the channel width of meandering rivers is continuously forced towards spatial oscillations ( $O(\unicode[STIX]{x1D708}^{2})$ effect). In summary, mutual feedbacks between width and curvature are primarily accounted for in the $O(\unicode[STIX]{x1D708}^{2})$ solution and in $O(\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FF})$ solution (Luchi et al. Reference Luchi, Zolezzi and Tubino2011). More details on the relation between the solution at different perturbation orders can be found in Zolezzi et al. (Reference Zolezzi, Luchi and Tubino2012b ).

3.2 Perturbation solution

Substituting the expansion (3.1) into the system (2.15) we derive an ordinary differential system in the transverse direction for each unknown at every order of the perturbation expansion, which can be conveniently written in matrix form for any system of order $O(\unicode[STIX]{x1D708}^{i}\unicode[STIX]{x1D6FF}^{j})$ ,

where the coefficients $a_{lk}^{ij}$ are analogous to those computed by Luchi et al. (Reference Luchi, Zolezzi and Tubino2011), while the coefficients $b_{l}^{ij}$ are the forcing terms arising from nonlinear interaction between lower-order terms. Each system of equations can be easily solved analytically once lateral boundary conditions are imposed (2.17); they take the form

The solution of (3.3) yields the transverse structure of the perturbations of flow field and bed topography at the different orders of approximation. Here we focus on the longitudinal velocity perturbations, which in the present approach (2.25) are responsible for bank migration. In particular, solutions $U_{10}$ and $U_{11}$ are antisymmetric and therefore involved in bend amplification, while $U_{20}$ and $U_{01}$ , both being symmetric, regulate autogenic width adjustments along meander bends.

Substituting the perturbation solutions into (2.26) we obtain the components of the lateral migration rate $\unicode[STIX]{x1D701}$ of the channel centreline and of the local rate $\unicode[STIX]{x1D709}$ of width change in the form,

where, according to (3.2), we have used the following decompositions

It is important to note that the antisymmetric perturbations of the flow field act on the centreline migration rate $\unicode[STIX]{x1D701}$ , while symmetrical perturbations either force or damp the local widening rate $\unicode[STIX]{x1D709}$ . Therefore, from a linear point of view, the curvature forced $\unicode[STIX]{x1D708}$ -perturbation produces effects on meander amplification, while the $\unicode[STIX]{x1D6FF}$ -perturbation, forced by width variation, affects the width changes. From the two-dimensional analysis of Repetto et al. (Reference Repetto, Tubino and Paola2002), the $O(\unicode[STIX]{x1D6FF})$ linear solution for the longitudinal velocity at the bank (the first term on the right-hand side of (3.5c )) turns out to invariably produce a self-suppressing effect, at least for the values of the wavenumber relevant for the present analysis. Therefore, the only mechanism able to produce width variation is the autogenic effect related to the second-order solution for curvature perturbation, $O(\unicode[STIX]{x1D708}^{2})$ .

Expansions (2.21) and (3.6) are finally incorporated into the planimetric evolution model (2.20). The resulting system of ordinary differential equations eventually reads

where the erosion coefficient $E$ has been absorbed into the dimensionless time, $T=2Et$ , consistently with (2.25), and

Equations (3.7a,b ) describe the temporal evolution of the centreline curvature, while (3.7c ) quantifies meander elongation as curvature and width coevolve. These equations are formally identical to those derived by Seminara et al. (Reference Seminara, Zolezzi, Tubino and Zardi2001) for constant width meanders. However, the terms on the right-hand side include the effect of width oscillations which is felt through the dependence on $\unicode[STIX]{x1D6FF}$ of the components of the migration rate $\unicode[STIX]{x1D701}_{1},\unicode[STIX]{x1D701}_{3}$ (3.5a,b ). Therefore, their solution is coupled with that of (3.7d ) describing the temporal evolution of width oscillations, which in turn is affected by curvature through the dependence of $\unicode[STIX]{x1D709}_{2}$ on the autogenic $O(\unicode[STIX]{x1D708}^{2})$ term appearing in (3.5c ). As a limiting case, system (3.7) can be reduced to the classical constant width case by setting $\unicode[STIX]{x1D6FF}=0$ . On the other hand, the evolutionary solution for a straight channel with periodic width oscillations can be obtained by setting $\unicode[STIX]{x1D708}=0$ .

The solution of (3.7) has been computed by means of a variable-coefficient ordinary differential equation system solver with fixed-leading-coefficient implementation, backward differentiation method and adaptive numerical time step (Brown, Byrne & Hindmarsh Reference Brown, Byrne and Hindmarsh1989) in order to achieve both stability and a good level of accuracy (Jones et al. Reference Jones, Oliphant and Peterson2001). The initial condition has been set in terms of a small amplitude sine-generated meander, whose initial wavenumber is $\unicode[STIX]{x1D706}_{0}$ , with given initial values of the reference flow parameters ( $\unicode[STIX]{x1D6FD}$ , $\unicode[STIX]{x1D703}_{0}$ and $d_{s}$ ), while $\unicode[STIX]{x1D711}$ has been set equal to zero. At each time step the complex coefficients that describe the distribution of curvature and width oscillations in (3.2) are computed through (2.7) in the following form

Recalling the definitions (2.8) and (2.22) of the parameters $\unicode[STIX]{x1D708}$ and $\unicode[STIX]{x1D6FF}$ , one can readily prove that $|C_{1}|=|b_{2}|=1$ .