#### Abstract

Linear and weakly nonlinear instability of shallow mixing layers is analysed in the present paper. It is assumed that the resistance force varies in the transverse direction. Linear stability problem is solved numerically using collocation method. It is shown that the increase in the ratio of the friction coefficients in the main channel to that in the floodplain has a stabilizing influence on the flow. The amplitude evolution equation for the most unstable mode (the complex Ginzburg–Landau equation) is derived from the shallow water equations under the rigid-lid assumption. Results of numerical calculations are presented.

#### 1. Introduction

Understanding the interaction between fast and slow fluid streams at river junctions and in compound channels is important for a proper description of mass and momentum transfer in shallow flows. In order to analyse the problem, hydraulic engineers apply depth-averaged shallow water equations [1] where either Chézy or Manning formulas are used to take into account the effect of bottom friction. These formulas contain a constant friction coefficient determined from empirical relationships [1] if the Reynolds number of the flow and surface roughness are given. There are cases, however, where the resistance force changes considerably in the transverse direction [2, 3]. One example of such a situation is flow in compound channels or rivers in case of floods. Hence, the analysis of instability characteristics of shallow mixing layers with variable friction is important from environmental point of view. In particular, contaminants and residues can accumulate in some parts of the flow due to instability. Thus, it is necessary to analyse factors affecting shallow flow instability and development of perturbations in a weakly nonlinear regime.

Three different approaches for the investigation of shallow water flows are suggested in [4]: experimental studies, numerical modelling, and stability analysis. Experimental data [5–7] indicate that bottom friction has a stabilizing influence on the flow. Temporal linear stability analysis of shallow mixing layers [8–12] shows that bed friction reduces the width of a mixing layer and stabilizes the flow. The development of perturbations can also be analysed from a spatial point of view [13–16]. Calculations show that bed friction reduces the spatial growth rate of perturbations.

There are other factors affecting stability characteristics of shallow mixing layers: (a) flow curvature, (b) presence of solid particles in a fluid stream, and (c) variable friction force in the transverse direction. These aspects of the linear stability problem are analysed in [13–16]. In particular, a stably curved mixing layer has a stabilizing influence on the flow, while an unstably curved mixing layer destabilises the flow. In addition, the presence of solid particles in the stream reduces the growth rate of perturbations.

Linear stability analysis is a useful tool for calculating critical values of the parameters characterizing the flow. However, linear theory cannot be used to analyse the development of perturbations in unstable regime. Assuming that the bed friction number is slightly smaller than the critical value (i.e., the flow is linearly unstable with a small growth rate), weakly nonlinear theories can be used to obtain an amplitude evolution equation for the most unstable mode [17–22]. The analysis in [17–21] shows that the amplitude evolution equation is the complex Ginzburg–Landau equation (GLE) under the rigid-lid assumption (free surface is treated as a rigid lid). As it is shown in [10], rigid-lid assumption works well for small Froude numbers. The Froude number is of order 0.2-0.3 in experiments [2] so that rigid-lid assumption is justified. The following cases of shallow flows are considered in [17–21]: (a) Rayleigh friction [19], (b) bottom friction modelled by the Chézy formula [21], and (c) generalization of the analysis in [21] for the case of slightly curved mixing layers and two-phase flows in [17–19]. Recent application of weakly nonlinear theory to shallow water flows without rigid-lid assumption [22] shows that the amplitude evolution equation for the most unstable mode is the complex Schrödinger equation. Since the solutions of the GLE vary from truly deterministic to almost chaotic [23] (depending on the values of the coefficients), in many cases, the GLE is used as a phenomenological equation for the analysis of spatiotemporal dynamics of complex flows. The coefficients of the equation are estimated using experimental data. Complex phenomena in fluid mechanics (see, e.g., [24]) can be modelled using the obtained equation.

In the present paper, we perform linear and weakly nonlinear analysis of shallow mixing layers under the assumption that the friction force changes in the transverse direction. Weakly nonlinear analysis shows that the amplitude evolution equation for the most unstable mode is the complex GLE. The equation is derived from the shallow water equations under the rigid-lid assumption. The coefficients of the equation are calculated in closed form using linear stability characteristics of the flow. Results of numerical calculations are presented. Practical applications of the proposed model are also discussed.

#### 2. Linear Stability Problem

Consider the system of shallow water equations under the rigid-lid assumption written in the form [25]where and are the depth-averaged velocity components in the and directions, respectively, is water depth, is the pressure, and is the friction coefficient depending on the transverse coordinate . The function is assumed to be of the formwhere is a constant and is an arbitrary differentiable “shape” function.

Since friction terms in (2)-(3) do not contain derivatives, the boundary conditions for the longitudinal velocity component are not needed (otherwise the problem would be overdetermined). Thus, the boundary conditions for the velocity components are essentially the same as in the case of inviscid flow (we need only zero normal velocity at the boundaries).

The boundary conditions are

The assumption of an unbounded layer in the transverse direction is widely used in the stability analysis of mixing layers and wakes (see, e.g., [26, 27]) and is adopted here.

Introducing the stream function by the relationsand eliminating the pressure from (1)–(3), we obtainwhere is the derivative of with respect to and .

Consider a perturbed solution to (7) of the formwhere is the base flow velocity. The role of the small parameter will be clarified later.

Substituting (8) into (7) and linearizing the resulting equation in the neighbourhood of the base flow , we obtainwhere

Using the method of normal modes, we represent the function in the formwhere is the amplitude of the normal perturbation and is the wavenumber.

Substituting (11) into (9), we obtainwhere is the bed friction number and is the half-width of the mixing layer.

The boundary conditions follow from (5) and have the form

Complex eigenvalues determine the linear stability of the base flow. The flow is linearly stable if all and unstable if at least one . Problems (12) and (13) are solved numerically by means of the collocation method based on the Chebyshev polynomials [15]. In the classical hydrodynamic stability theory (see, e.g., [26, 27]), the base flow is obtained as a simple one-dimensional steady solution of the equations of motion (e.g., a plane Poiseuille flow is obtained as a steady one-dimensional solution of the Navier–Stokes equations). It is impossible to find a simple analytical solution of the form for shallow water equations (1)–(3) due to the empirical character of the Chézy formula in (2) and (3). The paper [28] gives an example of a longitudinal velocity distribution in the form of a generalized power law. This profile does not have an inflection point (recall that the presence of an inflection point is the necessary condition for instability of inviscid flows [26]). The obtained formula for the longitudinal dispersion coefficient in [28] is based on the assumption of straight symmetrical channel and uniform flow. The resulting formula is multiplied by a “catch-all” revision coefficient taking into account many kinds of nonuniformities in natural rivers. However, this adjustment does not change the shape of the longitudinal velocity distribution. The mixing layer case (the case considered in our paper) has a completely different longitudinal velocity distribution [2]. The presence of a porous layer with much higher resistance than in the main channel results in a highly nonuniform velocity distribution resembling a hyperbolic tangent function with some asymmetry (however, the asymmetric profile is obtained after all the instabilities have been taken into account). The model profile has an inflection point. Thus, instabilities can occur. The hyperbolic tangent profile is widely used in stability analyses of shallow flows [8–12]. The goal of our investigation is to analyse instabilities of the flow, obtain the linear stability characteristics, and consider the development of perturbations in a weakly nonlinear regime.

Assuming that the velocity of undisturbed flow is and at and , respectively, the function is chosen in the formwhere all the dimensional quantities are represented with tildes. Using as the velocity scale and as the length scale, we rewrite (14) in the dimensionless formwhere is the velocity ratio. All calculations in the paper are performed for the case . The base flow velocity distribution is given in Figure 1. The domain of the flow is defined as follows: the transverse coordinate varies from to . A porous layer occupies the region , while the region corresponds to the main channel. Away from the shear layer, the base flow velocity approaches the undisturbed values and , respectively. We assume that the volume fraction is small (in such a case, as it is shown in [2] shallow water equations of the form (1)–(3) can be used to analyse the problem).

The shape function is selected in the formwhere is the ratio of the friction coefficients in the main channel and floodplain and is the factor reflecting the shape of the transitional zone from floodplain to the main channel.

The following arguments are used to select in (16). The friction coefficients for the undisturbed flow (away from the shear layer) in the main channel and floodplain are denoted by and , respectively. A stronger resistance force results in a smaller base flow velocity so that (15) and (16) are consistent. In addition, we would like to remove discontinuity in the friction force used in [2] and consider a more realistic case of a continuous resistance changing with respect to the transverse coordinate. Examples of marginal stability curves are shown in Figure 2 for two different values of . The flow is stable above the curves and unstable below the curves. Small circles in Figures 2–4 indicate the calculated points, while the solid lines represent the interpolating function. The flow becomes more unstable as increases. This is related to the fact that larger result in a steeper transition from the value of the friction coefficient in the main channel to the one in floodplain (16). The critical values of the parameters and , namely, and , respectively, correspond to the maximum of the marginal stability curves. Figure 3 plots the critical bed friction numbers versus . It is seen from the graph that the flow becomes more unstable as the parameter decreases (i.e., if the ratio of the friction coefficients in the main channel to that in floodplain becomes smaller). Thus, the increase in the ratio of the friction coefficients stabilises the flow.

#### 3. Weakly Nonlinear Analysis

Suppose that the bed friction number is slightly smaller than the critical value :

Thus, the small parameter in asymptotic expansion (8) can be defined by the formula

This relationship shows how close the bed friction number is to the corresponding critical value. Since the value of can be calculated for any set of experimental conditions, the corresponding value of also can be computed.

Using the method of multiple scales and following [29], we introduce the slow time and longitudinal variable by the relationswhere is the group velocity. Substituting (17) and (19) into (8) and collecting the terms containing , we obtain the system of equationswhere the operator is defined in (10). The function in (21) has the form

Equation (20) coincides with (9) and corresponds to the linearized problem where the function is assumed to be of the form (11). Suppose thatwhere and are, respectively, the critical wavenumber and speed of the most unstable perturbation in accordance with the linear theory, is the unknown amplitude function, and the abbreviation c.c. means complex conjugate. The goal of the weakly nonlinear analysis is to derive the amplitude evolution equation for the function . Derivation of the amplitude equation can be found in the Appendix. The Ginzburg–Landau equation has the form

The sign of the real part of (known as the Landau constant in the literature) plays an important role. The Landau constant has the “wrong sign” in [26]. Thus, finite amplitude saturation is not possible, and higher order terms (with respect to ) quickly become important so that (25) can be used for a very short time (in other words, practical application of (25) is very limited). It is shown in [17–21] (in contrast to [29]) that, for shallow water flows, the Landau constant in (25) has the “right sign” so that (25) can be used (and is successfully used in [17–21]) in order to describe some important features of shallow wake flows.

Equation (25) can be written in a more convenient form using the substitutions

Thus,

Some closed-form solutions of (27) are known in the literature [23]. One of the simplest solutions is the solution of the formwhere

Since (28) represents a periodic solution in both space and time, it is important to investigate the linear stability of such a plane wave. If solution (28) is stable, then periodic wave-like solutions can also be observable in experiments.

Assuming thatlinear stability of (28) can be analysed [23].

Substituting (30) into (27), we obtain the equation for . For the case of small , the stability condition has the formprovided that satisfies the inequality

Condition (31) is known as the Benjamin–Feir stability condition. If (31) is not satisfied, then plane wave solutions of the form (28) are unstable (and, therefore, are not observable in experiments).

#### 4. Numerical Results

The complex Ginzburg-Landau equation as the amplitude evolution equation for the most unstable mode is analysed in the previous section using asymptotic expansion (8). Since it is not known in general whether (8) is convergent, it would be nice to have a criterion which can be used to convince us that the Ginzburg–Landau model can adequately represent the dynamics of a fully nonlinear model at least at the initial stage of transition period when the base flow becomes linearly unstable. One such criterion is given in [31]. The criterion is as follows: if the growth rates of the unstable perturbation can be well approximated by a parabola in the whole range of unstable wavenumbers, then the GLE can be used to analyse the dynamics of the flow (at least in the beginning of the nonlinear regime). In order to test this assertion, we computed the growth rates for the range of unstable wavenumbers for the following values of the parameters of the problem: .

The results of calculations are shown in Figure 4, where *S* < *S*_{c} and the dots represent the calculated growth rates for three values of . The solid lines represent the best-fit parabolas. As can be seen from the graph, the calculated growth rates are well approximated by parabolas in all regions of unstable wavenumbers. Thus, we conclude that the GLE can be successfully used to analyse the dynamics of the flow above the threshold.

The results of the numerical computations of the linear stability characteristics and the coefficients of the GLE are shown in Table 1.

Table 2 shows numerical values of the coefficients and (26) for different values of . As can be seen from Table 2, condition (31) is satisfied for all cases considered. Hence, periodic solutions of the GLE are stable and can be observed in experiments. The Landau constant (the real part of ) is positive. Thus, finite amplitude saturation is possible.

Numerical solutions of the GLE (27) are presented below for different values of the parameters and and different initial conditions. The problem is formulated as follows: find the solution of (27) for the given boundary conditionsand the initial condition

Method of lines implemented in Mathematica is used for the numerical solution of problems (27), (33), and (34). The value of is used for all simulations below.

The first computation is performed for the case *c*_{1} = 1.3294 and *c*_{2} = 0.1252. The values of these parameters are taken from Table 2. The function in (34) is assumed to be a small random noise of order 0.01. The results are shown in Figure 5. Since we are interested in the analysis of spatiotemporal development of a perturbation, the modulus of the amplitude is plotted versus the longitudinal coordinate and time . The parameters of the problem satisfy (31) and (32) (in other words, are in the region of stability). Thus, the modulus of the amplitude reaches a constant value.

The second set of computations is performed for the case where *q* = 0.5. The results are shown in Figure 6. The value of in this case corresponds to the stable range of wavenumbers (both conditions (31) and (32) are satisfied in this case). As a result, perturbations do not propagate throughout the domain and finite amplitude saturation occurs.

Finally, we consider the case where the Benjamin–Feir stability condition (31) is not satisfied. The values of the parameters are taken from [21]: *c*_{1} = −0.799564 and *c*_{2} = 2.189654 (these parameters correspond to weakly nonlinear analysis of wake flows). Random noise of order 0.01 is used as the initial condition. The results are shown in Figure 7. As can be seen from the figure, stabilization of the amplitude does not occur in this case. Moreover, instability is spreading throughout the whole domain (see the parts of the graph for large times). This fact is confirmed by experiments (see, e.g., discussion of wake flow instability in [21] and references therein)—periodical structures do not appear in experiments in the weakly nonlinear regime.

These examples illustrate the well-known fact that the Ginzburg–Landau model is quite rich in terms of different solutions. Illustrative computations in Figures 5–7 show that both initial conditions and the values of the coefficients are responsible for spatiotemporal dynamics of the amplitude.

Direct comparison of linear and weakly nonlinear instability results with field data is not an easy task. In order to compare directly the results of weakly nonlinear theory, one has to use carefully designed experiments (or numerical experiments). Available experimental data (see, e.g., [2, 3]) give the velocity distribution in the longitudinal direction which is averaged over sufficiently large period of time (in this case, nonlinear effects have been already taken into account). The goal of weakly nonlinear theory is to describe the initial stage of the development of instability. Using weakly nonlinear calculations, we can answer the following questions: Should we expect finite amplitude saturation? Would a periodic solution be stable (and, therefore, observable in experiments)? Illustrative calculations for the solution of the boundary value problem containing GLE are consistent with experimental observations. For example, if the Benjamin–Feir instability condition (31) is not satisfied for wake flows, then periodic solutions of the GLE are unstable (and are not observable in experiments). This fact is supported by experimental observations for wake flows (see the references in [21]). On the other hand, our calculations show that condition (31) is satisfied for shallow mixing layers. Thus, periodic solutions of the GLE are stable (this fact is supported by experimental observations in [2, 3]).

#### 5. Conclusions

Linear stability of shallow mixing layers is investigated under the assumption that the resistance force varies in the transverse direction. Stability of the base flow depends on the ratio of the friction coefficients in the main channel to that in the floodplain. Weakly nonlinear theory is used to derive the amplitude evolution equation for the most unstable mode. It is shown that the amplitude evolution equation is the Ginzburg-Landau equation with complex coefficients. The coefficients of the equation are obtained in closed form available for calculations. Numerical results show that plane wave solutions of the GLE are stable.

#### Appendix

Substituting (24) into (21) and taking into account the structure of (23), we seek the solution to (21) in the formwhere is the complex conjugate of and are unknown functions of .

Collecting the terms proportional to , we obtain the equation for :

The boundary conditions have the form

Similarly, collecting the terms proportional to on the left-hand and right-hand sides of (21), we obtain the equation for the function :with the boundary conditions

Finally, collecting the terms proportional to , we obtainwith the boundary conditions

The left-hand side of (A.4) coincides with the left-hand side of (12). Thus, we have to solve nonhomogeneous equation (A.4) under the assumption that the corresponding homogeneous equation has a nontrivial solution. The Fredholm’s alternative [30] has to be used to guarantee the solvability of (A.4).

Let

The adjoint operator and adjoint eigenfunction are defined by the relation

The left-hand side of (A.9) is equal to zero since . Thus, the adjoint equation is defined by the formula

Integrating the left-hand side of (A.9) by parts and using the boundary conditions , we obtain

Thus, the adjoint equation has the form

The boundary conditions are

The adjoint eigenfunction is the solution to problems (A.12) and (A.13). In accordance with the Fredholm’s alternative, (A.4) is solvable if and only if the right-hand side of (A.4) is orthogonal to all eigenfunctions of the corresponding adjoint problems (A.12) and (A.13). Applying the solvability condition to (A.4), we obtain

Equation (A.14) defines the group velocity:where

Solving three boundary value problems (A.2)-(A.3), (A.4)-(A.5), and (A.6)-(A.7) numerically, we obtain the functions , and . The function (the second-order correction) is given by (A.1).

Let us consider the solution to (22) at the third order in . The function in (22) has the form

Equation (22) also has a solution if and only if the right-hand side of (A.18) is orthogonal to all eigenfunctions of the corresponding homogeneous adjoint problem. Applying the solvability condition to (22), we obtain

Equation (A.19) is the complex Ginzburg–Landau equation with the complex coefficients , and :where , and are complex coefficients, which can be computed using linearized characteristics of the flow.

The coefficients , , and are given by

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was partially supported by the Grant 623/2014 of the Latvian Council of Science.