Inception of Regular Valley Spacing in Fluvial Landscapes: A Linear Stability Analysis

Incipient valley formation in mountainous landscapes is often associated with their presence at a regular spacing under diverse hydroclimatic forcings. Here we provide a formal linear stability theory for a landscape evolution model representing the action of tectonic uplift, diffusive soil creep, and detachment‐limited fluvial erosion. For configurations dominated by only one horizontal length scale, a single dimensionless quantity characterizes the overall system dynamics based on model parameters and boundary conditions. The stability analysis is conducted for smooth and symmetric hillslopes along a long mountain ridge to study the impact of the erosion law form on regular first‐order valley formation. The results provide the critical condition when smooth landscapes become unstable and give rise to a characteristic length scale for incipient valleys, which is related to the scaling exponents that couple fluvial erosion to the specific drainage area and the local slope. The valley spacing at first instability is uniquely related to the ratio of the scaling exponents and widens with an increase in this ratio. We find compelling evidence of sediment transport by diffusive creep and fluvial erosion coupled with the specific drainage area equation as a sufficient mechanism for first‐order valley formation. We finally show that the predictions of the linear stability analysis conform with the results of numerical simulations for different degrees of nonlinearity in the erosion law and agree well with topographic data from a natural landscape.

We focus here on a minimalist landscape evolution model (LEM) representing the action of soil creep, runoff erosion, and tectonic uplift Bonetti et al., 2020). The model contains the least amount of complexity to describe valley initiation and pinpoints the corresponding feedback that induces this phenomenon over long time scales. There are comprehensive LEMs (e.g., Attal et al., 2008;Collins et al., 2004;Coulthard et al., 2013;Van De Wiel et al., 2007) that include spatiotemporal heterogeneity and extensive parametrizations for a wide array of geomorphic processes and serve as a powerful tool to link distinct processes and the resulting morphological changes. These complex numerical models, however, are not suitable for conducting a theoretical analysis of relevant processes and examining their feedback on the emergence of valleys. Not only such an analysis would be based on heavy parametrization and closure assumptions typical of semi-empirical formulations, but it would also need a large amount of knowledge of initial environmental conditions and distinct biophysical properties, making the mathematical framework too complex to support theoretical developments. Using a minimalist LEM allows us to perform a theoretical stability analysis that explains the first-order control of leading geomorphic processes on the emergence of such ubiquitous landscape patterns.

Brief Literature Review
To begin our review of previous investigations of landscape stability, it is useful to orient the reader on the extensive literature on LEMs and how various formulations differ in their descriptions of the coupled water and sediment dynamics. Regarding the water flow modeling in LEMs, the quasi-steady-state approximation is invoked typically by noting that the time scale of the evolving landscape is very long compared to the time scale of the water flow adapting to the modified land surface (Chen et al., 2014). The more comprehensive approaches adopt the full version of the shallow-water equations based on the balance of the inertial, pressure, gravitational, and frictional forces on the water flow (Fowler, 2011;Izumi & Parker, 1995;Smith, 2010). The diffusion-flow approximation dismisses the inertial effects over long time scales (Leopold & Maddock, 1953;Weinmann & Laurenson, 1979). Lastly, the minimalist form of water transport, called the normal-flow approximation, considers a balance between gravitational and frictional forces such that the water flow direction is parallel to the steepest direction of the land surface. Efforts by Gallant and Hutchinson (2011), Bonetti et al. (2018), andPorporato (2022) established that this transport formalism is analogous to the mathematical equation of the specific drainage area for the surface flow fed by a constant unitary rainfall rate (see the derivation in Section 2.1).
Regarding the modeling of long-term fluvial erosion, LEMs typically consider either transport-limited (TL) or detachment-limited (DL) conditions. Some works have also considered intermediate conditions between these two regimes (Davy & Lague, 2009;Pelletier, 2012). Under TL approximation, the fluvial erosion assumes the form of the divergence of the sediment flux which, in turn, is related to the shear stress of the surface flow (Hergarten, 2020;G. E. Tucker & Bras, 1998;Willgoose et al., 1991). In DL approximation, the erosion flux is directly related to the shear stress by flowing water, assuming that the surface resistance to incision is the restricting factor for the erosion rate rather than the hauling capacity of the flow to transport the eroded material (Ahnert, 1987;Howard, 1994). Hence, the mathematical form of DL fluvial erosion becomes a sink term in the LEM, typically expressed as a power-law function of the specific drainage area and local slope.
Within this context, the pioneering work by Smith and Bretherton (1972) presented the first stability analysis of symmetric hillslopes to small lateral perturbations employing a continuous model for water under the normal-flow assumption and TL erosion conditions. This study showed that concave-up portions of a hillslope are unstable to lateral perturbations. Nevertheless, the analysis did not predict a characteristic wavelength for the channel instability but rather an unbound increase in the growth rate for high-frequency perturbations. While this shortcoming has been attributed to the use of normal-flow approximation for the water continuity equation (Izumi & Parker, 1995;Loewenherz-Lawrence, 1994;Smith, 2010), Fowler (2011) suggested that the lack of wavelength selection in the analysis of Smith and Bretherton (1972) could be related to the simplifying assumption of constant coefficients in the governing equations for perturbations. The findings of our work, which keeps the spatial dependency of the coefficients in the governing Equation 13 for perturbations, support the viewpoint of Fowler (2011).
3 of 27 Parker (1995, 2000) conducted a stability analysis for DL fluvial erosion with a threshold value and quasi-steady surface water flow under the effect of inertial, gravitational, pressure, and frictional forces. In particular, Izumi and Parker (1995) described upstream-driven channelization with the channels initiating as the water flux reached a critical value for the DL erosion to occur. They employed an infinitely long flat and tilted hillslope downstream of a drainage divide as the base state for obtaining the solutions. Izumi and Parker (2000) focused on the downstream-driven channelization for a steady concave-down hillslope bounded at the end by a cliff, with the upstream flow state being Froude-subcritical and exactly Froude-critical at the lower boundary. The stability analysis under such special flow conditions and synthetic hillslopes showed the emergence of channel instability for the specified range of model parameters. Finally, Smith (2010) presented a comprehensive mathematical framework for characterizing the channel formation in the TL and DL erosion environments and the quasi-steady surface flow model adaptable to different formulations of normal-flow or diffusion-flow approximations. The initial hillslope was assumed here to be a steady planar profile with a constant slope over which small perturbations could evolve over time.
All the previous theoretical contributions just described considered perturbations on simple but somewhat artificial surfaces. The simple hillslope forms used in these studies facilitated analytical tractability to determine the appearance of well-defined channels, but-not being steady-state solutions of LEMs-had only limited bearing on natural landscape morphologies. In this regard, a more realistic starting point to investigate the conditions of valley formation was pursued by  and Perron et al. (2009). They described the evenly spaced valley formation for numerical solutions of LEM under DL fluvial erosion conditions. These works did not carry out a formal stability analysis, but using numerical simulations they showed that the relative timescale of fluvial erosion to soil creep controls the scale of valley spacing. Employing unchannelized solutions of LEMs with specific boundary conditions as a base solution for a formal linear stability formulation would help to formulate precise criteria for the channelization onset. A preliminary analysis along these lines was conducted by Bonetti et al. (2020) using a minimalist DL-LEM for the special case of unitary exponents of the specific drainage area and topographic slope in the DL erosion law. However, a more complete stability analysis that includes the effect of the nonlinear scaling exponents in the erosion law on the incipient valley formation for base-state solutions of the minimalist LEM is still missing and motivates the work here.

Goal of This Contribution
Within the context outlined above, we focus here on a minimalist LEM in DL erosion conditions and employ the specific drainage area as an approximation for the water discharge, assuming the water flow along the direction of steepest descent of the land surface (i.e., normal-flow hypothesis). We conduct a linear stability analysis of the unchannelized solutions of the governing equations to identify the critical condition under which an initially smooth surface assumes a morphology similar to observed regularly spaced first-order drainage basins. The DL erosion condition is adopted based on the arguments that the bed erosion for the first channelization over the hillslope and low-order valleys is controlled by the erosive power/shear stress of the overland flow rather than the flow capacity to transport the eroded sediments (Howard, 1994;Izumi & Parker, 1995. We consider two symmetric hillslopes along a linear ridgeline as an idealization of a long mountain ridge in a natural landscape to derive unchannelized base-state solutions of the governing equations. Differently from previous contributions, the mathematical forms of the unchannelized solutions are obtained by applying the boundary conditions of water and sediment fluxes in the governing equations and solving for the steady-state (the so-called base-state profile) rather than assuming an arbitrary initial form of the erodible surface.
The stability problem is solved by means of a spectral technique based on the Galerkin projection with numerical quadrature (Canuto et al., 2006), which has been shown to be particularly performant and well suited for morphological problems (Camporeale, 2015;. Employing this strategy, the impact of nonlinearities present in the erosion law on the hillslope stability and the incipient channelization is discussed as erosion gets intensified with respect to soil creep. The predictions of the stability analysis are compared with numerical simulations in a long rectangular domain and also with the topographic data of a natural landscape. The obtained results show that regularly spaced valleys emerge at a certain proportion of fluvial erosion and soil creep. From the water-flow modeling perspective, our results also show that the minimalist normal-flow hypothesis leads to a spatial wavelength preference on the channelization onset under the action of DL erosion and soil creep.

of 27
The article is structured as follows. In Section 2, we present the governing equations of the LEM, along with domain geometry and boundary conditions employed in this study. From the dimensional analysis of the presented LEM, we obtain a single dimensionless quantity that characterizes overall model dynamics. We derive the linearized perturbed equations that are recast in terms of a third-order differential eigenvalue problem. Results from the linear stability analysis are discussed in Section 3. We show the control of the power-law exponents in the DL erosion on the erosion-to-creep threshold for first surface instability and incipient valley spacing. A comparison between stability analysis predictions and results from numerical simulations is carried out. The findings of the stability analysis are also evaluated against topographic data from a natural landscape. In Section 4, we provide the limitations and assumptions in our model and further show how these theoretical results can validate the accuracy of numerical solvers for LEMs. Conclusions and future research directions are finally discussed in Section 5.

Linear Stability Analysis
This section begins with a description of the governing equations of the LEM under the DL fluvial erosion approximation. We then explain the landscape geometry and boundary conditions used to obtain the base-state solutions (smooth landscape) with two opposite and symmetric hillslopes along a long drainage divide in the middle. The stability problem is posed by considering small arbitrary sediment redistribution or weak perturbations over the initially smooth landscapes. All perturbations are assumed to have very small amplitude compared to the initial landscape profiles; hence, the governing equations for perturbations are written in the linearized form. Using dimensional analysis, a non-dimensional parameter,  , is obtained that describes the overall control of erosion, creep, and uplift on the model solutions and perturbations growth. A single differential eigenvalue problem is derived for perturbation dynamics, where the non-dimensional form depends only on  for given erosion law exponents.

Governing Equations
The coupled dynamics of the landscape elevation and surface water fields can be written in general as Equation 1 expresses the temporal evolution of the elevation field z under the action of tectonic uplift U, sediment flux due to soil creep f c , and the flux transported by fluvial erosion, f e . Soil creep is used to represent a combined effect of various biophysical processes that result in the slow movement of soil over the hillslope. Surface and subsurface processes inducing this motion include animal burrowing, falling trees, wetting/drying of the upper soil layer, and freezing/thawing cycles of the pore water in the subsurface (Carson & Kirkby, 1972;Gabet et al., 2003). The combined effect of these movements smooths the topography so that the flux is written in an average sense as a diffusion term, f c = −D c ∇z, where D c is a coefficient based on the combined efficiency of different soil creep processes (Culling, 1963). In the DL condition, the amount of erosion flux is assumed proportional to the shear stress by the flowing runoff over the surface such that the erosion term is approximated as a sink term ∇ ⋅ ≈ ′ |∇ | , where ′ is an erosion coefficient, q is the specific runoff or the surface flow rate, and m and n are positive exponents (Howard, 1994;Whipple & Tucker, 1999).
In Equation 2, R represents a runoff-producing rainfall rate, that is, the amount of precipitated water contributing to runoff production q in the direction of n. Under the observation that the changes in surface water height h occur immediately compared to the time scale of the evolving landscape, a quasi-steady-state approximation (∂h/∂t = 0) is considered here. Assuming the water flow in the direction of steepest descent of the surface (n = −∇z/|∇z|) and a time-averaged runoff-producing rainfall rate R 0 , the water discharge q is proportional to the specific drainage area (= ∕ 0) . As a result, Equation 2 becomes the governing equation for the specific drainage area (Bonetti et al., 2020;Porporato, 2022). Employing this proportionality between q and a, the sink term for DL erosion can be written as K e a m |∇z| n , where = ′ 0 . The above conditions lead to the final form of the governing equations as 5 of 27 Through the coupling between a and z, this minimalist LEM captures the essential long-term dynamics of landscape evolution. Fluvial erosion acts as a sink term, while soil creep acts as a diffusion term in Equation 3. Erosion excavates sediments at locations where the accumulation of the specific drainage area is high. This yields a higher surface gradient at those locations with a further increase in a, thus enforcing the increased erosion and flow accumulation again. This feedback loop between the emerging topography and the accumulated specific drainage area tends to carve a preferential path over time, if the surface smoothing effect by the diffusive creep is not sufficient, with the progression toward the dissected landscape.

Boundary Conditions and 1D Smooth Morphology
The smooth landscape geometry considered here is the steady state solution of the previous model in conditions of no channels (i.e., no y dependence), consisting of a long mountain ridge of width l x with the drainage divide in the middle and two symmetric hillslopes on either side, as shown in Figure 1a. The x-axis points along the hillslope and the y-axis denotes the direction of the ridgeline/drainage divide. The surface has the highest elevation at the divide with a monotonic drop along either hillslope depending only on x. We consider the hillslope bases at the same fixed elevation level (taken as zero for reference) and no water or sediment flux movement across the drainage divide in the middle. These boundary conditions are consistent with those adopted in earlier studies on the analysis of 1D hillslope morphology (Bonetti et al., 2019;Loewenherz, 1991;Smith & Bretherton, 1972).

of 27
Such unchannelized solutions at steady-state can be obtained by using dz/dx = a = 0 at x = 0 and z = 0 at x = l x /2 in Equations 3 and 4 for the 1D transect. With elevation declining monotonically on either side of the divide, the steady-state solution for specific drainage area is simply the relation a 0 = x with x ∈ [0, l x /2]. Namely, it increases linearly with the distance from the ridgeline, as shown in Figure 1b. The subscript 0 is used here to represent the base state. The steady-state solution for the smooth elevation field z(x) depends on the value of exponents m and n. This solution can be obtained analytically only for m = n = 1, where it takes the form of a generalized hypergeometric function Bonetti et al., 2019Bonetti et al., , 2020, while it has to be obtained numerically for generic exponents m and n, as discussed in Section 3. We refer to this solution as z 0 (x) in the stability analysis formulation.

Linearized Perturbed Governing Equations
Having established the base-state solutions, we can analyze the conditions of DL erosion, creep, and uplift for which an unchannelized landscape like the one portrayed in Figure 1 loses its stability to small perturbations and give rise to the patterns of periodically spaced valleys. For a smooth landscape that has attained a steady state, a very small and arbitrary spatial redistribution of sediments would slightly modify the smooth topography. If the landscape evolution under the governing geomorphological processes dampens any arbitrary disturbance over time, the topography reverts to its earlier state. Physically, one would speculate that, as the extent of the erosion increases relative to the smoothening effects of soil creep, there should be a point when the smooth surface becomes unstable to these infinitesimal perturbations. The perturbation that grows fastest at the instability onset would alter the state of the unchannelized landscape to the one with valleys spaced at a typical spatial scale corresponding to the wavelength of the fastest-growing perturbation.
A normal-mode linear stability analysis provides a formal way to systematically detect the inception of the surface instability on the base-state solutions for every possible spatial perturbation that leads to the formation of first-order valleys at a specific length scale. In this methodology, such infinitesimal arbitrary perturbations are resolved into normal modes (e.g., sines and cosines). Since each spatial mode satisfies the linearized system individually, their superposition can describe the growth or decay of weak perturbations over the base-state solutions (Drazin & Reid, 2004). We refer to the following references for an extensive description of this approach and its applications in various other physical systems (Chandrasekhar, 2013;Cross & Hohenberg, 1993;Koch & Meinhardt, 1994;Vlase et al., 2019).
With infinitesimal arbitrary perturbations in the base-state solutions, the modified elevation and specific drainage area fields can be written as ( ) = 0( ) +̃( ) and ( ) = 0( ) +̃( ) . Here z 0 (x) and a 0 (x) are the unchannelized 1D solutions at steady-state discussed in Section 2.2. ( ) and ( ) denote the weak perturbations over the unchannelized solutions. We assume here homogeneous boundary conditions for the weak perturbations, namely = 0 at x = l x /2 and ∕ = = 0 at x = 0. The mathematical expressions for and are written as̃= where ψ(x) and ϕ(x) represent perturbation amplitudes varying along the hillslope with angular wavenumber k in the y-direction and initial growth rate σ (c.c. refers to complex conjugation). Depending on σ being higher or lower than zero, the perturbation of a particular wavenumber k respectively grows or decays over time.
An illustration of the modified elevation field with the weak perturbation form taken in Equations 5 and 6 is shown in Figure 2. The perturbed state and the initial smooth topography can be easily distinguished based on the horizontally projected streamlines and contour lines in Figures 1b and 2b. For the smooth landscape, contour lines run parallel to the hillslope base, and the streamlines initiate from the divide in the middle and run side by side till the hillslope base ( Figure 1b). In the perturbed state shown in Figure 2b, contour lines have the periodicity set by the spatial wavelength of the perturbation, which results in the periodic convergence and divergence of streamlines starting from the divide in the middle. The fastest growing perturbation of given wavenumber k reaching the critical balance of erosion to creep for channelization would result in the emergence of equally spaced valleys at a distance 2π/k in the landscape.

of 27
As shown in Appendix A, substituting the forms of perturbations from Equations 5 and 6 for the modified z and a fields in governing Equations 3 and 4 and linearizing, the coupled governing equations for perturbations ψ(x) and ϕ(x) become is the steady-state unchannelized topographic slope. The solutions to these coupled equations for perturbations provide the growth rates σ for weak perturbations of different values of wavenumber k under the distinct balance of the governing processes.

Dimensional Analysis
To quantify the overall control of creep, erosion, and uplift on the model solutions and the corresponding growth of the weak perturbations (Equations 7 and 8), we non-dimensionalize the system before proceeding toward the solution. There are three primary dimensions in the present problem, namely horizontal directions, vertical direction, and time. The width of the ridge l x is taken as the horizontal length scale  . The vertical scale  is taken as 2 and time scale  is defined as 2 . Using these scales  ,  , and  , the following non-dimensional quantities are obtained: Based on these quantities, the non-dimensionalized forms of governing Equations 3 and 4 read Figure 2. Schematic diagram presenting the perturbed state of the landscape used in the normal-mode analysis and the homogeneous boundary conditions. The weak perturbation has been exaggerated for better visualization. (a) A representative 3D surface z is displayed, where x-axis/y-axis denotes the direction along the hillslope/ridgeline. The perturbation with wavenumber k corresponds to the spatial wavelength (= 2 ∕ ) . (b) The horizontal projection of the surface is shown with streamlines (in blue) perpendicular to the projected contour lines (in brown). The projected streamlines converge at the emerging equally spaced valleys and diverge at the corresponding interlocked ridges.

of 27
This analysis demonstrates that a single non-dimensional quantity  controls the overall system dynamics based on global parameters (K e , D c , and U) and boundary conditions (width of the ridge, l x ) for given exponents m and n in the erosion term.
 provides a measure of dynamic similitude, which means that if two distinct model configurations have the same  value for prescribed m and n, they are dynamically equivalent . Physically speaking,  represents the relative proportion of DL fluvial erosion and diffusive soil creep to transport the sediment inflow in the landscape by uplift rate. This index is similar to the global Reynolds number Re that expresses the balance of inertial to viscous forces in fluid flows with one characteristic length scale (Kundu et al., 2015;Madlener et al., 2009). An increase in  describes the transitions in the landscape topography from an unchannelized surface to first-order non-branched valleys to complex patterns of ridges and valleys with a cascade of branching Bonetti et al., 2020). The growth of surface instabilities with finer and finer branching in landscapes with rising  values is reminiscent of an increase in flow instabilities and the shift from laminar to turbulent flow conditions as Re increases Porporato, 2022).
From Equation 10, we can infer the overall impact of changes in model parameters and boundary conditions on the obtained topographies. For example, either an increased magnitude of K e (diminished surface resistance to fluvial erosion and/or improved runoff-producing rainfall rate) or a decreased value of D c would raise the  value, showing an enhanced tendency of the LEM topography to branch due to relatively higher erosion component. Similarly, an increase in l x with unchanged model parameters, that is, a larger topography with identical geomorphic process rates has a higher tendency to become unstable. Interestingly, the differences in the uplift rate do not influence the balance of erosion and soil creep for n = 1 (i.e., invariant  value); however, an increased uplift augments the erosion component for n > 1, while it enhances the diffusive creep component of the LEM for n < 1.
The shift in the impact of uplift rate U on  based on the value of n compared to unity has a physical underpinning. The uplift rate is the only source term for surface elevation in the LEM and primarily sets the vertical scale  of the LEM solutions. Diffusive creep and DL fluvial erosion merely transport the sediment across and outside the landscape. For example, with the boundary conditions considered here, the steady-state landscape would be a flat surface if the uplift rate were zero, as there is no sediment influx to counterbalance the losses. Analyzing the scaling of diffusive creep and erosion sink terms with the vertical scale, one notices that the creep term scales linearly with vertical scale  and the erosion sink term scales as  (see Equation B1 in Appendix B). For an increase in U, the response of the diffusion term is always linear, while the erosion term responds either sub-linearly for n < 1 or super-linearly for n > 1. Hence, the global balance of creep and erosion intensity with a boosted uplift rate shifts toward either process depending on the value of exponent n being smaller or greater than unity, as encapsulated in the definition of  in Equation 10. A more detailed discussion of changes in  , m, and n on the LEM solutions is presented in Appendix B, with corroboration of numerical simulation results.
For the above-defined scales  ,  , and  , we introduce the following dimensionless quantities: =  , =  , ̂=  , ̂=  , and ̂0 = 0   . Using these quantities, the governing Equations 7 and 8 for weak perturbations read̂= 9 of 27 The coefficients in the above equations contain spatial dependencies ( and ̂0 ), which restrict the applicability of analytical treatment to estimate growth rates for weak perturbations. Here we employ the spectral technique to compute the solutions to this stability problem, as described in the following section.

The Eigenvalue Problem
To obtain the solutions for the evolution of the perturbations using the spectral technique, we recast the reference system from to (= 4̂− 1) and keep the solution domain between −1 and 1, so that the Legendre polynomials could be used as the basis functions in the spectral solver. By applying this reference change and substituting the value of from Equation 12 in Equation 11, a single differential equation in terms of ̂ can be obtained as The homogeneous boundary conditions for perturbations can be re-written as ̂ (̂ = −1) =̂ ′′ (̂ = −1) =̂ ′ (̂ = 1) = 0 . We refer to Appendix C for the derivation of the above equation as well as boundary conditions in terms of ̂(̂) . The expressions for all coefficients are provided in Table C1 (Appendix C).
Equation 13 with the imposed boundary conditions forms an eigenvalue problem, where non-zero solutions ̂(̂) exist for unique (eigen)values of the growth rate . This system can be solved to compute the growth rate for perturbation ̂(̂) of wavenumber ̂ with increasing  values, or say, increased erosion component of the LEM compared to diffusive creep. By gradually increasing  , a critical value of this dimensionless global quantity,  , can be found for which at least one of many possible perturbations starts growing with a positive rate .  marks the crucial balance of erosion to creep beyond which the smooth landscape is unstable to small perturbations. By tracking the wavenumber k cr with the highest positive growth rate at  , the spacing between emerged first-order valleys λ cr = 2π/k cr can be computed. Hence, the required proportion of fluvial erosion and diffusive creep ( ) , as well as the incipient valley spacing (λ cr ), can be obtained by replicating this approach for different degrees of nonlinearities in exponents m and n.
To proceed toward a solution, we converted the differential problem of Equation 13 into an integral form. This is usually referred to as a weak formulation of the problem due to a reduction in the differentiability constraint of the solution. The weak formulation was then solved by utilizing a spectral technique based on the Galerkin projection with numerical quadrature (Canuto et al., 2006). We employed the algorithm proposed by Swarztrauber (2003) to compute quadrature points and weights for the numerical quadrature. To guarantee an acceptable spectral accuracy, we used 200 points between −1 and 1 for the present results. A detailed explanation of the developed methodology and the spectral solver is provided in Appendix D.

The Emergence of First-Order Valleys
The unchannelized slope ̂ 0(̂ ) as well as its first and second derivatives for different values of  are needed to finalize the form of non-constant coefficients and solve the eigenvalue problem posed in Equation 13; see also Table C1. These expressions are analytically attainable only for the case of unitary exponents m and n, where the unchannelized slope and its derivatives take the form of Dawson functions (see Equation 15 below). For any other values of m and n, these derivatives have to be obtained numerically. ̂0 (= | 0∕ |) can be computed by first recasting the 1D form of Equation 9 in terms of ̂0 at steady-state aŝ which then can be integrated numerically with appropriate boundary conditions for any m and n values.
We solved here the differential Equation 14 with initial value ̂ 0(̂ = −1) = 0 . Once the numerical solution of ̂0 was obtained, the form of ̂0 ′ was computed by using Equation 14 at the discrete quadrature points. ̂0 ′′ was 10 of 27 then obtained by using the second-order accurate central difference of ̂0 ′ at the interior quadrature points and first-order accurate finite difference at the boundary points.

Verification of the Spectral Solver
We first performed a code verification to ensure that the spectral solver (using the numerical form of slope and its derivatives) correctly solves the stability problem without any programming error during the solver development (Oberkampf & Roy, 2010;Roache, 1998). For that, the results using numerical integration of ̂0 (and its derivatives) were compared with the linear stability results employing the analytical solution for ̂0 for the case of m = n = 1, which iŝ0 where  (⋅) is the Dawson function Bonetti et al., 2019Bonetti et al., , 2020.
Results from the linear stability analysis are shown for ridge width l x = 100 m in Figure 3 using numerical integration of Equation 14 for the base-state solutions. In panel a, each curve corresponds to a particular value of  and depicts the growth rate σ for lateral perturbations of different wavenumbers. For low values of  , all spatial perturbations have a negative growth rate, which indicates that the diffusive creep is strong enough compared to the fluvial erosion component in these conditions to fill the emerging surface instability of varying wavenumbers.
As the relative importance of erosion grows in the model solutions, a critical proportion of erosion to creep is observed, for which the fastest growth rate becomes slightly greater than zero for an intermediate wavenumber.
The red curve indicates this critical balance,  ≈ 37 , where the maximum of the curve designates the characteristic valley spacing λ cr ≈ 41 m. Any weak perturbation of lower or higher wavenumber decays at this juncture, which signifies an underlying tendency in the minimalist LEM for scale selection of the emerging valleys. The LEM presents here the type I linear instability characterized in Cross and Hohenberg (1993), similar to the Orr-Sommerfeld stability problem for the plane Poiseuille flow. The instability prediction using numerical solutions

For
 >  , there is a band of wavenumbers with a positive growth rate, showing the unstable state of the smooth landscape. In this regime, the peak of growth rate curves shifts toward the right, revealing the emergence of narrower valleys with an increase in  beyond  . This result agrees well with the observation of reduced spacing between first-order valleys with a relative increase in DL erosion compared to diffusive creep in the numerical simulations performed by  in a rectangular domain.
The marginal/neutral stability (black) curve is shown in Figure 3b to demonstrate the transition of an unchannelized hillslope from a stable to unstable state as the critical value of  is reached and the unstable wavenumbers (the red region) appear beyond  . The vertical red dotted line marks the fastest-growing wavenumber k cr at the channelization threshold  , which is indicated by the horizontal dotted line. The dashed red line in the unstable regime shows a reduction in the valley spacing with the growing value of  . After verifying the spectral solver for the special case of exponents m = n = 1, we studied the effect of non-unitary exponent values on the characteristic valley spacing λ cr at the first instability  , as discussed in the following sections. Comparison of the linear stability predictions with the numerical simulations of the complete LEM is done for the ridge width l x = 100 m.

The Influence of Different m Values
Values of exponents m and n describe the coupling of the specific drainage area and local slope in the fluvial erosion mechanism of the LEM. In modeling studies of landscape evolution, the value of n is typically assumed to be unity with a range taken between 0.67 and 1.67 (Perron et al., 2009;Seidl et al., 1992). The value of ratio m/n is reported between 0.35 and 0.8 in the context of equilibrium stream profiles under DL erosion and uniform uplift from digital elevation models, field and map studies (Bonetti et al., 2019;Slingerland et al., 1998;Snyder et al., 2000;Tarboton et al., 1991;G. Tucker & Whipple, 2002). The ratio m/n equal to 0.5 is generally used as the base case of the Optimal Channel Network (OCN) theory due to its close resemblance with scaling laws obtained in fluvial landscapes with negligible diffusive soil creep, that is,  → ∞ (Banavar et al., 1997;Hooshyar et al., 2020;Rinaldo et al., 2014;Rodriguez-Iturbe & Rinaldo, 2001). Whipple and Tucker (1999) and Lague (2014) provide a comprehensive review of this power-law relationship and the typical values of m and n derived from either shear stress or unit stream power law along bedrock channels, along with the evidence from field studies from tectonically active regions.
We discuss here the role of the exponent m on the emergence of first-order valleys for n = 1, while the non-unity values of n are examined further in Section 3.3. Figure 4a displays the marginal stability curves obtained for eight values of m between 0.125 (red) and 1 (blue), where the corresponding horizontal lines represent the channelization threshold  and the vertical lines mark the fastest-growing wavenumber k cr at this critical juncture. Figures 4b and 4c displays the dependency of the channelization threshold and emerging valley spacing on the value of m. Specifically, as m decreases, an increase in the critical  value,  , is observed together with the formation of narrower incipient valleys.
From the definition of  provided in Equation 10, it can be seen that the overall system's behavior is governed primarily by the ratio of the erosion coefficient (K e ) to the soil creep coefficient (D c ) for a mountain ridge of fixed width l x and exponent n = 1. The inset of panel b displays the above result in terms of the critical ratio of K e to D c , ∕ ( =  ∕ +1 ) , required to initiate the valley formation in natural landscapes under unchanged boundary conditions. ∕ grows by a factor of almost 100 as m reduces by a factor of 8 from 1.0 to 0.125. This required increase in the critical erosion-to-creep intensity for initiating valley formation as m approaches zero reveals the significance of non-local feedback conveyed by specific drainage area in the erosion mechanism and the ensuing self-organization of the fluvial landscape.

Numerical Simulations for Generic m and n
We compared the predictions of the linear stability analysis with the instance of the first channelization by using the numerical algorithm introduced in Anand et al. (2020) Figure 5a. This behavior is an extension of the trend observed for the case of varying m for n = 1 in Section 3.2. On the contrary,  rises with an increase in n at a fixed value of m. A good agreement between the stability predictions and the numerical simulations of the model is observed for different values of m and n explored here. These results can be further employed to express the required changes in the model parameters to achieve the critical balance of erosion to diffusive creep based on the property of dynamic similitude at the same  . Fixing the soil creep coefficient D e , uplift rate U, and exponent n for the landscape of given ridge width l x , Equation 10 reveals that a higher value of the erosion coefficient K e with decreasing values of exponent m is always needed in this scenario as  scales as × .
The effect of varying exponents in the DL erosion law on the scale of the first dissection over a smooth landscape is shown in Figure 5b. A decreasing value of m for a fixed n value shows the formation of narrower valleys at the first instability. The exponent n has little bearing on the preferential channelization scale at a given m. Interestingly, the variation of λ cr as a function of m for distinct n values appears to be a shifted version of the same curve. Re-scaling these results based on the ratio m/n reveals a remarkable collapse of the plots in Figure 5b as a single logarithmic function of the ratio of two exponents m/n for the explored parameter range in Figure 5c. This result provides quantitative evidence of the first ridge/valley wavelength selection based on the ratio of exponents coupling DL fluvial erosion with specific drainage area (m) and local slope (n), corroborated by the results of numerical simulations.  Predictions from the stability analysis agree fairly well with the occurrence of first channelization and incipient valley spacing obtained in the steady-state solutions from numerical modeling ( Figure 5). The slight difference in  and λ cr values for the numerical model and the stability analysis hints at the nonlinear interactions (higher-order terms in the governing equations of the perturbations) disregarded in the linear stability formulation that, despite being small, are present in the numerical simulations of the nonlinear LEM. Furthermore, the finite resolution of the discrete domain and the underlying grid orientation affect the simulation results of pattern-forming systems like the minimalist LEM. These factors result in a small distortion of the solution geometry mapped onto the discrete domain geometry. Unit grid spacing used here for the simulations on the rectangular grid could also lead to dissimilarities in the estimation of valley spacing compared to the linear stability predictions.
The presented analysis on the effect of the DL erosion exponents on the incipient valley formation in a long rectangular domain agrees well with the observations regarding the emerging first-order valleys in a square domain using numerical simulations in Bonetti et al. (2020). For example, the simulated landscapes in Figure 3 of Bonetti et al. (2020) show channelization and subsequent branching at higher  values as n increases (0.7, 1.0, and 1.3) at fixed m = 0.5. Similarly, the appearance of narrower primary valleys has been noted in the study for lower m values at a given exponent n.

Comparison With Regular Valley Spacing in a Natural Landscape
We compared predictions from the linear stability analysis with observations of first-order valley formation in a natural landscape dominated by diffusive creep and fluvial erosion. The landscape examined here is a portion of Gabilan Mesa in California characterized by a Mediterranean climate and oak-savanna (lightly forested grassland) vegetation cover and was previously investigated in ; Perron, Dietrich, and Kirchner (2008)  14 of 27 The displayed terrain has NE-SW trending channels (green) and evenly spaced intervening hillslopes along the ridges (brown) as shown in Figures 7a and 7b, with the distance between the two principal channels setting the width of the long ridge to be roughly 550 m and the valley spacing along the ridges around 163 ± 11 m noted in ; Perron et al. (2009). The spacing values were measured using 2D Fourier spectra for the landscape elevation after removing the spatial trend, which assumes the mean and variance of the provided signal to be roughly constant (Beamish & Priestley, 1981;. The external geometrical constraints for a long ridge surrounded by two main channels at roughly 550 m spacing resembles the domain geometry of a long ridge with a constant width l x employed in the stability formulation. In addition, the presence of two parallel flowing channels at the base of the ridge also agrees well with fixed elevation boundary conditions at the hillslope base used in the stability analysis here. Assuming exponent n = 1, the values of D c /K e = 124 ± 3 and m = 0.35 ± 0.003 were computed by using the shapes of hilltops and stream profiles for the given topography in Perron et al. (2009). Employing these values of the parameters and ridge width l x = 550 m with relative uncertainties l x and n assumed to be 2.5%, we estimated the value of  to be 40.4 ± 7.3. We conducted the linear stability analysis for m = 0.35 ± 0.003 and n = 1.0 ± 0.025 and tracked the first channelization instance along with the dominant valley spacing in the calculated  range. The stability analysis results predict the value of  ≈ 44 ± 3.5 , which falls in the estimated  range for the landscape. The dominant valley spacing is computed to be 175 +6 −38 m, which is also in line with the measured spacing around 163 ± 11 m in Perron et al. (2009). Figure 7c shows the stability analysis result for average values of the parameters m = 0.35, n = 1, and l x = 550 m.
The satisfactory agreement between the first-order valley spacing obtained from the stability analysis with those acquired by the topographic measurements of the landscape suggests that the linear stability formulation of the minimalist LEM captures well the feedback between the competing diffusive creep and fluvial erosion for the

Discussion
Previous theoretical studies on landscape stability analysis and the onset of first-order drainage basins employed simplified forms of base-state surfaces in the perturbation equations. The shapes of such surfaces were selected to maintain analytical tractability but had no direct linkage with the actual LEM equations. For this reason, not only a proper validation of the LEM solutions was unfeasible, but it also remained unclear whether the regular valley formation observed in the numerical simulations of the LEM was due to the artificial diffusion induced by grid discretization. As we now discuss in detail, results from the present work fill this gap and serve as a validation benchmark for numerical LEM solvers. Before doing that, it is important to assess the assumptions and limitations of our modeling framework.

Limitations of the Present Analysis
In the minimalist LEM used here the normal-flow approximation, utilizing the specific drainage area field a, assumes that flow lines coincide with the streamlines of the topographic field. If inertial and diffusion effects also influence the water discharge, streamlines of the underlying topography and water flow lines do not coincide. Hence, a clear distinction between the two is required in those conditions.
The specific drainage area a (and total drainage area A employed in spatially discrete LEM formulations (Howard, 1994;G. E. Tucker & Hancock, 2010)) is a morphometric variable that serves as a proxy for surface water flow and performs well to represent the effects of erosion on landscape evolution over geological time scales. Theoretically, a is described at a point while A is defined for a finite width of the contour line Bonetti et al., 2018;Gallant & Hutchinson, 2011). Evidently, any LEM based only on the drainage area field cannot differentiate between channelized and overland flow as there is no underlying physical mechanism for 16 of 27 river dynamics. Ridges and valleys are also defined based on the convergence and divergence of streamlines in the landscape geometry, with valleys characterizing the domain locations of negative plan (or contour) curvature and ridges as the areas of flow divergence, represented by positive plan curvature (Figure 2). It is crucial to bear in mind these points when analyzing the results in comparison with natural landscapes as well as for application purposes.
The present analysis is limited to long rectangular domains with fixed elevation boundary conditions. This was chosen to have base-state solutions that are relatively simple while being an exact solution for the specific drainage area a and the surface elevation z. The idealization of the domain geometry enabled us to disentangle the complex interactions of the internal physics of the model from the external constraints posed by more realistic boundary conditions. Extensions to other regular domains are possible, such as a smooth topography on a circular geometry, for which the analytical expression of a is r/2 (r being the radius of the domain). The imposed boundary conditions and the trial functions in the spectral solver need to be modified accordingly to approximate the linearized perturbations and predict the emerging modes of periodic valleys in respective regular domains.
The linear stability theory shows here the first-order periodic valley emergence under the varying proportion of creep, DL erosion, and uplift under small sediment redistribution in a smooth landscape. Although comparison with data from a natural landscape is very promising, we emphasize that these predictions are theoretically accurate only before the nonlinear interactions become prominent with further expansion in the erosion component, that is, the finite size of the emerged valleys becomes dominant in the overall dynamics of pattern formation. Any extension of these results beyond this regime is purely based on empirical evidence. This is an important point to consider when the present analysis is utilized for further applications in natural landscapes, which have estimated  value far beyond the critical threshold  . To this aim, an extension of the presented formulation to a fully nonlinear stability analysis could greatly improve our understanding of these finite-sized nonlinear interactions and their role in the perturbations growth and the preference for dominant valley spacing.

Testing LEM Numerical Solvers
Modeling studies simulating landscape evolution under the varying balance of DL fluvial erosion and diffusive creep showed the presence of first-order valleys at roughly uniform spacing Perron et al., 2009). However, these studies were based on numerical simulations, which are necessarily approximated. Our work based on linear stability theory provides compelling evidence of first-order valley emergence with a good agreement between simulation results and stability predictions at small  values. The stability analysis demonstrates that the shift from a smooth topography to a channelized one occurs in the LEM at a critical proportion of DL erosion to creep transport and that, at this juncture, the solutions exhibit a preference for the spatial scale of the first dissection in the linear regime. Such a length scale strongly depends on the ratio of exponents m and n in the DL erosion law.
Results from the linear stability analysis can serve as benchmarks for testing the accuracy of the numerical solvers for landscape evolution simulations. Before the first instability occurs, the mean elevation profile of the steady-state simulation for a long rectangular domain should conform with the unchannelized 1D elevation profile. In particular, Equation 14 here provides the variation of the slope along the smooth 1D transect as a function of  , a result that can be readily obtained for the given value of the exponents m and n. The accuracy of a numerical solver can be assessed by checking the  value at which the simulation solution starts deviating from the 1D solution due to channelization and comparing it with the  value obtained from the stability analysis. As an example, we compared the theoretical expression of unchannelized 1D slope with the slope of the mean elevation profile obtained from the numerical simulations for different  values. To this purpose, we numerically integrated Equation 14 for m = n = 1 and calculated the slope at the hillslope base (̂0 (̂= 1) ) . Following Anand et al. (2020), a was numerically approximated as A/Δx, where Δx is the grid resolution. We used D∞ and D8 flow-direction methods for computing A (O'Callaghan & Mark, 1984;Tarboton, 1997). The LEM simulations were performed on a long rectangular domain (width = 100 m and length = 500 m, same as in Section 3.3) using these two flow-direction methods and fixed elevation boundary conditions. We compared the boundary slope of the mean elevation profile along the length (ignoring the last 100 m of the domain) from the simulations with the solution given by Equation 14 for increasing values of  . Figure

10.1029/2022JF006716
17 of 27 symbols), the mean boundary slope starts deviating around  ≈ 32 with the emergence of first channel instability, which is close to the linear stability prediction of  ≈ 37 (Section 3.1). The boundary slope starts deviating at  ≈ 90 for the LEM solver using the D8 method (filled violet symbols), which indicates that the transition from smooth to dissected landscape is not captured well by the simulations based on the D8 flow-direction method.

Conclusions
The linear stability analysis of the governing equations of a minimalist DL-LEM allowed us to quantify the role of DL erosion and diffusive creep on the evenly spaced valley formation. The use of the spectral method made it feasible to compute solutions to the posed stability problem in the form of a differential equation where non-constant coefficients elude analytical tractability Canuto et al., 2006). The flexibility provided by the spectral method can be extended further to quantify the effects of erosion thresholds, spatially varying parameterization, etc., on the first channelization.
Results showed that first-order valleys with spacing λ cr emerge at a specific proportion of fluvial erosion and soil creep given by the critical value of the non-dimensional index  . We obtained the dependency of λ cr and  on the exponents m and n in the erosion mechanism. In particular, a reduction in m for a fixed value of n increases the  threshold, indicating that a higher proportion of the fluvial erosion to the diffusive soil creep in the LEM is required to carve the hillslope for channelization as the relative importance of the specific drainage area in the erosion mechanism diminishes. Conversely, the  threshold for the incipient ridge/valley formation rises with an increase in the value of n for a particular m value. Only the ratio of exponents m/n influences the selection of a characteristic valley spacing λ cr with narrower valleys emerging for the declining values of m/n. A close agreement between the linear stability analysis and numerical simulations of the LEM in a long rectangular domain was observed for the inception of the regularly spaced valleys at a certain erosion threshold. Predictions from the stability analysis were further validated here by using topographic data from a natural landscape.
One of the main results of our linear stability theory is that the counteracting effects of diffusive soil creep and DL fluvial erosion coupled with the specific drainage area equation give rise to a minimalist model capable of regular first-order valley formation. Obtained solutions to the posed stability problem demonstrate that preserving the spatial dependency of the base-state solutions and the coefficients of the final eigenvalue problem allows a characteristic wavenumber selection based on the erosion law form, providing an explicit connection between the nonlinear erosion feedbacks and the spectral signature of channelization in natural landscapes.

of 27
Extension of the current formulation of the linear stability analysis to more complex transport models is feasible. A subject of particular interest for future work would be the first-order valley formation under TL erosion, where the sink term in Equation 3 is replaced by the divergence term for the sediment transport based on a similar power-law relationship to the specific drainage area and local slope. The present work on valley formation by fluvial erosion can also be applied to the models describing glacial landscape evolution (Deal & Prasicek, 2021;Prasicek et al., 2015, and references therein). A strong parallel between the DL fluvial erosion model and the minimalist glacial erosion model (see the similar forms of Equations 3 and 4 in the present work and the governing Equation 22 in Deal and Prasicek (2021)) offers the potential use of the developed theory for understanding the landscape evolution governed by glacial and a mixed glacial-fluvial sediment transport environment.

Appendix A: Linearized Perturbed Equations and Boundary Conditions
The modified elevation and specific drainage area fields with weak perturbations can be written as where z 0 (x) and a 0 (x) are the steady-state unchannelized solutions; and denote perturbation fields.
Using Equation A1, the updated topographic gradient vector becomes is the unchannelized local slope, i is the unit vector in x−axis direction, and j is the unit vector in the direction of y−axis. Employing this form of the gradient, the linearized expression for the updated topographic slope is written as Employing Equations A1, A2, and A4, the governing Equation 3 for the elevation field z(x, y, t) can be updated to = ∇ 2 0 + ∇ 2 − ( 0 + ) where ( 0 +̃) = 0 + −1 0̃ for small perturbation . Writing (A6) Using Equations A3 and A4, the unit vector in the direction of steepest descent of the updated elevation field is We employ the mathematical expressions for ( ) and ( ) as = ( ) exp ( + ) + c.c., where ψ(x) and ϕ(x) denote perturbation amplitudes varying along the hillslope with angular wavenumber k in the y-direction and initial growth rate σ (c.c. refers to the complex conjugation). Substituting these in Equations A6 and A8, we write the coupled equations for ψ(x) and ϕ(x) as gives Under the assumption that S 0 (x) behaves linearly in the limit x → 0, we get the boundary condition 2 2 ( = 0) = 0 .

Appendix B: Non-Dimensionalization of the LEM
The governing Equations 3 and 4 of the model have three primary dimensions, namely horizontal directions, vertical direction, and time. The treatment of vertical and horizontal directions as independent dimensions is justified based on the fact that the elevation field is obtained from the sediment budget and it represents volume per unit of ground area. It can also be written as the mass per unit area by supposing a constant sediment density Porporato, 2022). The elevation field in the vertical direction, hence, is physically distinguishable from the horizontal directions (corresponding to the so-called directional dimensional analysis) that measure the spatial extent of the domain and establish the boundary conditions for the model. Since we consider here an infinitely long strip of finite width (l y ≫ l x ), one horizontal length scale enters the dimensional analysis of the physical law. Representing the horizontal length scale with  , the vertical scale with  , the time scale with  , and denoting the non-dimensional forms of the variables by overhat (⋅) , the dimensionless form of the governing equations reads

of 27
We can now determine the form of these scales to reach the non-dimensional controlling parameter that characterizes the overall system behavior. The selection of l x as  is straightforward as l x is the dominant scale for the semi-infinite geometry and establishes the boundary condition in the model formulation. Equation B1 discloses the linear scaling of the diffusion term with the vertical scale  , and the scaling of the erosion term as  . Using the combination of parameters in the diffusion term, we write  = 2 . From the LHS of Equation B1, selecting  as ∕ would make this term without a non-dimensional parameter. Using these defined scales  ,  , and  , a dimensionless combination of model parameters and boundary conditions appear in the fluvial erosion term as Using numerical simulations of the complete landscape evolution model (LEM), we verified the validity of the presented dimensional analysis for different values of the parameters and exponents m and n in the LEM and their effects on the steady-state morphologies. We simulated a very long rectangular domain (l x = 100 m and l y = 5 × l x m) with fixed elevation boundary conditions at four sides such that only one typical length scale (l x ) remains relevant to compute  in Equation B5. The values of model parameters and exponents used for the nine cases presented here are reported in Table B1.
In Figure B1, the steady-state solutions are shown for the case of exponents m = 0.5 and n = 1, where the value of  does not depend on the uplift rate U. The value of U varies by four orders of magnitudes in these three cases, which does not alter the balance of creep to erosion component in the model ( = 50) and the solutions do not change with the presence of first-order valleys in all three cases. These steady-state solutions appear as re-scaled copies of each other, which is apparent when their non-dimensional forms (̂= 2 ) are noted. Figure B1. Steady-state topographies obtained from numerical simulations of the complete landscape evolution model in a rectangular domain with m = 0.5 and n = 1 and three sets of parameters that preserve the same  value equal to 50 (see Table B1). The dimensional form of the elevation field is shown in black-colored labels and the non-dimensional form is shown in blue-colored labels, indicating that the presented solutions are re-scaled copies of each other. The simulations were performed using the numerical algorithm developed in Anand Figure B2. Decreasing the value of uplift rate U by eight orders of magnitudes from panel a to panel c increases the intensity of the fluvial component compared to the diffusive creep in the LEM, which is quantified by the increased value of  from 5 to 500 in these cases. The steady-state landscape profile varies from an unchannelized case at  = 5 before the first instability with first-order valleys and few secondary branching at  = 50 to a complex pattern of ridges and valleys at  = 500 . Figure B3 displays three steady-state solutions for exponents m = 0.5 and n = 1.25, when the erosion term scales super-linearly with the vertical scale and increasing value of U results in more channelization. Uplift rate U is increased by twelve orders of magnitude in these three cases, which changes the non-dimensional  from 5, 500, to 5,000 and explains the obtained channelization cascade in these cases.  Table B1.

Appendix C: Eigenvalue Problem Formulation
Using  (= l x ),  where the prime (′) refers to the derivative with respect to . The expressions for coefficients are specified in Table C1 Finally, the boundary conditions for ̂(̂) in the changed reference variable arê Figure B3. Steady-state results obtained from the numerical simulations over a rectangular domain with m = 0.5 and n = 1.25.  value increases from 5 in panel a to 500 in panel b to 5,000 in panel c, which explains the increased channelization in the steady-state landscapes (please refer to Table B1 for the parameter values used in these results). The value of U increases by twelve orders of magnitude from panel a to panel c. The dimensional form of the elevation field z is shown in black-colored labels and the non-dimensional form is marked in blue-colored labels.
ANAND ET AL.

Appendix D: Weak Formulation and Galerkin Discretization
Equation C3 along with boundary conditions mentioned in Equations C4-C6 constitute a eigenvalue problem, which is solved here by transforming the final equation into an integral form (weak formulation). The dependency on has been omitted hereafter in the expressions for ease of notation.
A weak formulation is obtained by multiplying both sides of Equation C3 by a generic L 2 -test function v i (with i ∈ 1, N) and integrating over the interval (−1, 1) as ( where ( , )∶= ∫ 1 −1 (̂ ′ ) (̂ ′ ) ̂ ′ defines the inner product between two functions. The numerical approximation of inner products in the above equation can be computed by the interpolatory Legendre-Gauss quadrature formula, which approximates the integration of a generic function f in the domain [−1, 1] through the use of weights w k computed at discrete (Gauss-Lobatto) nodes as In the numerical solver developed for this work, and w k are computed using the algorithm provided by Swarztrauber (2003).
Based on previous works on spectral solutions of eigenvalue problems in shear flows (Shen, 1994), we seek a solution of ̂ in the form̂ where α j are the unknown coefficients of the linear expansion and u j reads

5(̂ )
̂0 2 (̂+ 1) 2 Note. The prime (′) refers to the derivative with respect to In the above expressions, L j represents the Legendre polynomial of degree j. So, u j (±1) = 0 for j ≥ 1 with u −1 (−1) = u 0 (−1) = 0. The additional functions u −1 and u 0 have been added to the basis to accommodate the non-vanishing boundary conditions.
The relationship between the coefficients for the imposed boundary conditions can be obtained as Applying this relation among the coefficients, the modified left-hand and right-hand matrix entries read ′ = 0 + −1 + ′ = 0 + −1 + ( ∈ [1 ]) . (D22) The algebraic system, ′ =̂′ , now consists of N equations in N unknowns (α j , j ∈ [1, N]), which can be solved as a generalized eigenvalue problem to compute the growth rate (̂) for different values of ̂ and  .

Data Availability Statement
The Python code developed for the linear stability analysis is available at https://github.com/ShashankAnand1996/ LEM_Stability_Analysis. The details of the solver used for the numerical simulations is described in Anand et al. (2020) and the well-commented Python code is available on GitHub https://github.com/ShashankAnand1996/ LEM (Anand, 2022).