1. Introduction
River meanders are amongst the most common terrestrial landforms (Leopold et al. Reference Leopold, Wolman, Wolman and Wolman1957). Research on meandering has been motivated by scientific curiosity (Einstein Reference Einstein1926); by civil engineering considerations such as flood control, the preservation of land and infrastructure, artificial channel construction, bank stabilization and sediment dredging (Bridge Reference Bridge1993; Güneralp et al. Reference Güneralp, Abad, Zolezzi and Hooke2012); and by economic interests, especially of the petroleum industry in ancient meandering rivers whose deposits have since formed oil reservoirs (Sun et al. Reference Sun, Meakin, Jøssang and Schwarz1996).
Meandering rivers are dynamic landforms that migrate across floodplains. Their morphodynamic evolution is governed by complex interactions among sediment, flow, bed topography and channel curvature. Meandering does not result from obstacle avoidance: given a flat, slightly tilted plain of appropriately cohesive sediment, meanders will develop absent natural obstacles. Under the right conditions, the centreline of the erodible channel migrates laterally, and near-bank hydrodynamics lead to higher streamwise velocity near the outer bank of a curve and preferential erosion that results in meandering. Periodic oxbow cutoff events occur when meanders become sufficiently large that the channel intersects with itself, causing sudden reductions of length and sinuosity. When a cutoff occurs, the abandoned loop may be filled through the deposition of sediment carried by overbank flows associated with floods, though not all oxbows are filled by that process. A second important cutoff type, a chute cutoff, occurs when a river cuts a new channel that bypasses its own existing meander loop. Multiple chute cutoffs may initiate braiding. While this last phenomenon is not the subject of the current paper, it is widespread. See Kleinhans (Reference Kleinhans2010) for a discussion of braided rivers and chute cutoffs. An excellent summary of the basic dynamics of slender meandering rivers is provided in the review article by Seminara (Reference Seminara2006).
1.1. Modelling of meandering rivers
Bend instability has long been identified as one of two principal mechanisms by which meanders develop; under this paradigm, an initially straight channel with small perturbations will develop into meanders (Seminara Reference Seminara2006). In 1981, a landmark publication by Ikeda, Parker & Sawai (IPS) proposed a set of governing equations that describe the evolution of a meandering channel centreline as a function of channel curvature. In a series of four papers in 1989, Johannesson and Parker revisited the modelling of secondary flow in curved rivers and amended the treatment of the secondary flow presented by IPS. Johannesson & Parker (Reference Johannesson and Parker1989a) studied the basic mechanisms for secondary flow, while Johannesson & Parker (Reference Johannesson and Parker1989b) corrected errors in the IPS formulation and included the correct form for the influence of secondary flow on momentum transport in the IPS model of meander dynamics. Surprisingly, this paper found that the original IPS partial differential equations for the river planform required only a modified definition for a single coefficient 
$A$ to implement the corrections needed in IPS. Thus the original linear stability theory (with the new definition for 
$A$) carries forward and dictates that relatively long-wavelength meanders should grow and migrate, and that a dominant wavelength should emerge amongst the meanders (Ikeda et al. Reference Ikeda, Parker and Sawai1981). Both of these predictions are consistent with field observations and simulation-based works. Next, Johannesson & Parker (Reference Johannesson and Parker1989c) generalized their treatment for more complex bottom topography than in the IPS theory (from previous restriction to lateral beds of fixed transverse slope at any river cross-section). This allowed for prediction of additional phenomena including alternating channel bars and overdeepening near river bends. In their fourth paper, Parker & Johannesson (Reference Parker and Johannesson1989) gave a comparison of several meandering models from five different research groups, citing similarities and differences with respect to treatment of the fluid velocity and bed topography.
Enhanced modelling of bed topography is important because an eroding river bed provides a second independent mechanism which may lead to meandering in rivers. Blondeaux & Seminara (Reference Blondeaux and Seminara1985) developed ‘bar’ theory, noting that a straight channel with a perturbed bed will develop bedforms, or bars, at the floor of the channel. Their theory dictates that a critical bar wavelength emerges, at which maximum growth occurs, and that those bars should migrate downstream. From linear theory for bars, a ‘neutral stability’ wavelength also emerges, which corresponds to a bar growth rate of zero. Reference Blondeaux and SeminaraBlondeaux & Seminara recognized that at some critical width-to-depth ratio the ‘neutral stability’ bar wavelength corresponds to the ‘maximum growth’ wavelength from bend theory such that a resonance phenomenon occurs. In a sinuous channel, the bars that emerge are quasisteady in the sense that they neither grow nor propagate, except by the slow bend-erosion process (Blondeaux & Seminara Reference Blondeaux and Seminara1985; Monegaglia et al. Reference Monegaglia, Zolezzi, Tubino and Henshaw2017).
In Reference Zolezzi and Seminara2001, Zolezzi & Seminara provided a notable advance in the modelling of the meandering of slender rivers as subject to the effects of both centreline curvature and changes in bottom topography. They considered full expansions for both the streamwise and transverse variation in velocity and bedform to provide detailed solutions for the linearized regime. Seminara et al. (Reference Seminara, Zolezzi, Tubino and Zardi2001) built on this extended model, presenting a procedure for numerical simulations of nonlinear equations and river planfom evolution. Seminara et al. (Reference Seminara, Zolezzi, Tubino and Zardi2001) provided theory describing meander behaviour under both subresonant and superresonant conditions (where the transition between the two depends on the width of the river as compared with the resonant width). They provide detailed discussion of the differences associated with upstream and downstream influence in the two regimes. Further discussion of these phenomena is in the review article cited above Seminara (Reference Seminara2006).
Camporeale et al. (Reference Camporeale, Perona, Porporato and Ridolfi2007) give an excellent overview of the hierarchy of the meandering dynamics models described above using a common notation to facilitate comparisons. They use the shorthand IPS, JP and ZS to describe the models of Ikeda et al. (as corrected by Johannesson & Parker (Reference Johannesson and Parker1989b)), of JP in Johannesson & Parker (Reference Johannesson and Parker1989c) and of ZS in Zolezzi & Seminara (Reference Zolezzi and Seminara2001), respectively. They further generated numerical results from each of the three models and compared the resulting planforms; additional comparisons were given in Camporeale et al. (Reference Camporeale, Perona, Porporato and Ridolfi2005).
The three linear-velocity models differ in how they compute the transverse velocity profile for the streamwise velocity in the channel. We note that all references to IPS in the current paper refer to the IPS equations with corrections introduced by Johannesson & Parker (Reference Johannesson and Parker1989b). As noted, the JP and ZS models include more general modelling of the bottom topography and include proper treatment of the full bend/bar mechanisms. These give better modelling in the subresonant regime and are the only appropriate models to use in studying rivers which may feature both subresonant and superresonant meandering conditions – as IPS is able to model only subresonant conditions. All three models have been used extensively by researchers, see e.g. Sun et al. (Reference Sun, Meakin, Jøssang and Schwarz1996) and numerous citations in Camporeale et al. (Reference Camporeale, Perona, Porporato and Ridolfi2007) with the widespread early use of the IPS model leading to its inclusion in numerous later articles for comparison purposes. Numerous efforts (Camporeale et al. Reference Camporeale, Perona, Porporato and Ridolfi2007) have demonstrated that the ZS formulation is capable of modelling upstream migrating/downstream facing meanders (i.e. superresonant conditions) as well as downstream migrating/upstream facing meanders (i.e. subresonant conditions). At this time, the JP and ZS models should be regarded as the standard models for capturing the physics of meandering dynamics for slender rivers. In the present paper, we focus our attention on the IPS model, despite its inability to capture the full physics in all regimes. Our motivation for this choice is that it facilitates a formal mathematical formulation which addresses issues of mathematical uniqueness and proper specification of boundary conditions in finite domains, and provides a proper platform for development of robust, stable, high-accuracy numerical simulations. These fundamental aspects of modelling for meandering river dynamics have not been fully addressed in previous work with any of the three slender river models – IPS, JP, ZS. These issues are more fully discussed in §§ 1.3 and 1.4 below. Given the validity of the IPS model in capturing the physics of subresonant meandering for rivers, our novel results for the subresonant dynamics in finite domains add to our understanding of meandering rivers in these conditions. Possible steps for the future extension of the current mathematical formulation to the JP and ZS models is discussed briefly in § 1.4 below.
Camporeale, Perucca & Ridolfi (Reference Camporeale, Perucca and Ridolfi2008) investigated the statistical differences in river planforms between IPS and ZS models in the subresonant regime when many neck cutoffs occur and found that cutoffs acted to reduce the statistical differences between IPS and ZS planforms. Frascati & Lanzoni (Reference Frascati and Lanzoni2009) discuss Camporeale's results and agreed that the dynamic cutoff mechanisms likely smooth out differences associated with IPS and ZS level approximations. They did, however, find subtle differences in statistical measures of length scales generated by the different approximations. In addition to their discussion of the three velocity models discussed above, Camporeale et al. (Reference Camporeale, Perona, Porporato and Ridolfi2007) comment on various limitations in all of the models including restriction to the shallow-water approximation with two-dimensional flow fields. They further discuss three-dimensional numerical efforts to model the hydrodynamics as well as analytical efforts for the same.
In nearly all of the classic meandering models, a highly simplified model is adopted to predict the dynamics of the migration of the river banks. By far the most common (also adopted by IPS, JP and ZS) is a simple erosion law wherein the normal velocity of the bank is proportional to the local depth-averaged velocity near the bank with an empirical erosion coefficient 
$E$. This erosion model first popularized by Ikeda et al. (Reference Ikeda, Parker and Sawai1981) was independently suggested by Hasegawa (Reference Hasegawa1977); in recognition, it has been referred to as the HIPS model in Parker et al. (Reference Parker, Shimizu, Wilkerson, Eke, Abad, Lauer, Paola, Dietrich and Voller2011). In a major effort to improve the modelling of the bank dynamics, a large group of experimentalists, numerical analysts and modellers (Parker, Shimizu, Wilkerson, Eke, Abad, Lauer, Paola, Dietrich and Voller) set out to develop a new framework to more realistically model the complex dynamics of the migration of the eroding bank and the depositing bank. See Parker et al. (Reference Parker, Shimizu, Wilkerson, Eke, Abad, Lauer, Paola, Dietrich and Voller2011). The new formulation of these authors was intended to replace the HIPS model of erosion with more realistic dynamics and to better evaluate the dynamics of changes in river width. Modern numerical simulations which build on the bank models of Parker et al. (Reference Parker, Shimizu, Wilkerson, Eke, Abad, Lauer, Paola, Dietrich and Voller2011) include Asahi et al. (Reference Asahi, Shimizu, Nelson and Parker2013) and Weisscher, Shimizu & Kleinhans (Reference Weisscher, Shimizu and Kleinhans2019) who demonstrated that sustained dynamic inflow perturbation at the upstream entry and bar-floodplain conversion are required for dynamic meandering in a domain of finite length.
1.2. Experimental studies
Fundamental gaps remain in our understanding of the mechanisms that drive meandering. Despite the proliferation of fieldwork-derived data sets and ever improving modelling and simulation efforts, experimental efforts to reproduce dynamic, self-maintaining meandering channels have been largely unsuccessful. Efforts to reproduce meandering rivers in a laboratory setting date back to 1945 (Reference FriedkinFriedkin). Whereas braided rivers are relatively easy to produce in the laboratory (Friedkin Reference Friedkin1945), dynamic, self-maintaining meandering rivers with cutoffs have proven notoriously challenging. For our interests, the most notable aspect of most experimental studies is that nearly all used an initially straight channel as the initial condition, along with a fixed upstream feed (Smith Reference Smith1998) that might be set at a static angle (Friedkin Reference Friedkin1945; Jin & Schumm Reference Jin and Schumm1987; Gran & Paola Reference Gran and Paola2001; Peakall, Ashworth & Best Reference Peakall, Ashworth and Best2007; Tal & Paola Reference Tal and Paola2007; Braudrick et al. Reference Braudrick, Dietrich, Leverich and Sklar2009; Kleinhans et al. Reference Kleinhans, van Dijk, van de Lageweg, Hoendervoogt, Markies and Schuurman2010; Tal & Paola Reference Tal and Paola2010).
In Reference van Dijk, Lageweg and Kleinhans2012, van Dijk et al. departed from the standard fixed injection point and used a dynamic injection point traversing the channel. In so doing, they achieved the first high-sinuosity migrating meanders with scroll bars and (infrequent) cutoffs (van Dijk et al. Reference van Dijk, Lageweg and Kleinhans2012; Kleinhans et al. Reference Kleinhans, Schuurman, van Dijk, van de Lageweg, Castelle and Ruessink2013). Following van Dijk et al. (Reference van Dijk, Lageweg and Kleinhans2012), simulation-based work by Schuurman et al. (Reference Schuurman, Shimizu, Iwasaki and Kleinhans2016) demonstrated that a necessary condition for sustained meandering is a perturbation at the upstream boundary. These results are significant, as they suggest that the fixed inlet of previous experimental efforts may be the most significant reason for the notorious difficulty of creating dynamic, self-sustaining meandering channels in flumes. As noted by Reference van Dijk, Lageweg and Kleinhansvan Dijk et al., their results are consistent with the framework of theory developed by Lanzoni & Seminara (Reference Lanzoni and Seminara2006). From the asymptotic models for slender rivers, Reference Lanzoni and SeminaraLanzoni & Seminara demonstrated that bend instability in many natural systems is a convective instability which will convect in only one direction. For ‘sub-resonant’ rivers, that direction is downstream. A flume with a fixed injection point upstream should therefore result in a straight channel propagating downstream (Schuurman et al. Reference Schuurman, Shimizu, Iwasaki and Kleinhans2016), as has been seen in many experimental flume experiments (Friedkin Reference Friedkin1945; Smith Reference Smith1998; Braudrick et al. Reference Braudrick, Dietrich, Leverich and Sklar2009). Following these efforts, Weisscher et al. (Reference Weisscher, Shimizu and Kleinhans2019) provided numerical results demonstrating that the period of the inflow perturbation affects how rapidly meanders develop from a straight channel.
The need for a dynamic upstream boundary for sustained meanders in finite domains has now been shown in both experimental and numerical work for rivers in subresonant conditions. It has been shown numerically that the perturbation frequency determines observed wavelengths near the inlet, while for nonlinear amplitude, wavelengths evolve to the periodic maximum-growth wavelength as the meanders migrate downstream. The potential for more elaborate meandering dynamics to manifest in experiments with dynamic forcing offers the possibility for valuable comparisons between experiments and modern numerical simulations for finite width rivers and more realistic bank dynamics (Asahi et al. Reference Asahi, Shimizu, Nelson and Parker2013; Weisscher et al. Reference Weisscher, Shimizu and Kleinhans2019). Parker et al. (Reference Parker, Shimizu, Wilkerson, Eke, Abad, Lauer, Paola, Dietrich and Voller2011) introduced a new system for modelling complex bank dynamics and supported their models with detailed observations from field studies. Weisscher et al. (Reference Weisscher, Shimizu and Kleinhans2019) conducted detailed simulations using more complex bank dynamics and investigated the interaction of hydrodynamic effects with geological and biological mechanisms for stabilization of emerging floodplains. These simulations were conducted to offer insight for further experiments in their laboratory. Given the critical role of sustained dynamic perturbations at the system boundaries, there is significant interest to investigate different modes of boundary perturbation to find the most versatile scheme amenable to experimental implementation.
1.3. Numerical methods
Turning to standard practices in numerical algorithms for meandering dynamics, multiple authors have noted that even small changes in boundary conditions, initial conditions or cutoff algorithms may lead to vastly different channel forms (Crosato Reference Crosato1990; Howard Reference Howard1992; Sun et al. Reference Sun, Meakin, Jøssang and Schwarz1996; Sun, Meakin & Jøssang Reference Sun, Meakin and Jøssang2001). Authors who documented efforts to implement the nonlinear equations have further cited issues with numerical instabilities and the need for averaging from point to point, smoothing, numerical filtering and fitting polynomials (Crosato Reference Crosato1990; Howard Reference Howard1992; Sun et al. Reference Sun, Meakin, Jøssang and Schwarz1996, Reference Sun, Meakin and Jøssang2001; Camporeale et al. Reference Camporeale, Perona, Porporato and Ridolfi2007; Frascati & Lanzoni Reference Frascati and Lanzoni2009; Shimizu et al. Reference Shimizu2020). Reference CrosatoCrosato is one of the few authors to address in detail the numerical artefacts these practices introduce, and to emphasize the need for testing numerical algorithms. See Crosato (Reference Crosato1990) (figure 26, p. 72) for a sample of three runs testing planform convergence vs time step, and (figure 27, p. 74) for a filtering/smoothing test. In further work, Crosato (Reference Crosato2007) presented convergence studies showing that, for three conceptually different meander models, the choice of numerical filter affects the size, form and migration rate of the meander channels. It follows that efforts to characterize planform statistical behaviour from numerics are likely complicated by unaccounted-for numeric artefacts. Crosato (Reference Crosato2007) presents additional tests of smoothing and filtering in meander studies and notes that ‘Meander migration models include an as yet poorly investigated source of numerical errors related to the computation of the channel curvature, which are amplified by the procedure of adding and deleting grid points as the river planform evolves’. In fact, throughout all of the numerical studies on meandering rivers with IPS, JP and ZS models, there is not a single case where convergence results are presented which demonstrated clear error estimates and that showed convergence with increasing spatial and temporal resolution. To be clear, individual pieces of an algorithm may be higher than first-order accurate and tested carefully (e.g. see time step selection, or high-order polynomial fits in Sun et al. Reference Sun, Meakin, Jøssang and Schwarz1996), but subsequent steps in the time step (averaging, filtering, smoothing, interpolating, quadratures) likely render the final results first order at best.
Together, these studies point to a need for investigation of rapid, efficient stable numerical methods to eliminate or reduce the numerical uncertainties in model predictions and to accurately represent the physics of the governing differential equations. Given the care and effort devoted to the derivation, modelling and evaluation of empirical approximations in the IPS, JP and ZS models, one would hope to obtain the cleanest vision of the physical impacts, unobscured by numerical artefacts or excess smoothing required for stability. The two more sophisticated models in particular rely on high-order derivatives for their improved modelling of the physics, and it is possible that some of the advantages of the enhanced models are lost to deficiencies in smoothing/filtering/averaging in the numerical implementations.
In addition to the numerical uncertainty introduced by the limitations in numerical algorithms described above, there is the issue that no author has provided a clear definitive statement of the boundary conditions required for a unique solution of the meandering dynamics equations. A note should be made of the valuable insight given by Reference Lanzoni and SeminaraLanzoni & Seminara on the dramatic impact of subresonant and superresonant river dynamics on the number of boundary conditions required at upstream and downstream boundaries.
Despite the valuable guidance on the number and placement of boundary conditions, Reference Lanzoni and SeminaraLanzoni & Seminara did not prescribe specific sets of boundary conditions which would guarantee unique solutions for initial-boundary-value problems. Nearly all authors place the boundaries of the systems far upstream and far downstream so that conditions at the boundaries will not be dominant in the middle of the channel where they focus their analysis. This approach is not applicable if one is studying the dynamics in a domain of finite length. We emphasize that there has been no study to date demonstrating converged solutions in a non-periodic finite domain. The specification of proper boundary conditions is a necessity for successful resolution of this issue. Reference should be made here to Crosato (Reference Crosato1990) whose convergence studies come the closest to this standard.
1.4. Focus of the present paper
The major focus of our paper is the dynamics of meandering rivers in finite domains in contrast to strictly periodic domains or domains of near infinite extent. Our goal is to set out a formal mathematical model with a rigorous statement of boundary and initial conditions such that a well-posed mathematical problem is presented for which uniqueness results and stability theory can be demonstrated in the case of the linear system. With a well-posed initial-boundary-value problem, the powerful tools for numerical solution of partial differential equations may be applied to develop robust, stable numerical algorithms for both linear and nonlinear cases with documented numerical convergence results and properties.
Models and simulations in finite domains are of interest for three primary reasons: (i) all laboratory results are inherently limited to flumes of significantly bounded length and subject to boundary conditions imposed by inlets and downstream exit conditions of the flumes; (ii) river engineering projects often involve managing short finite stretches of rivers and predicting the ensuing bank dynamics over periods of a few years to decades; (iii) finite domains with well-posed boundary and initial conditions provide an ideal environment for developing, testing and validating numerical algorithms for meandering dynamics. Such an environment could provide a platform for benchmarks against which researchers can test demonstrate the accuracy, stability and robustness of their numerical algorithms.
In the present paper, we have limited our focus to the IPS equations (as corrected in Johannesson & Parker Reference Johannesson and Parker1989b). The mathematical details of our IPS model are strictly consistent with Johannesson & Parker (Reference Johannesson and Parker1989b) with equations, physical parameters and non-dimensionalizations briefly described in § 2. IPS is ideal for a clean, clear construction of the theory presented here as well as for developing stable numerical methods that require neither smoothing nor interpolation. While the restriction to the IPS models limits application of the present results to rivers in subresonant flow conditions, there are a broad range of naturally occurring rivers for which these conditions apply. Moreover, we find value in comparing subresonant experimental work with this first boundary-specific theory without the complications introduced by consideration of upstream influence in the experiment by exit conditions which may be difficult to fully control or characterize. Nevertheless, more modern models such as JP and ZS (as per the discussion in § 1.1) more fully capture the physics of river dynamics. We believe that subsequent work to extend the contributions of the present paper for IPS to those more complex models would constitute a significant advance for improved modelling and simulation of the full physics of meandering river dynamics. We note that a wealth of experiments for meandering rivers are still required for validation of the many empirical approximations in even our best models, as is discussed by the broad range of scientists contributing to Parker et al. (Reference Parker, Shimizu, Wilkerson, Eke, Abad, Lauer, Paola, Dietrich and Voller2011).
We note that future extensions of the current theory to superresonant conditions and JP or ZS models would constitute a major expansion of the present effort. It appears feasible employing standard mathematical approaches. In the present effort, the IPS model yields a single third-order partial differential equation (PDE) in one dimension. The JP and ZS equations lead to systems of higher-order PDEs. The procedure for analysing such systems follows techniques in the classic monograph Courant & Hilbert (Reference Courant and Hilbert1989), Vol 1. The systems of high-order PDEs are converted to systems of first-order PDEs, and the characteristics of the systems dictate the required boundaries boundary conditions and path to uniqueness results appropriate for subresonant and superresonant systems, respectively. For simulations with combinations of both subresonant/superresonant domains, different PDEs with distinct characteristic equations arise in different regions with standard matching of boundary data required at the interfaces. The present restriction to subresonant conditions and the IPS model makes the specification of proper boundary conditions more straightforward. With this subresonant restriction, the required boundary conditions immediately follow, as shown in the analysis presented in § 3.
We now focus our attention on the contributions of the present paper and impact for IPS computations. In § 3, we find that the boundary conditions required for a unique solution of the initial boundary value include specification of the channel position and slope at the inlet plus one additional boundary condition at the outlet. Notably, a simple specification of outlet position or slope is not sufficient, but rather a specific combination involving slope and some choice of position and curvature is required. This is an important point. It is not just the number of boundary conditions at inlet and outlet that matter – but their specific form. If one chooses from among the proper set for uniqueness, then the results are completely insensitive to boundary condition type and different types will lead to the same results if consistent boundary condition data are provided. A brief example illustrates this concept. If we solved a problem in Stokes flow for pressure driven flow between parallel plates in a finite rectangular domain, uniqueness tells us that specification of velocity on four sides plus the average pressure guarantees uniqueness. We could also use velocity boundary conditions on three sides plus pressure and shear stress at the outlet and find the exact same solution. Similarly, specification of normal velocity and tangent shear stress on two boundaries with velocity and pressure/shear stress on the other boundaries also guarantees the same solution. The key message is that all examples select from a set of boundary condition types which is sufficient for uniqueness. Selecting boundary condition data for problems two or three which are consistent with the solution for the same data from problem one guarantees that all problems give identical answers. This result which eliminates concerns about sensitivity to boundary conditions if the boundary condition set is sufficient is why we prize the standard of sufficiency for uniqueness theory. Generally, selection of proper boundary conditions guarantees this result for the linearized system but also yields excellent results in nonlinear systems. Conversely, using boundary conditions for nonlinear systems which are not even proper in the linearized limit guarantees difficulties in the nonlinear case.
The lack of convergence results and higher-order simulation algorithms for meander dynamics stands in contrast to the results for other fluid dynamics simulations based on standard techniques for PDEs. In our research efforts, our goal is to provide the well formulated initial-boundary-value problems with proper boundary conditions which makes such quality numerical algorithms easily transferable to meandering dynamics simulations. In the present paper, we focus on the formulation of the initial-boundary-value problem, the specification of boundary conditions, the conditions for uniqueness, results for linear stability theory in finite domains and guidance for experimentalists. In a companion paper, we shall show how the proper formulation enables the development of stable, accurate numerical simulations with proper convergence theory and demonstrate convergence results for nonlinear river meander planforms without filtering, smoothing or averaging steps. The numerical algorithm and results in the companion paper are based on the reported efforts in Weiss (Reference Weiss2016).
We turn briefly to the issues to be addressed for finite domains. First with respect to experiments, we note that dynamic forcing at the flow inlet can be either in displacement, as in the van Dijk et al. (Reference van Dijk, Lageweg and Kleinhans2012) experiment, or in entry angle or both. The temporal variation can include time periodic variation, temporal growth, a combination with complex exponent or any preferred stepwise or continuous time dependence amenable to experiments. The well-posed initial-boundary-value problem for the meandering dynamics developed in § 3 allows us to present, in § 4, the appropriate form for linear stability theory in finite domains which predicts the channel dynamics not just in response to a perturbed initial shape, but also to the sustained dynamic boundary perturbations. The linear stability theory results provide valuable guidance to experimentalists for the meander dynamics arising as a function of the form of the upstream driving perturbations. We observe that finite perturbations in inlet angle of entry may prove more accessible to experimental design than the finite amplitudes used in the van Dijk et al. (Reference van Dijk, Lageweg and Kleinhans2012) experiments.
In the final section of the paper § 5, we present comparisons of results from the linear stability theory and field observations. The final results at the end of § 5 give original scaling arguments that show how geological parameters govern the growth of meanders. In particular, we find that the rate of meandering increases with river width, downreach slope and bank erodibility; however, counter to intuition, it decreases with increasing volumetric flow rate.
At the end of this paper, we have provided Appendix A along with a MATLAB script, which provides predictions in dimensional variables based on linear stability theory for both finite and periodic domains.
2. Mathematical formulation
The formulation of the governing equations for a meandering river dynamics in this section follows the classic formulation introduced by Ikeda et al. (Reference Ikeda, Parker and Sawai1981) and with corrections by Johannesson & Parker (Reference Johannesson and Parker1989b). We do not add any new features to this classic set of equations but give a brief resume of the derivation without details to highlight the approximations going into the model equations. The overall governing equations for high Reynolds number flow in meandering rivers are the Navier–Stokes equations together with the continuity equation for incompressible flow. Following the standard procedure, one introduces the Reynolds-averaged equations to model the turbulence with appropriate simple closures for the Reynolds stress. Next, one assumes the length scale for changes in the lateral 
$x, y$ direction is much greater than the length scale for changes in the vertical 
$z$ direction. A depth average equation is obtained by taking averages along the 
$z$-axis, integrating from the bed of the river to the surface of the river. We introduce mean and fluctuating quantities 
$\boldsymbol {u}=\bar {\boldsymbol {v}}+ \boldsymbol {u''}$.
Here, we define the mean velocity vector 
$\bar {\boldsymbol {v}}$, which has two non-zero velocity components (in the 
$x$ and 
$y$ directions). We observe that there is non-zero momentum flux associated with vertical fluctuations in velocity, and define a new stress term 
$\tau _{DA}$ representing the mean momentum flux arising from both the turbulent fluctuations and the new depth-averaged velocity stresses. One must introduce a new closure for 
$\tau _{DA}$ in terms of mean field quantities
\begin{equation} \tau_{DA}=\tau_{RS}+\overline{-\rho \boldsymbol{u}'' \boldsymbol{u}''}. \end{equation}We arrive at the standard shallow-water equations
\begin{equation} \frac{\partial {}}{\partial {t}}(\rho \bar{\boldsymbol{v}})+\text{div}(\rho \bar{\boldsymbol{v}} \bar{\boldsymbol{v}} )={-} \rho {g \boldsymbol{\nabla} \tilde{\xi}}+\frac{\tau_{DA}}{\tilde{h}}, \end{equation}
where 
$\bar {\boldsymbol {v}}$ is the depth-averaged velocity with two components (in the 
$x$ and 
$y$ directions), 
$\tilde {h}$ is local depth, 
$\tilde {\xi }$ is water-surface elevation, 
$\tilde {\eta }$ is bed elevation such that 
$\tilde {h}=\tilde {\xi }-\tilde {\eta }$, 
$\rho$ is fluid density and 
$g$ is acceleration due to gravity.
It is convenient to expand the left-hand side and express the shallow-water equations as
\begin{equation} {\rho \left( \frac{\partial {\bar{\boldsymbol{v}}}}{\partial {t}}+{\bar{\boldsymbol{v}}} \boldsymbol{\cdot} \boldsymbol{\nabla} {\bar{\boldsymbol{v}}} \right) ={-} \rho {g \boldsymbol{\nabla} \tilde{\xi}}+ \frac{\tau_{DA}}{\tilde{h}}}, \end{equation}
and the depth-averaged continuity equation is 
$\tilde {h}\boldsymbol {\nabla }\boldsymbol {\cdot } \bar {\boldsymbol {v}}=0$.
The morphodyamic time scale is much longer than the hydraulic time scale, and so we consider a quasi steady state, which allows us to drop the time derivative from the shallow-water equations
\begin{equation} \rho \left( {\bar{\boldsymbol{v}}} \boldsymbol{\cdot} \boldsymbol{\nabla} {\bar{\boldsymbol{v}}} \right) ={-} \rho {g \boldsymbol{\nabla} \tilde{\xi}}+ \frac{\tau_{DA}}{h}. \end{equation}Following standard practice for the meandering river models, one may choose curvilinear coordinates, based on a curve running along the river centreline and derive (1a)–(c) from Ikeda et al. (Reference Ikeda, Parker and Sawai1981) with corrections from Johannesson & Parker (Reference Johannesson and Parker1989b).
Given the two-dimensional shallow-water equations, the final step is to assume an asymptotic expansion for slender rivers in which the channel half-width 
$b$ is assumed much less than the centreline radius of curvature 
$r_c$. Thus we define 
$\nu ={r_c}/{b}\gg 1$. This approach is due to Ikeda et al. (Reference Ikeda, Parker and Sawai1981) with corrections by Johannesson & Parker (Reference Johannesson and Parker1989b). The process is as follows. Based on the slender-channel assumption, one defines local flow in terms of a reach-averaged tangential velocity (averaged over a streamwise length scale) plus a primed quantity associated with curvature-induced perturbations 
$\tilde {u}=\tilde {U}+u'$. Finally, after considerable algebra, one defines a linearization neglecting nonlinear terms in perturbation variables such as the velocity 
$u'$ and aspect ratio 
$\nu$.
The final form of the momentum governing equations given in the notation of Ikeda et al. (Reference Ikeda, Parker and Sawai1981) is
\begin{equation} { \frac{\partial {u'}}{\partial {\tilde{s}}}+\frac{2}{\tilde{H}}C_fu'= \tilde{b}\left[ -\tilde{U} \frac{\partial {\kappa'}}{\partial {\tilde{s}}} +C_f \kappa' \left( \frac{\tilde{U}^{3}}{gH^{2}}+A\frac{\tilde{U}}{\tilde{H}} - \frac{\tilde{U}}{\tilde{H}}\right) \right].} \end{equation}
Here, 
$C_f$ is an empirical friction coefficient arising from the turbulence closure and 
$A$ is an empirical term known as the scour factor based on the closure for 
$\tau _{DA}$, 
$\kappa (s)$ is the local centreline curvature as a function of arc length introduced when transforming into curvilinear coordinates.
Following the corrections by Johannesson & Parker (Reference Johannesson and Parker1989b), 
$A = A'+A_s$, where 
$A'$ is a measure of the transverse slope and varies between zero (an initial straight channel) and 
$A'_{max}$ (maximum river sinuosity). We use 
$A'$ as defined by Chen & Duan (Reference Chen and Duan2006) and 
$A_s$ as defined by Johannesson & Parker (Reference Johannesson and Parker1989b).
Many authors have made contributions to extend the modelling of meandering rivers to include additional physical phenomena as briefly reviewed in the introduction. We do not consider the most complex or sophisticated models in our present study but focus attention on the classic IPS equations of Ikeda et al. (Reference Ikeda, Parker and Sawai1981) as corrected by Johannesson & Parker (Reference Johannesson and Parker1989b).
With (2.5), we have the required momentum equation for the streamwise perturbation velocity 
$u'$ in the river channel. The remaining modelling approximations concern the rates of deposition and erosion on the banks of the river, which then give an equation for the dynamics of the river centreline. The erosion equation (Ikeda et al. Reference Ikeda, Parker and Sawai1981) derives from the heuristic observation that the magnitude of observed bank erosion is directly proportional to 
$u'$, the perturbation velocity due to curvature where 
$u'$ is equivalent to the difference between the streamwise velocity at the river centreline and the streamwise velocity at the bank. Meandering rivers commonly maintain a constant width, and it is assumed that deposition (as well as the erosion) is proportional to 
$u'$. Let 
$\tilde {x}$ and 
$\tilde {y}$ be the coordinates of the river centreline, and 
$N_{(\tilde {x})}, N_{(\tilde {y})}$ be the 
$\tilde {x}$ and 
$\tilde {y}$ components of a unit vector normal to the river centreline in the plane of the river. Finally, one introduces an empirical erosion coefficient 
$E$ which is a dimensionless coefficient covering a wide range of dependencies based on the composition of the eroding banks, the deposition process, local bank irregularities, etc. With these definitions, the erosion equation is written
\begin{equation} \left[\begin{array}{c} \tilde{x}_{\tilde{t}}\\ \tilde{y}_{\tilde{t}} \end{array} \right] = \left[\begin{array}{c} EN_{(\tilde{x})}u'\\ EN_{(\tilde{y})}u' \end{array} \right]. \end{equation}Note in this equation one must adopt a consistent definition for the sense of the unit normal such that it always points from a designated bank to the opposite bank and does not flip direction. A cross-product of the streamwise unit tangent vector and the unit vector perpendicular to the plane of the meandering river serves this purpose.
2.1. Non-dimensionalization
To facilitate a clear physical understanding of the momentum (2.5) and erosion (2.6) equations, we convert to dimensionless variables, which reduces the number of variables and clarifies fundamental scaling relationships.
We start with the non-dimensionalization of equation (7) in Ikeda et al. (Reference Ikeda, Parker and Sawai1981), which is identical to our (2.5) with all dimensional variables. In this equation, 
$\tilde {U}$ is downstream velocity; 
$A$ is scour factor; 
$\tilde {H}$ is depth; 
$g$ is gravity; 
$\kappa '$ is curvature; 
$C_f$ is a friction coefficient; 
$\tilde {s}$ is arc length; 
$\tilde {b}$ is the river half-width, 
$u'$ is differential velocity. Before considering non-dimensionalization, it proves useful to present the equation
\begin{equation} \frac{\partial {u'}}{\partial {\tilde{s}}}+\underbrace{\left(\frac{2C_f}{\tilde{H}}\right)}_{1/L^{*}}u'= \tilde{b}\left[-\tilde{U} \frac{\partial {\kappa'}}{\partial {\tilde{s}}} +\tilde{U}\underbrace{\left(\frac{2C_f}{\tilde{H}}\right)}_{1/L^{*}} \underbrace{\frac{1}{2}\left( \frac{\tilde{U}^{2}}{g\tilde{H}}+A -1\right)}_{\alpha} \kappa' \right]. \end{equation}
We define the dimensionless parameter 
$\alpha$,
\begin{equation} { \alpha=\frac{1}{2} \left( \frac{\tilde{U}^{2}}{g\tilde{H}}+A - 1 \right)=\frac{1}{2} \left( Fr^{2}+A - 1 \right) .} \end{equation} The secondary flow induces streamwise vorticity, which throws the locus of highest streamwise velocity to the outside bank and enables meander behaviour. In this work, we refer to 
$\alpha$ as the streamwise velocity parameter.
Let 
$H_0$ be a characteristic depth of the river, 
$U_{0}$ be a characteristic streamwise velocity of the river and 
$C_{f0}$ a characteristic friction coefficient. Consistent with standard practice, we define a characteristic length 
$L^{*}$ and characteristic velocity 
$U^{*}$ as
\begin{equation} { L^{*}= \frac{H_0}{2C_{f0}} ; \quad U^{*}={U_{0}} .} \end{equation}Based on a simple momentum balance between fluid inertia and turbulent friction, the characteristic length is the distance downstream required for the velocity profile to recover from disturbances owing to effects of an upstream perturbation or bend.
Throughout this paper, we approximate 
$C_f=C_{f0}$, and hence assume the friction coefficient is constant over the entire river bed. Furthermore, in this paper, and until indicated otherwise (in the nonlinear section of the companion paper), we approximate 
$\tilde {H}=H_0$ and 
$\tilde {U}=U_{0}$, although 
$\tilde {H}$ and 
$\tilde {U}$ are actually functions of sinuosity, 
$\mathcal {S}$.
We also define the following non-dimensional variables: 
$u$ (cross-width velocity perturbation), 
$s$ (arc length), 
$\kappa$ (curvature)
\begin{equation} u=\left(\frac{L^{*}}{\tilde{b}U^{*}}\right){u'} ; \quad s=\frac{\tilde{s}}{L^{*}} ; \quad \kappa={L^{*}}\kappa'. \end{equation}
For reference, we observe the dimensional variables 
$u'$, 
$\tilde {s}$ and 
$\kappa '$ appearing on the right-hand side of these equations, and note these equations are also convenient for returning to dimensional variables.
For consistency with the notation of Ikeda et al. (Reference Ikeda, Parker and Sawai1981), we will use primes and tildes as dimensional quantities and unmarked lowercase and uppercase variables as non-dimensional quantities. (Earlier dimensional variables in this chapter, including 
$\rho$ and 
$g$, did not use primes or tildes, but will not appear explicitly henceforth.) In particular, primes indicate cross-stream perturbations; tildes indicate other dimensional quantities; uppercase letters indicate values that have been cross-stream averaged. Lowercase, unmarked variables are all non-dimensional. These definitions, when applied to (2.7), give the following dimensionless momentum equation:
\begin{equation} { \frac{\partial {u}}{\partial {s}}+u={-}\frac{\partial {\kappa}}{\partial {s}} +{\alpha}\kappa }, \end{equation}
with the dimensionless parameter 
$\alpha$ given by (2.8) above.
Next, we address the non-dimensionalization of the erosion equation (2.6).
We define the non-dimensional variables 
$x$, 
$y$, 
$u$ and 
$t$ such that
\begin{equation} x=\tilde{x}/L^{*}; \quad y=\tilde{y}/L^{*};\quad u=\left(\frac{L^{*}}{\tilde{U}}\right)\frac{u'}{\tilde{b}};\quad t= \tilde{t}/t^{*}, \end{equation}
where the characteristic time scale, 
$t^{*},$ is taken to be
\begin{equation} { t^{*}=\frac{L^{*}}{U^{*}}\frac{1}{E (\tilde{b}/L^{*})} .} \end{equation}
Recall that 
${E}$ is a dimensionless proportionality constant between the normal bank migration rate and the excess near-bank downstream fluid velocity 
$u'$ (Ikeda et al. Reference Ikeda, Parker and Sawai1981). Thus the time scale is based on the transverse erosion velocity of the river centerline.
We write the normal vector in terms of arcwise derivatives of 
$x$ and 
$y$, yielding
\begin{equation} \left[\begin{array}{c} N_{({{x}})}\\ N_{({{y}})} \end{array} \right] =\left[ \begin{array}{c} -{y}_{{s}}\\ {x}_{{s}} \end{array} \right]. \end{equation}Altogether, these definitions give us the following non-dimensionalized erosion equation:
\begin{equation} { \left[ \begin{array}{c} \dfrac{\partial {x}}{\partial {t}}\\ \dfrac{\partial {y}}{\partial {t}} \end{array} \right] = \left[ \begin{array}{c} -y_su\\ x_su \end{array} \right] .} \end{equation}3. Mathematical theory
The governing equations for meandering rivers are
\begin{gather} \frac{\partial {u}}{\partial {s}}+u={-}\frac{\partial {\kappa}}{\partial {s}} +{\alpha}\kappa \end{gather}
\begin{gather} \left[\begin{array}{c} x_t\\ y_t \end{array} \right]= \left[ \begin{array}{c} -y_su\\ x_su \end{array} \right] \end{gather}where
\begin{equation} \alpha =\tfrac{1}{2} \left( Fr^{2}+A - 1 \right). \end{equation}3.1. Linearization
The nonlinearities in these equations enter through the curvature and derivative of curvature terms
\begin{gather} \kappa=y_sx_{ss}-x_sy_{ss}, \end{gather}
\begin{gather}\kappa_s=y_sx_{sss}-x_{s}y_{sss}. \end{gather}
To simplify the analysis of the basic mathematical character of the equations, we consider a linearized form for small amplitude perturbations 
$y$ measured from a straight line which we choose along the 
$x$ axis. In the limit as amplitude 
$\to 0$: 
$s \to x$, 
$x_s \to 1$, 
$x_{ss}$, 
$x_{sss} \to 0$ and 
$y_s$, 
$y_{ss}$, 
$y_{sss} \to y_x$, 
$y_{xx}$, 
$y_{xxx}$. With these definitions, 
$\kappa =-y_{xx}$ and 
$\kappa _s=-y_{xxx}$.
Our linearized equations become
\begin{gather} u_x+u= y_{xxx}-\alpha y_{xx}, \end{gather}
\begin{gather}y_t=u. \end{gather}Combining (3.3) and (3.4), we obtain this nonlinear third-order PDE, which was described in Ikeda et al. (Reference Ikeda, Parker and Sawai1981)
\begin{equation} { y_{tx}+y_t=y_{xxx}-\alpha y_{xx} .} \end{equation}In the previous section, we considered the question of appropriate initial and boundary conditions for (2.11) and (2.15). With the linearization (3.5), we find that the problem is amenable to a direct proof to determine the conditions required for uniqueness.
3.2. Problem definition
Consider the third-order PDE (3.5) on the interval 
${0} \leqslant x \leqslant \ell$ with 
$\alpha > 0$ constant.
We define the following:
\begin{align} \text{PDE: }\ &y_{tx}+y_t=y_{xxx}-\alpha y_{xx} \quad\text{over}\ {0} \leqslant x \leqslant \ell \end{align}
\begin{align} \text{IC: }\ & y = f(x)\text{ at }t=0 \end{align}
\begin{align} \text{BC: }\ & {m}_{11} y +{m}_{12} y_x+{m}_{13} y_{xx}=g_1(t) \text{ at }x=x_1 \end{align}
\begin{align} & {m}_{21} y +{m}_{22} y_x+{m}_{23} y_{xx}=g_2(t) \text{ at }x=x_2 \end{align}
\begin{align} & {m}_{31} y +{m}_{32} y_x+{m}_{33} y_{xx}=g_3(t) \text{ at }x=x_3. \end{align}Here, we have considered a system of three linear boundary conditions with no a priori constraints on the boundary condition locations.
For any two solutions which satisfy the initial and boundary conditions, the difference between the two solutions is zero over 
$x$ at the initial time and over all 
$t$ at the boundary condition positions. For uniqueness, it is sufficient to demonstrate that these homogeneous conditions guarantee 
$y \equiv 0$. We now determine the conditions on the 
$m_{ij}$ coefficients, as well as the locations 
$x_i$ required for this condition.
For convenience, we define 
$w=y_x+y$ and note that this leads to 
$y_x = w-y$, 
$y_{xx} = w_x-y_x$, 
$y_{xxx}= w_{xx}-y_{xx}$.
The governing equation (3.5) can now be written as
\begin{equation} w_t=w_{xx}-(1+\alpha)w_x+(1+\alpha)w-(1+\alpha)y. \end{equation} We start by multiplying both sides by 
$w$
\begin{equation} ww_t=ww_{xx}-(1+\alpha)w w_x+(1+\alpha)w^{2}-(1+\alpha)wy \end{equation}or
\begin{equation} \frac{\partial {}}{\partial {t}} \left( \frac{1}{2} w^{2} \right)=\frac{\partial {}}{\partial {x}}(w w_x)-(w_x^{2})-(1+\alpha) \frac{\partial {}}{\partial {x}}\left( \frac{1}{2} w^{2} \right)+(1+\alpha) w^{2}-(1+\alpha) wy. \end{equation}
We substitute the last term of the right-hand side with the expression for 
$w$ in terms of 
$y$ and write
\begin{equation} (1+\alpha) wy =(1+\alpha) (y_x+y)y=(1+\alpha)\frac{\partial {}}{\partial {x}}\left(\frac{1}{2} y^{2} \right)+(1+\alpha)y^{2}. \end{equation}Substituting this result into (3.12) and rearranging the terms
\begin{align} {\frac{\partial {}}{\partial {t}} \left( \frac{1}{2} w^{2} \right)}&=\frac{\partial {}}{\partial {x}}(w w_x)-(1+\alpha) \frac{\partial {}}{\partial {x}}\left( \frac{1}{2} w^{2} \right)-(1+\alpha) \frac{\partial {}}{\partial {x}}\left( \frac{1}{2} y^{2} \right)-w_x^{2}\nonumber\\ &\quad -(1+\alpha) y^{2}+(1+\alpha) w^{2}, \end{align}
and we define 
$Q$ as 
$1/2$ of the 
$L_2$ norm of 
$w$
\begin{equation} Q(t)=\underbrace{\int^{\ell} _{0} \left( \tfrac{1}{2} w^{2} \right) {{\rm d}\kern0.06em x}}_{\geqslant 0}. \end{equation}
The value of 
$Q(t)$ must be 
$\geqslant 0$ by simple examination of the form of (3.15). From this point forward, our goal is to demonstrate that 
$Q\equiv 0$, which would imply that 
$w \equiv 0$, and, with consideration of the boundary, that 
$y\equiv 0$, proving uniqueness. To do so, we only need to show that 
$Q \leqslant 0$.
We integrate (3.14) with respect to 
$x$ over (0, 
$\ell$) to obtain
\begin{align} \frac{\partial {Q}}{\partial {t}}&=(ww_x)_{0} ^{\ell} -(1+\alpha)\left( \frac{1}{2} w^{2} \right)^{\ell} _{0} -(1+\alpha)\left( \frac{1}{2} y^{2} \right)^{\ell} _{0} -\int_{0} ^{\ell} (w_x^{2}+(1+\alpha)y^{2})\,{{\rm d}\kern0.06em x}\nonumber\\ &\quad +2(1+\alpha)Q. \end{align}Equation (3.16) is of the form
\begin{equation} \frac{\partial {Q}}{\partial {t}}=g(t)+2(1+\alpha)Q, \end{equation}where
\begin{equation} g(t)=(ww_x)_{0} ^{\ell} -(1+\alpha)\left( \frac{1}{2} w^{2} \right)^{\ell} _{0} -(1+\alpha)\left( \frac{1}{2} y^{2} \right)^{\ell} _{0} -\int_{0} ^{\ell} (w_x^{2}+(1+\alpha)y^{2})\,{{\rm d}\kern0.06em x}. \end{equation}The solution for the (3.17) is
\begin{equation} Q(t)=Q_0\exp({2(1+\alpha)t})+\int_{0} ^{t} g(t) \,{\rm d}t. \end{equation} We know by the homogeneous initial condition that 
$Q_0 = 0$. It remains only to show that 
$g \leqslant 0.$ The integral term in (3.18) was already shown to be non-positive. Under the appropriate conditions, the first three terms in 
$g(t)$ will be shown to be non-positive as well, in which case 
$Q\equiv 0$ and so 
$w \equiv 0$. This implies that 
$\partial u/\partial x+u\equiv 0$, which implies that 
$u=u_0{\rm e}^{-x}$ and 
$u_0=0$ from the boundary condition at 
$x=a$.
Thus we have the desired result: 
$y\equiv 0$ and the solution is unique.
3.3. Homogenous boundary conditions
The final step is to determine boundary conditions which make the terms
\begin{equation} (ww_x)_{0} ^{\ell} -(1+\alpha)\left( \tfrac{1}{2} w^{2} \right)^{\ell} _{0} -(1+\alpha)\left( \tfrac{1}{2} y^{2} \right)^{\ell} _{0} \end{equation}strictly less than zero.
There are two cases of interest.
(i) For periodic boundary conditions in
$y$, all terms are identically zero and we have uniqueness.(ii) For uncoupled boundary conditions at
${0}$ and $\ell$, we require the sum of terms at each end must independently be non-positive.
For convenience, let 
$\zeta =(1+\alpha )/2$, and rewrite (3.20) accordingly.
We then require
\begin{gather} ww_x-\zeta w^{2}-\zeta y^{2} \geqslant 0\quad \text{at }x={0} , \end{gather}
\begin{gather}ww_x-\zeta w^{2}-\zeta y^{2} \leqslant 0\quad \text{at }x=\ell. \end{gather}3.4. Sufficient boundary conditions
For a simple case of sufficient boundary conditions, choose two upstream boundary conditions at 
$x={0}$ specifying the slope and position of the curve
\begin{equation} w=0,\quad y=0\quad \text{at }x={0}, \end{equation}
and a downstream boundary condition at 
$x=b$ of the form
\begin{equation} w_x+hw=0\quad\text{at }x=\ell, \end{equation}
where 
$h>-\zeta$.
Equation (3.22) becomes
\begin{equation} -hw^{2}-\zeta w^{2} -\zeta y^{2}\leqslant 0, \end{equation}which is always satisfied. Equation (3.24) is the same as
\begin{equation} \frac{\partial {}}{\partial {x}}\left(y+y_x\right)+h(y+y_x)=0. \end{equation}
Therefore, the downstream boundary condition equation (3.22) may be written in terms of 
$y$ as
\begin{equation} y_{xx}+(1+h)y_x+hy=0. \end{equation}
Note the downstream conditions 
$w_x = 0$ and 
$w=0$ also guarantee uniqueness. These correspond to the limits 
$h=0$ and 
$h\to \infty$, which give simpler downstream boundary conditions which also guarantee uniqueness
\begin{equation} y_x+y=0\quad or \quad y_{xx} + y_x=0. \end{equation}3.5. General boundary conditions
We define a vector 
$v=(y,w,w_x)$ and rewrite the expression on the left-hand side of (3.21) and (3.22) as a quadratic form 
$\psi$ in terms of 
$v$
\begin{equation} \psi=m_{ij}v_iv_j, \end{equation}and
\begin{equation} M=\left( \begin{array}{ccc} -\zeta & 0 & 0 \\ 0 & -\zeta & 1/2\\ 0 & 1/2 & 0 \end{array} \right). \end{equation}
The desired requirements are now 
$\psi \geqslant 0$ at 
$x={0}$ and 
$\psi \leqslant 0$ at 
$x=\ell$.
We note that 
$M$ is real and symmetric; hence its eigenvalues are real. In particular, 
$\lambda _1 =-\zeta$, 
$\lambda _2 = -\zeta /2-d$ and 
$\lambda _3 =-\zeta /2+d$ where 
$d=\frac {1}{2}\sqrt {\zeta ^{2}+1/4}$
The associated eigenvectors are
\begin{gather} {u}_1=(1, 0, 0), \end{gather}
\begin{gather}{u}_2=(0, \lambda_2, 1/2)/r_2, \end{gather}
\begin{gather}{u}_3=(0, \lambda_3, 1/2)/r_3, \end{gather}
where 
$r_2 = \sqrt {\lambda _2^{2}+1/4}$ and 
$r_3=\sqrt {\lambda _3^{2}+1/4}$.
An arbitrary vector 
$v$ may now be written in terms of the eigenvectors as
\begin{equation} v=\phi_1 u_1 +\phi_2 u_2 +\phi_3 u_3, \end{equation}
where the scalar coefficients 
$\phi$ are functions of 
$y$ and its derivatives. Taking the inner product and noting orthogonality,
\begin{gather} \phi_1=v_1 \cdot u_1 = y, \end{gather}
\begin{gather}\phi_2=v_2 \cdot u_2 =\left[ \lambda_2 w +\tfrac{1}{2}w_x \right]/r_2, \end{gather}
\begin{gather}\phi_3=v_3 \cdot u_3 =\left[ \lambda_3 w +\tfrac{1}{2}w_x \right]/r_3. \end{gather}From (3.29)
\begin{equation} \psi=(M v ) \cdot v, \end{equation}we find that
\begin{gather} Mv=\phi_1 M u_1 +\phi_2 M u_2 + \phi_3 M u_3, \end{gather}
\begin{gather}Mv=\lambda_1 \phi_1 u_1 +\lambda_2 \phi_2 u_2 + \lambda_3 \phi_3 u_3. \end{gather} We may now write 
$\psi$ in terms of the eigenvalues 
$\lambda$ and coefficients 
$\phi$ as
\begin{equation} \psi=(M v ) \cdot v =\lambda_1\phi_1^{2}+\lambda_2 \phi_2^{2}+\lambda_3 \phi_3^{2}. \end{equation} Note that the first two eigenvalues are negative while the third is positive. Setting 
$\psi =0$ gives the surface of an elliptic cone (figure 1) in 
$\phi$ space with the cone along the 
$\phi _3$ axis.
Figure 1. Elliptic cone.
To satisfy the boundary conditions at 
$x=b$, 
$\phi <0$, we may choose any plane which intersects the origin and lies completely outside the cone. Let the plane be described by director 
$d$ which is perpendicular to the chosen plane
\begin{equation} d_1\phi_1+d_2\phi_2+d_3\phi_3=0. \end{equation}
Substituting for the value of 
$\phi$ gives the boundary condition for 
$(y,w,w_x)$ or, equivalently, for 
$(y,y_x,y_{xx})$ at 
$x=b$.
Next, for the boundary condition at 
$x={0}$, we require that 
$\psi \geqslant 0$, which requires a vector 
$e$ which lies within the cone. Note that we know that the vector 
$e$ must have a non-zero component 
$e_3$. Without loss of generality, we define two planes whose intersection is the line coincident with 
$e$ whose perpendicular vectors are respectively
\begin{equation} (0, e_3, -e_2)\quad\text{and}\quad (e_3,0,-e_1). \end{equation}
It is trivial to show that each of these vectors is perpendicular to 
$e$ and hence that 
$e$ lies in both planes.
We then specify the two required boundary conditions at 
$x={0}$ as
\begin{equation} e_3\phi_2-e_2\phi_3=0, \end{equation}and
\begin{equation} e_3 \phi_1-e_1 \phi_3=0. \end{equation} The coefficients of 
$e$ are unrestricted as long as the vector lies within the elliptical cone described above.
3.6. Summary
In summary, the theory determines that there are three sets of appropriate boundary conditions. All of these have various levels of restrictiveness, and are sufficient for uniqueness.
(i) Periodic boundary conditions.
(ii) Simple boundary conditions.
(3.46)\begin{gather} y=f_1\quad \text{at}\ x=0 \end{gather}(3.47)\begin{gather}y_x=f_2\quad \text{at}\ x=0 \end{gather}(3.48)\begin{gather}y_{xx}+(1+h)y_x+hy=f_3 \quad \text{at} \ x=L. \end{gather}(iii) General boundary conditions.
As described in the previous section where
(3.49)\begin{gather} \lambda_2={-}\zeta/2-d \end{gather}(3.50)and\begin{gather}\lambda_3={-}\zeta/2+d, \end{gather}(3.51)with\begin{equation} d=\tfrac{1}{2}\sqrt{\zeta^{2}+1/4}, \end{equation}(3.52)\begin{equation} \zeta=\frac{1+\alpha}{2}. \end{equation}
The key result of this uniqueness proof is that, except in the case of periodic boundary conditions, three boundary conditions must be specified. In a finite domain, the boundary condition will affect the behaviour of the system. We shall find in later sections that the downstream boundary has minimal effect on the behaviour of the channel. The two upstream boundary conditions, however, such as the displacement and slope at 
$x={0}$, have a significant effect.
4. Linear stability theory
In this section, we begin with a brief reprisal of the linear stability theory for periodic domains given by Ikeda et al. (Reference Ikeda, Parker and Sawai1981). These results will prove useful for comparison with our results for finite non-periodic domains. While the finite-domain linear theory is straightforward, it could not be developed previously owing to lack of knowledge of appropriate boundary conditions.
For the general case we follow standard practice and consider solutions of the form
\begin{equation} y={B }{{\rm e}^{\gamma t}{\rm e}^{rx}}, \end{equation}
where 
$\gamma$ and 
$r$ may be complex coefficients. For the spatially periodic case, we will restrict ourselves to coefficients 
$r={\rm i}k$ which are strictly imaginary.
Substitution into the governing equation 
$y_{tx}+y_t=y_{xxx}-\alpha y_{xx}$ yields an algebraic equation relating 
$\gamma$ and 
$r$
\begin{equation} r^{3}-\alpha r^{2}-\gamma r-\gamma=0. \end{equation} We solve for 
$\gamma$
\begin{equation} {\gamma=\frac{-{\rm i}k^{3}+\alpha k^{2}}{1+{\rm i}k}.} \end{equation} We can decompose 
$\gamma$ into its real (
$\gamma _R$) and imaginary (
$\gamma _i$) parts
\begin{equation} \gamma_{R}=\frac{\alpha-k^{2}}{1+k^{2}}k^{2} ,\quad \gamma_{i}={-}\frac{\alpha+1}{1+k^{2}}{\rm i}k^{3}. \end{equation} The function 
$y$ takes the form
\begin{equation} y={B }\exp({\gamma_R t})\exp({{\rm i}k(x-ct)}) \quad where \ c={-}\gamma_i /k, \end{equation}
which gives a wave moving downstream with wave speed 
$c$ and temporally growing amplitude with growth rate 
$\gamma _R$.
It follows that long waves are unstable; the condition for stability with 
$\gamma _{R}< 0$ is that the wavenumber is 
$k > k_{neutral}$ where
\begin{equation} k_{neutral}=\sqrt{\alpha}, \end{equation}while the condition for maximum growth rate is
\begin{equation} k_{max}=\sqrt{-1+\sqrt{1+\alpha}} .\end{equation} From a physical viewpoint, the rate of change of curvature 
$\kappa _s,$ which is represented by the 
$y_{xxx}$ term in the linearized equation, is stabilizing whereas the curvature 
$\kappa$ associated with the 
$y_{xx}$ term is destabilizing.
For a value 
$\alpha =5$, typical for geophysical systems, the growth rate and wave speed as functions of wavenumber are plotted in figures 2 and 3.
Figure 2. Temporal meander amplitude growth rate as a function of wavenumber. For dimensional units, 
$\gamma _R / t^{*}$ will yield units of 1/time; 
$k/L^{*}$ will yield units 1/length.
Figure 3. Wave speed (
$c = \gamma _i/k$) as a function of wavenumber. Asymptote at 
$\alpha +1$. For dimensional units, 
$c L^{*}/t^{*}$ will yield units length/time; 
$k/L^{*}$ will yield units 1/length.
Next we consider a generalization of the linear stability theory. For spatially non-periodic solutions both 
$r$ and 
$\gamma$ may be complex.
Substituting the functional form (4.1) into the PDE yields a cubic equation for the roots 
$r$ and we designate three specific roots as 
$r_1,\ r_2,\text { and }r_3$ which satisfy the cubic equation
\begin{equation} r^{3}-\alpha r^{2}-\gamma r-\gamma=0. \end{equation} We introduce non-periodic solutions with complex values for the exponential coefficients 
$r$. We write the perturbation solution as
\begin{equation} y={\rm e}^{\gamma t}({B }_1{\rm e}^{r_1x}+{B }_2{\rm e}^{r_2x}+{B }_3{\rm e}^{r_3x}). \end{equation}
For a given complex growth coefficient 
$\gamma$ there are three complex roots 
$r_1, r_2,\ r_3$; the full mapping of these roots for the entire complex plane of 
$\gamma$ is straightforward (via numerical solution) but not very illuminating or instructive for physical realizations. We focus here on two special cases: (i) solutions where at least one root, say 
$r_1 = {\rm i}k$ is strictly imaginary such that the spatially periodic solutions fall within the family and (ii) solutions where the growth rate is purely imaginary, say 
$\gamma = {\rm i} \, \omega$. This second class of solutions is of interest in driving system boundaries with a temporally periodic disturbance and observing spatial growth of disturbances. Taking up class (i) first – with specified wavenumber 
$(r_1)_i =k$, the complex growth rate 
$\gamma$ is known from (4.3) in periodic theory. The factor 
$r- r_1$ may be divided out of the cubic polynomial and the resultant quadratic equation solved directly for the roots 
$r_2$ and 
$r_3$. For convenience we note the coefficients in the quadratic equation
\begin{gather} r^{2}+c_1r+c_0=0, \end{gather}
\begin{gather}c_0={-}\gamma+r_1^{2}-\alpha r_1,\quad c_1=r_1-\alpha. \end{gather} We find three solutions. The first associated with 
$r_1={\rm i}k$ is our strictly periodic solution; the second associated with 
$r_2$ having a negative real part decays exponentially from the inlet over length scales of the order of a wavelength 
$2 {\rm \pi}/k$; the third 
$r_3$ affecting the outlet has a positive real part. While this nominally implies exponential spatial growth downstream, the imposition of the downstream boundary condition shows that this third solution is significant only in a small region at the system outlet and has negligible effects upstream for any plausible values of the parameter 
$\alpha$. In effect, it is rapid exponential decay in the upstream direction. The numerical coefficients 
${B }_1$, 
${B }_2$ and 
${B }_3$ will be determined given the exact form of the non-homogeneous boundary conditions. As noted 
${B }_3$ is largely irrelevant given the rapid decay upstream. By choosing a driven amplitude or slope at the inlet to a flume with proper complex 
$\gamma$ one could continuously create a train of waves which are spatially periodic away from the inlet and outlet. Unfortunately, given the rapid propagation velocity, it will be difficult to observe significant growth of the periodic wave forms within the length of channels of practical dimension.
The next class of interesting solutions is that of temporally periodic solutions with temporal coefficient 
$\gamma ={\rm i}\omega$. Numerical values for the three roots for the spatial coefficients 
$r_1$, 
$r_2$, 
$r_3$ as a function of 
$\omega$ are shown in figures 4 and 5. For convenience, the real and imaginary parts are plotted separately.
Figure 4. Spatial growth per unit length: real part (Re) of roots 
$r_1$, 
$r_2$ and 
$r_3$ (as per (4.8)) for 
$\alpha =5$ and as a function of driving frequency 
$\omega$. For dimensional units, 
$R(r)/L^{*}$ will yield units of 1/length; 
$\omega /t^{*}$ will yield units 1/time.
Figure 5. Emerging wavenumbers: imaginary part (Im) of roots 
$r_1$, 
$r_2$ and 
$r_3$ (as per (4.8)) for 
$\alpha =5$ and as a function of driving frequency 
$\omega$. For dimensional units, 
${\rm Im}(r)/L^{*}$ will yield units of 1/length; 
$\omega /t^{*}$ will yield units 1/time.
The functional form for 
$y$ for each of the roots 
$r$ gives a function
\begin{equation} y={B }\exp({r_R x})\exp({{\rm i} r_i (x-ct)}) \quad where \ c={-}\omega/r_i, \end{equation}
which gives a wave moving downstream with wave speed 
$c$ and spatially growing or decaying amplitude with growth rate 
$r_R$.
Of the three roots, the root 
$r_3$ always has large positive real part. This is the part of the solution associated with the influence at the downstream boundary and has negligible affect away from that region. The root 
$r_2$ has negative real part for all 
$\omega > 0$ and represents the decay of the solution away from the upstream boundary. Note that the real part of 
$r_3$ is much larger than that of 
$r_2$, implying upstream influence at the downstream boundary is much more localized. The root 
$r_1$ is the most interesting for the stability analysis. For high frequencies above 
$\omega =11.2$, the root has negative real part with downstream decay. For high frequencies, all disturbances quickly decay and no significant downstream propagation of the boundary perturbation is observed. For low frequencies 
$\omega <11.2$, 
$r_1$ had positive real part. The boundary perturbations show exponential spatial growth as they propagate downstream. An expanded plot of the real part of 
$r_1$ is shown in figure 6. From this plot, one can easily see the frequencies for neutral stability and for the frequency for maximum growth rate. The corresponding wave speeds, which depend on the imaginary part of 
$r_1$, are shown in figure 7.
Figure 6. Spatial growth per unit length (real part of 
$r_1$) as a function of driving frequency 
$\omega$ at 
$\alpha =5$. For dimensional units, 
$r_1/L^{*}$ will yield units of 1/length; 
$\omega /t^{*}$ will yield units 1/time.
Figure 7. Wave speed (c) as a function of driving frequency 
$\omega$ (
$c=-\omega /{\rm Im}(r_1$)) at 
$\alpha =5$. For dimensional units, 
$cL^{*}/t^{*}$ will yield units of length/time; 
$\omega /t^{*}$ will yield units 1/time.
The actual values of the coefficients 
${B }_1$ and 
${B }_2$ and the form of the solution in the vicinity on the inlet will depend on the type of the boundary forcing imposed. For illustration purposes, we include examples for zero displacement and oscillatory slope in figure 8, and for an oscillatory displacement and fixed zero entry slope in figure 9. In each case, we show a sub-critical frequency with spatial growth and a super-critical frequency with spatial decay. In all cases, the frequency 
$\omega$ is chosen near 
$11.2$ to limit the spatial growth per wavelength and allow easier visualization. Higher and lower frequencies would show drastic decay and growth respectively. We note that the oscillatory slope and oscillatory displacement give similar planforms, except very close to the inlet. We further note that, in a real physical system, the oscillatory slope could be achieved with an inlet that rotates in angle around the point of entry, i.e. in two dimensions, so that it always remains in the same plane as the flume. Alternatively, a rotating paddle in the inlet channel before the test section can give alternating maximal flow rates on different sides of the channel which will naturally lead to rotating channel inlet angles in the erodible bed.
Figure 8. Solution to the linearized meander-morphodynamics equation (4.8) where 
$\gamma = {\rm i} \omega$; 
${B }_1=1$, 
${B }_2=-1$ and 
${B }_3=0$ such that angle of entry is driven and displacement is fixed. (a) Spatial growth (
$\omega =10.7$ and 
$\alpha =5$). (b) Spatial decay (
$\omega =12.1$ and 
$\alpha =5$).
Figure 9. Solution to the linearized meander-morphodynamics equation (4.8) where 
$\gamma = {\rm i} \omega$; 
${B }_1=1$, 
${B }_2=-r_1/r_2$ and 
${B }_3=0$ such that displacement is driven and angle of entry is fixed. (a) Spatial growth (
$\omega =10.7$ and 
$\alpha =5$). (b) Spatial decay (
$\omega =12.1$ and 
$\alpha =5$).
The wave speed for oscillatory boundary perturbations shown in figure 7 provides an interesting comparison with the plot of spatial growth rate shown in figure 6. The earlier figure shows a maximum spatial growth around 
$\omega =3.1$ while the critical frequency for neutral growth is at 
$11.2$. The wave speed is a monotonically increasing function of 
$\omega$ in this range, and the wave speed appears to increase without limit at increasing 
$\omega$ – waves are highly dispersive. Given these speeds, the first wave one might observe far downstream might not be associated with the frequency for maximum spatial growth rate, but rather depend on the distribution of initial amplitudes for the slowly growing but fastest propagating modes. The growth in asymptotic wave speed for large 
$\omega$ is slower than linear, however, that high frequency behaviour is largely irrelevant. For all frequencies above 
$\omega =11.2$, the perturbations are exponentially decaying with spatial distance and wave speeds above that frequency are academic. Moreover, the behaviour as a function of 
$\omega$ is dependent on the imposed frequency at the driven boundary, and hence high frequencies of oscillation are of little interest and would rarely be utilized. One should be cautious in drawing too many conclusions from this analysis. The general equations are highly nonlinear and in a companion paper, we shall see how the nonlinearity rapidly modifies the downstream evolution from the results of linear theory.
In practice, these eigenvalues are extremely poorly conditioned, and so we rely on numerical methods in the examples that follow.
4.1. Transient solution
In the section above, we focused on the important results from linear stability theory – showing both spatially periodic/temporally growing and temporally periodic/spatially growing solutions. Each of the solutions we discussed could be expressed in the simple form of exponential functions with complex valued coefficients. To study further the case of driven boundaries, and to compare with direct experimental observations, one would have to consider the transient solution based on a mismatch between the initial condition and the periodic wave forms. An example would be an oscillatory driven inlet slope or displacement in a flume with an initially straight river channel. For linear PDEs, one would traditionally break the solution into two parts. One part as discussed in the previous section would cover the oscillatory driving, while the second part based on eigenfunction expansions would address the transient solutions. Unfortunately, for the present third-order PDE, the numerical solution for the eigenvalues is extremely poorly conditioned and owing to the convective and dispersive nature of the solutions, numerical solutions for the eigenfunctions and appropriate eigenfunction expansions are also poorly conditioned. As a consequence, it is much more effective to study the transient effects by direct numerical solutions. In the companion to this paper, we will be presenting a robust numerical procedure for solving the full nonlinear equations and we will leave a discussion of numerical algorithms for that paper. We will observe that numerical solutions of the linearized equations can be obtained from those algorithms or from a simplified approach for the linear PDEs. A full discussion of transient effects in the nonlinear regime will be included in the companion paper. However, it is important to give some transient results for the linear regime in this section as they introduce some important features into the dynamics.
Figure 10 shows an example of the effect of the transient solution. If we consider an oscillatory slope at the inlet at high frequency, one obtains a spatially decaying oscillatory wave form essentially as seen in figure 10(b). If the initial condition does not agree with the oscillatory function evaluated at 
$t=0$, the initial condition produces a transient solution and the total solution at early instants of time is dramatically different from the purely oscillatory part as seen in figure 10(a). After several periods of oscillation, the transient part propagates out of the domain, leaving behind only the oscillatory part as seen in figure 10(b). In the linear regime, the transient and oscillatory solutions do not interact, and what is left behind is the purely oscillatory solution unmodified by the passage of the transient solution. This is true independent of the value of 
$\omega$. The interaction of transient and driven boundary solutions becomes of major importance in the nonlinear regime.
Figure 10. Driving frequency 
$\omega =12.1$ and 
$\alpha =5$. (a) Transient and oscillatory solution at time = 6 periods. (b) Oscillatory solution after transient solution has propagated out of the system at 
${\rm time}=60$ periods.
While figure 10 showed the comparison of the oscillatory and transient solutions in a single frame with superposition over multiple instants of time, we next consider the comparison over several individual time instants to illustrate both the growth and propagation of the transient solution.
In figure 11, we start with an initially straight channel; in effect the initial condition for the transient solution is the negative of oscillatory function evaluated at 
$t=0$; where the oscillatory part has upstream boundary oscillation at a frequency of 
$\omega =10.3$. In this plot, we see how the transient and oscillatory solutions evolve over 70 periods of oscillation. Overall, the wavelengths observed from the oscillatory driven boundary will be those resulting from the imaginary parts of the roots of the cubic equation of last section for the specified value of 
$\omega$. The perceived wavelengths associated with the transient solution will formally be determined by the form of the eigenfunctions of the boundary-value problem for the PDE. For relatively long domains with many periods of the wave form, such as presented here, these wavelengths will be close to those predicted by spatially periodic linear theory. Note that the transient solution which propagates out of the domain takes the form of a wave packet with a local position of large amplitude disturbance with amplitude growing exponentially in time. In the nonlinear model, large amplitudes will lead to drastic modification of the waveform, as has been demonstrated in numerous previous studies. These effects are studied in detail in the numerical studies in our subsequent paper.
Figure 11. Propagation at spatially growing frequency 
$\omega =10.3$. Panels (a) and (b) are identical except that the 
$y$-axis range is fixed on the left and allowed to drift on the right. Therefore we can see the oscillatory solution clearly on the left side and the transient solution clearly on the right.
4.2. Entrance effects
In the preceding sections, we have discussed three broad classes of solution: spatially periodic solutions, temporally periodic solutions and transient solutions which arise owing to the form of the initial conditions. In this section, we will consider how these different solutions will affect hypothetical experiments in flumes of finite length. We begin with an example to illustrate the danger of using classic spatially periodic theory in infer the behaviour in a finite flume where entrance effects may be significant. In figure 12, we consider the case of a perturbed initial condition corresponding to a small amplitude sine wave. In this run, we set 
$\alpha =5$ and a wavenumber of 
$k=1.5$, which corresponds to spatial growth (although not as fast as at the maximum growth 
$k=1.2$) and downstream wave propagation (faster than at 
$k=1.2$). We wish to illustrate the time-dependent growth which a true spatially periodic domain would produce with this initial condition. To mimic those conditions, we merely need to drive the upstream end of the channel with the appropriate time-dependent angle consistent with a periodic wave propagating from upstream. When we drive the upstream boundary with the proper condition from spatially periodic theory, we see significant temporal growth in amplitude exactly in accordance with classic stability theory. Next, we consider the exact same initial condition, but we use an entry boundary condition where we clamp the upstream boundaries such that 
$y(x=0,t)=0$ and 
${{\rm d} y}/{{\rm d}\kern0.06em x}(x=0,t)=0$. This reflects the true set-up in an experimental flume with fixed inlet. In this case, the portion of the wave far from the inlet (perhaps two wavelengths downstream) appears unchanged but the clamped inlet yields a nearly flat line propagating downstream. This plot shows that one cannot rely on predictions from spatially periodic theory to estimate observations in an experimental flume. One would need a significantly longer channel and will only be able to observe for short time periods before the wave exits the flume. As a point of comparison, for an experiment with a reasonable channel depth of 
$\tilde {H}=0.015\ {\rm m}$ and friction coefficient 
$C_f=0.015\ {\rm m}$, the characteristic length would be 
$L^{*}=0.5\ {\rm m}$ such that figures 12(a) and 12(b) would represent flumes of a 10 m length.
Figure 12. Entrance effects for a system at 
$\alpha =5$ and initial perturbation at a wavenumber of 
$k=1.5$. Thick black line indicates most recent time step. Panel (a) demonstrates periodic upstream boundary conditions, where the inlet location is fixed, but angle of entry is oscillating. Panel (b) demonstrates homogeneous upstream boundary conditions (inlet fixed in both location and angle of entry).
To emphasize this point further, we repeat the same test with the same initial and boundary conditions, but restrict ourselves to a domain with 
$0 \leqslant x \leqslant 10$ as opposed to the domain 
$0 \leqslant x \leqslant 20$ of the previous example. As we will demonstrate in our detailed review of experimental results in our subsequent paper, 
$0 \leqslant x \leqslant 10$ more closely resembles the domain common to actual flume experiments. Figure 13 shows a very dramatic effect with almost all of the growth suppressed in the shorter channel despite the fact that the initial condition remains very close to the wavenumber for maximum growth in spatially periodic linear theory. Returning to our dimensional example, these plots would correspond with flumes of a 
$5\ {\rm m}$ length.
Figure 13. Entrance effects inside shorter region for a system at 
$\alpha =5$ and initial wavenumber of 
$k=1.5$. Dotted line indicates initial condition, whereas thick black line indicates most recent time step. Panel (a) demonstrates periodic upstream boundary conditions, where the inlet location is fixed, but angle of entry is oscillating. Panel (b) demonstrates homogeneous upstream boundary conditions (inlet fixed in both location and angle of entry).
4.3. Linear theory: variations in the parameter $\alpha$
In the IPS theory of meandering rivers, we have found that there is a single dimensionless parameter 
$\alpha$ that accounts for the relative strengths of secondary and irrotational flow. In previous sections, we have given illustrations of various meandering phenomena assuming a value 
$\alpha =5$ which is comparable to the observed values in most real geophysical systems. There is significant variation in the value of 
$\alpha$ in natural systems, however, and it is of interest to investigate the effect of different values of 
$\alpha$. In addition, values in laboratory flumes may show significant departures from this number. We shall see that the value of 
$\alpha$ may have a significant impact on the scaling of experimental equipment – both with respect to length scales and time scales.
To illustrate the impact of natural variation on the spatial growth rate as a function of frequency, we turn to figure 14, which gives the real part of 
$r_1$ at various 
$\alpha$ values. We see that for a given 
$\alpha$, low frequency driving is associated with spatial growth and high frequencies with spatial decay. The maximum real part of 
$(r_1)$ gives the maximum spatial growth that can be achieved. In figure 14, we see that as 
$\alpha$ decreases the frequency for maximum spatial growth also decreases. Thus at lower values of 
$\alpha$, a lower driving frequency is needed to achieve maximum spatial growth. In a laboratory flume, a lower value of 
$\alpha$ would thus require a longer duration experiment to achieve the same number of periods.
Figure 14. Spatial growth (
${\rm Re}(r_1$)) as a function of driving frequency 
$\omega$ for various values of 
$\alpha$.
For the lowest values of 
$\alpha$ in figure 14, it is difficult to see the portions of the curve at low frequency. For clarity, we pick the value 
$\omega _{neutral}$ which gives neutral stability and normalize the abscissa as 
$\omega /\omega _{neutral}$. On the ordinate axis we normalize with the maximum growth rate at that 
$\alpha$. With these changes, we replot in figure 15 to show the changes in the shape of the curves with changing 
$\alpha$. It is interesting to see the shift in the position of the maximum growth with changing 
$\alpha$; this occurs at roughly 
$1/3$ the value of 
$\omega _{neutral}$ for high 
$\alpha$ but drops to less than 
$1/10$ the neutral frequency at the smallest 
$\alpha$.
Figure 15. Normalized spatial growth as a function of normalized driving frequency 
$\omega /\omega _{neutral}$ for various values of 
$\alpha$.
Next, we examine the scaling of the wavenumber as a function of 
$\alpha$ which is given by the imaginary part of 
$r_1$ plotted in figure 16. We observe that the wavenumber increases as 
$\alpha$ increases. The figure also shows the frequencies for neutral stability and maximum spatial growth rate on each curve. Given the clustering of the symbols, it is challenging to determine how the length scale changes with changing 
$\alpha$. For greater clarity, we replot in figure 17 showing the wavelength vs frequency on log–log axes. As before, the circles on each curve show the frequency for maximum growth at that 
$\alpha$. Thus, for each circle plotted, one can read off both the frequency (horizontal axis) and wavelength (vertical axis) for the maximum spatial growth at that 
$\alpha$. In figure 18, we show wave speed as a function of driving frequency. Figures 17 and 18, together, give a convenient record for scaling length and time in the design of experiments with changes in 
$\alpha$. There is a monotonic trend with larger length scales and longer time scales associated with decreasing 
$\alpha$.
Figure 16. Wavenumber (
${\rm Im}(r_1$)) as a function of driving frequency 
$\omega$ at various values of 
$\alpha$. Squares mark conditions for maximum spatial growth. Triangles mark conditions for neutral stability. For each line, values left (right) of the corresponding triangle correspond to regions of spatial growth (decay).
Figure 17. Emerging meander wavelength as a function of driving frequency 
$\omega$ at various values of 
$\alpha$. Circles indicate conditions for maximum spatial growth.
Figure 18. Wave speed as a function of driving frequency 
$\omega$ at various values of 
$\alpha$. Squares mark conditions for maximum spatial growth. Triangles mark conditions for neutral stability.
While figure 17 shows the qualitative trends, it proves convenient to calculate the values of frequency and wavenumber for maximum spatial growth as well as the spatial growth rate achieved. These three quantities are plotted as functions of 
$\alpha$ in figure 19, figure 21 and figure 20, respectively. Observe that the frequency for maximum growth scales roughly as the square of 
$\alpha$.
Figure 19. Driving frequency associated with maximum spatial growth 
$\omega _{max}$ as a function of 
$\alpha$.
Figure 20. Spatial growth 
$\operatorname {Re}(r_1)$ as a function of 
$\alpha$.
Figure 21. Wavenumber 
${\rm Im}(r_1)$ as a function of 
$\alpha$.
It is also instructive look at the real and imaginary parts of 
$r_1$ as a function of 
$\alpha$, and observe that, as 
$\alpha$ drops to low values (figure 20), the spatial growth rate decreases dramatically. Similarly, in figure 21, the wavenumber increases, resulting in larger meanders.
There is one final parameter which is of interest for varying values of 
$\alpha$. For spatial growing oscillatory functions, it is helpful to define an amplification factor, which is the growth in amplitude over one wavelength of oscillation. This is given by the real part of 
$r_1$ times the wavelength or 
${2{\rm \pi} \mathbb {R}(r_{1})}/{\mathbb {I}(r_1)}$. For small amplification factor, there is little change in shape from one wavelength to the next, and a spatially growing wave will exhibit a dynamics similar to its spatially periodic counterpart exhibiting temporal growth. By contrast, for large amplification factor, the dynamics may be significantly different between the two types of instability. The amplification factor is plotted as a function of 
$\alpha$ on semilog axes in figure 22. There are modest changes in amplification factor over normal ranges of 
$\alpha$ but it does begin to drop significantly at very small values of 
$\alpha$.
Figure 22. Amplification factor: ratio of spatial growth to wavenumber in log–log as a function of 
$\alpha$.
5. Comparison with field observations
In previous sections, we have focused on mathematical results and scaling analysis in terms of dimensionless parameters. In this final section, we briefly review some of the field data on natural systems, showing the merit of the simple IPS equations as a fundamental first model for meandering rivers. In addition, we add an original scaling result for the speed of meander propagation.
From section (2.1), we use the characteristic length 
$L^{*}\sim {H_0}/{2C_{f,0}}$ which arises from a balance between the downstream change of momentum and turbulent frictional forces at the river bed. It may be viewed as the distance downstream that a fluid remembers the influence of an upstream perturbation (Edwards & Smith Reference Edwards and Smith2001).
We choose dimensional units (marked by an asterisk) and non-dimensionalize all length scales with respect to 
$L^{*}$, e.g. 
$s^{*} = s\, L^{*}$, etc. To obtain predictions for the most likely meander wavelengths, we need field estimates for the dimensionless parameter 
$\alpha$ or estimates of the Froude number and the scour factor and bed slope. Finally, with known 
$\alpha$, the dimensionless wavenumbers or wavelengths from linear theory may be converted to dimensional units and compared with field observations.
One of the most studied collections of natural meanders lies along the Beatton River in Canada (Hickin & Nanson Reference Hickin and Nanson1975; Nanson & Hickin Reference Nanson and Hickin1983; Parker & Andrews Reference Parker and Andrews1986). Relevant parameter values from Parker & Andrews (Reference Parker and Andrews1986) are given in table 1 yield a value 
$\alpha = 4.59$. In dimensional terms, classic linear stability gives the predicted wavelengths for maximum growth rate and for neutral stability as 1703 and 928 m, respectively. (Parker & Andrews (Reference Parker and Andrews1986) report 890 m as neutral stability, but the estimate was made before a correction was made to the meander-morphodynamics equation by Johannesson & Parker Reference Johannesson and Parker1989b.) Figures from Hickin & Nanson (Reference Hickin and Nanson1975) show a variety of meander shapes from which one may estimate wavelengths of 600, 900, 1500 and 1600 m. These are quite close to the predicted range given the number of assumptions in the model and uncertainty in determining precise wavelengths from irregular meanders.
Table 1. Beatton River parameters (Parker & Andrews Reference Parker and Andrews1986).
Given a reasonable agreement with length scales for meander length, the next question concerns time scales, specifically, what is the lateral migration speed of meanders as a function of geophysical parameters? Here, there is no direct prediction in the literature, but fortunately simple arguments yield a scaling which provides useful information. From section (2.1), an appropriate choice for a reference time scale is
\begin{equation} { t^{*}=\frac{L^{*}}{U^{*}}\frac{1}{E (\tilde{b}/L^{*})} .} \end{equation} The characteristic fluid velocity is defined by the mean river velocity 
$U^{*}={U_0}$. Recall that the characteristic migration velocity of the channel centreline is based on the difference in velocity across the river from one bank to the other. This can be written as the product of the width of the river 
$\tilde {b}$ and a shear rate scaling as 
$U^{*}/L^{*}$. Multiplying by the erosion coefficient yields the dimensional migration velocity
\begin{equation} { u_{migration}^{*} \sim \left(\frac{U^{*}}{L^{*}}\right)\tilde{b}E \sim 2 g^{2/3} I^{2/3} C_f^{1/3} Q^{{-}1/3}\,\tilde{b}^{4/3}E. } \end{equation} This result gives us physical intuition into the rate at which meanders migrate. For example, consistent with intuition, increasing the gradient 
$I$ or the erosion rate, 
$E$, causes an increase in the migration rate. Contrary to intuition, however, it follows from this equation that increasing volumetric flow rate 
$Q$ causes a decrease in migration rate. Another interesting feature is that the width of the river has no impact on the predicted wavelength of meanders (for asymptotically slender rivers) but the lateral migration velocity scales as 
$\tilde {b}^{4/3}$; e.g. an increase of 
$10 \times$ in the river width will lead to meanders moving more than 
$20 \times$ faster.
We find that this scaling result is broadly consistent with experimental data. From Parker & Andrews (Reference Parker and Andrews1986), we find the values presented in table 1 for the Beatton River. This gives us a predicted 
$u_{migration}$ of 
$0.13\ {\rm m}\ {\rm yr}^{-1}$. According to Nanson & Hickin (Reference Nanson and Hickin1983), the true observed rate of migration of the Beatton River is 
$0.16\text {--}0.70\ {\rm m}\ {\rm yr}^{-1}$. Note that 
$E$ is a heuristic value that can be off by a factor of 
$1000$. Therefore, while our scaling equation provides valuable insight, the agreement shown here may be fortuitous.
Funding
The authors would like to thank ExxonMobil Upstream Research Company for support of this research through sponsored research agreement EM07944. T. Sun was helpful in the initial phases of this project. G. Parker provided valuable feedback throughout the work.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Included MATLAB script
We have included a script (weisshigdon2021_JFMscript.m) intended as a guide for experimentalists in the design of flume-based meander experiments. The script accepts as its input a .xlsx file. We have included a sample input file (inputfile.xlsx) with data formatted appropriately for use with the script. The script outputs tables (displayed to the command line and stored as .txt files) and plots. The tables include user inputs, calculated values, dimensionalization parameters, meander predictions for behaviour far downstream from the inlet (as derived from IPS as corrected by Johannesson & Parker Reference Johannesson and Parker1989b) and meander predictions for behaviour close to the inlet (as derived from our own finite-domain theory).
To run the code from a MATLAB prompt, use weisshigdon2021_JFMscript(filename, plotflag) where filename is the name of the input file and plotflag is a Boolean value that enables plotting when set to 1 and suppresses plotting when set to 0. As an example:
> weisshigdon2021_JFMscript(inputfile.xlsx,1)
Note that the input file requires the following experimental design parameters (where 
$\tilde {}$ indicates dimensional units): 
$\tilde {H}$ (channel depth in metres), 
$\tilde {U}$ (average flow velocity in metres/second), 
$I$ (flume slope), 
$\tilde {b}$ (channel half-width in metres), 
$\tilde {L}$ (flume length in metres) and 
$t$ (experimental time in hours). It requires the following sediment parameters: 
$\tilde {D}_{50}$ (average sediment diameter in metres) and 
$\tilde {\rho }_s$ (sediment density in 
${\rm kg}\ {\rm m}^{-3}$). It further requires 
$E$ (erodibility). The input file has two optional inputs, which, if included, will generate additional calculations. If the experimental design includes an upstream perturbation, the perturbation frequency 
$\tilde {\omega }$ (frequency in 1/h) may be included. If that upstream perturbation is applied by means of an inlet that moves via displacement (i.e. perpendicular to the direction of flow) rather than by changing the angle of entry, the displacement amplitude 
$\tilde {B}$ (amplitude in metres) can be entered as well.
The dimensionalization parameters in the output include characteristic length 
$L^{*}$ (2.9a,b) and time 
$t^{*}$ (2.13), as well as non-dimensional length 
$L$ and time 
$t$ (2.12a–d).
Open access
