text:
# Integral methods for shallow free-surface flows with separation

## 1 Introduction

In this paper we develop a simple quantitative method to describe flows with a free surface which can undergo large distortions. Our method is capable of handling flows whose velocity profile may become far from parabolic — even including separation and regions of reverse flow. We are concerned with the case when the fluid layer is thin. For low Reynolds number flows the lubrication approximation can be used with great success (see e.g. ). For high Reynolds number flows without separation an inviscid approximation and the shallow water equations  are widely used. For moderate Reynolds numbers where these limiting approximations are invalid it is important to take both inertial and viscous effects into account in a consistent way, and yet one would like to keep the model simple enough to be tractable. In this paper we show that integral methods, like the ones developed by von Karman, can handle a class of such problems successfully. To be concrete we develop the method in the context of two physical examples: the circular hydraulic jump and the flow down an inclined plane. Both geometries support jump- or kink-like solutions with abrupt changes in the surface shape and internal velocity profiles. Analytical solutions for such flows are extremely difficult to obtain, and simple approximate theories that capture the phenomena are invaluable.

The two flows are studied in separate sections, and an introduction is provided in the beginning of Secs. 2 and 3, respectively. In Sec. 2 we develop the theory for the circular hydraulic jump. We first study the boundary layer approximation to the full Navier-Stokes equations, and reduce it to a simple set of equations by averaging over the thickness. Stationary solutions are obtained by solving a two-point boundary value problem for a system of only two ordinary differential equations. The solution is compared to previous experiments, showing good agreement. Taking advantage of the simplicity of the reduced equations, it is possible to obtain analytic approximations for the stationary solution. Two “outer” solutions connected by an “inner” transition region are studied separately and we obtain a relationship analogous to the shock condition in the classical shock theory, but within our viscous model.

The flow down an inclined plane is then studied in Sec. 3. We use the same strategy as in Sec. 2 to derive a simple model for the two-dimensional flow. One family of solutions found in this model is kink-like traveling wave solutions that occur, e.g., when the flow rate is suddenly changed. Their velocity profiles along the inclined plane are found to stay close to parabolic even when a variable profile is assumed. There is another family of solutions with a sudden change in the surface that would correspond to the circular hydraulic jump in case of the radial geometry. These solutions can be interpreted as the stationary hydraulic jump, created behind a sluice gate in a river, even though turbulence is not included in the model. The flow downstream of the jump approaches a simple stationary flow, but the flow upstream is an expanding flow with a linear growth in thickness. The velocity profile departs considerably from parabolic near the jump.

It is not easy to analyze the stability of the solutions with jumps obtained in Secs. 2 and 3, even in the linear geometry. Instead, in Sec. 4, we include time-dependence in the models and study the dispersion relation for the stationary flow with constant thickness. A well-established concept in the inviscid theory is to classify flows as super- and subcritical when the thickness is small and large, respectively. They do not have obvious counterparts, however, when viscosity is included. By looking carefully at the dispersion relation in the long and medium wave regime, we can classify the stationary flows into these two categories in our viscous model. The model shows spurious divergencies in the short wavelength region which we do not know how to overcome at present. This makes the model unsuited for direct time-dependent simulations. A short paper describing some of the main results has appeared earlier .

## 2 The circular hydraulic jump

### 2.1 Introduction to the problem

When a jet of fluid hits a flat horizontal surface, the fluid spreads out radially in a thin, rapidly flowing layer. At a certain distance from the jet a sudden thickening of the flow takes place, which is called the circular hydraulic jump. This is commonly seen, e.g., in the kitchen sink, but it is also important as a coating flow and in jet-cooling of a heated surface . In these practical flows with typically high Reynolds numbers, disturbances often make the jump non-stationary and distorted. In controlled laboratory experiments corresponding to a more moderate Reynolds number, an apparently stationary, radially-symmetric flow can be achieved. Such experiments was carried out by C. Ellegaard et al. and the results have been published elsewhere . We thank them for providing us with data and pictures. A schematic view and a video image of the circular jump are shown in Fig. 1.

In these experiments the hydraulic jump is formed on a flat disc with a circular rim. The rim height $`d`$ can be varied, and is an important control parameter. Since the rim is located far from the impinging jet with the diameter of the disc around 36cm, it does not affect the jump except that it changes the height of the fluid layer $`h_{ext}`$ exterior to the jump. The jump still forms even when $`d=0`$, but a larger $`d`$ makes $`h_{ext}`$ larger and, therefore, the jump stronger. Typically, $`h_{ext}`$ exceeds $`d`$ by 1-2mm. The surface profiles for varying $`d`$ are shown in Fig. 2. An interesting transition in the flow structure has been observed  as $`d`$ is varied. For $`d=0`$, it was noticed before  that the jump contains an eddy on the bottom, called a separation bubble, whose inner edge is located very close to the position of the abrupt change on the surface, as illustrated in Fig. 3(a). Such a hydraulic jump is referred to as a type I jump. While $`d`$ remains small, this jump is stable, but as $`d`$ is increased further, a wave-breaking transition occurs  which results in a new state of the flow. In this type II state, the flow has an additional eddy, called a roller or a surfing wave, just under the surface as shown in Fig. 3(b).111 If $`d`$ is increased even further, the jump “closes” as seen in Fig. 2.  This state resembles a broken wave in the ocean, but is still apparently laminar. On reducing $`d`$, the type I pattern reappears, and there is almost no hysteresis associated with this transition. The transition from type I to II often leads also to breaking of the radial symmetry. An intriguing set of polygonal jumps  are created rather than the circular one. In this paper we shall concentrate on the type I flow which already poses considerable difficulties. We hope to be able to generalize our approach in the future to be able to handle the transition to the type II flow.

Considering how simple and common the circular hydraulic jump appears to be, it is surprising that a satisfactory systematic theory does not exist. The approach considered as “the standard” for the study of hydraulic jumps is to combine the inviscid shallow water equation with Rayleigh’s shocks . In the beginning of the century Lord Rayleigh treated  a discontinuity in a one-dimensional linear flow geometry. Such a structure is usually called a river bore if it is moving and a hydraulic jump if it is stationary and is created due to, e.g., variations in the river bed. His approach was based upon the analogy between the shallow water theory and gas theory . He assumed that, across such a shock, the mass and momentum flux are conserved but not the energy flux.

In a coordinate system moving with the shock, the flow velocity $`v_1`$ and height $`h_1`$ upstream of the jump as well as $`v_2`$ and $`h_2`$ downstream of the jump are taken to be positive constant values. Then, conservation of mass flux $`Q`$ across the jump is given by

$$v_1h_1=v_2h_2=Q.$$
(1)
Conservation of momentum flux is

$$h_1\left(v_1^2+\frac{1}{2}gh_1\right)=h_2\left(v_2^2+\frac{1}{2}gh_2\right).$$
(2)
These shock conditions lead to the relation

$$\frac{h_2}{h_1}=\frac{1}{2}(1+\sqrt{1+8F_1^2})=\frac{2}{1+\sqrt{1+8F_2^2}}$$
(3)
where $`F_1=\sqrt{v_1^2/gh_1}=(h_c/h_1)^{3/2}`$ is the upstream Froude number, $`F_2=\sqrt{v_2^2/gh_2}=(h_c/h_2)^{3/2}`$ the downstream Froude number, and

$$h_c=(Q^2/g)^{1/3}$$
(4)
is called the critical height. It is easy to see that $`h_c`$ is always between $`h_1`$ and $`h_2`$, and that $`F_1>1>F_2`$ if $`h_1<h_c<h_2`$, and $`F_1<1<F_2`$ if $`h_1>h_c>h_2`$. In other words the jump connects a supercritical flow with $`F>1`$ on the shallower side ($`h<h_c`$) to a subcritical flow with $`F<1`$ on the deeper side ($`h>h_c`$). Since the Froude number measures the ratio of the fluid velocity $`v`$ and the velocity of linear surface waves $`\sqrt{gh}`$, it means that, in the moving frame, the flow moves more rapidly than the surface waves on the shallower side, but moves slower on the deeper side — in a precise analogy with the gas theory . Further, it is found that the upstream $`h_1`$ must be supercritical by considering the change in the energy flux across the jump :

$$Q_{e2}Q_{e1}=\frac{gQ(h_2h_1)^3}{8\pi h_1h_2}$$
(5)
where $`Q_e`$ denotes the energy flux. Since the energy must be dissipated through the jump, i.e. $`Q_{e2}Q_{e1}<0`$, rather than generated, it is required that $`h_1<h_2`$. The origin of the dissipation is usually attributed to the turbulent motions at the discontinuity and surface waves carrying energy away from it.

It is possible to apply this theory, combined with an assumption of the potential flow, for describing the circular hydraulic jump. However, it leads to incorrect estimates  of the radius of the jump $`R_j`$. Most notably, $`R_j`$ is predicted to be sensitive to the radius of the impinging jet which should be greatly influenced by radius and height of the inlet nozzle where liquid comes out. In experiments  such a strong tendency was not observed. Instead, it has been found that $`R_j`$ scales with the flow rate $`Q`$ with a certain power, and it supports a model in which viscosity plays an important role. Watson  constructed a model of the flow consisting of the inviscid and viscid regimes, and solved the viscid part assuming a similarity profile. By connecting to the specified external height $`h_{ext}`$ via a Rayleigh shock, he obtained a prediction for the radius of the jump which compares favorably with the measurement , as we explain in Sec. 2.5. In his model the viscous layer starts from the stagnation point at $`r=0`$ on the plate and quickly reaches the surface at a small $`r`$. There is a fairly long stretch from this $`r`$ to $`R_j`$ in which the flow is fully viscous.222 This assumption is confirmed by recent laser-doppler measurements of the velocity profile before the jump . Thus, the assumption made by  that the jump occurs at the point where the growing viscous layer touches the surface and the flow becomes fully developed, is incorrect.  Thus, one could neglect the inviscid region and assume a fully viscid flow everywhere in order to derive a simpler model. This assumption was made by Kurihara  and Tani  who started from the boundary layer equations developed by Prandtl . They took an average of the equations over the thickness, assuming also a similarity velocity profile. It resulted in a single ordinary differential equation for the stationary jump. This theory was elaborated in  who realized that the flow outside the jump would naturally lead to a singularity at a large $`r`$. By identifying this singularity with the outflow over the rim of the plate, the flow outside the jump could be uniquely specified. By introducing a Rayleigh shock, the jump radius and its parameter dependence was calculated and compared to measurements. The model predicted the observed $`R_j`$ reasonably well, as we review in Secs. 2.2–2.5.

Obviously, treating the jump as a discontinuity provides us with no information on the internal structure of the jump region such as the type I to II transition of the flow patterns. It also seems inconsistent to assume a Rayleigh shock when viscous loss occurs in the whole domain. Why do we assume an extra energy loss at the “jump” where the flow is stationary and apparently laminar? It seems possible to attribute the energy dissipation entirely to laminar viscous forces, and to construct a viscous theory which produces a smooth but kink-like surface shape without the need for a discontinuity. Nevertheless, such a description must overcome a difficulty arising from the Goldstein-type singularity  of the boundary layer equations in the vicinity of separation points. This singularity is thought to be an artificial one created by truncation of higher derivatives from the Navier-Stokes equations. It also arises in the “usual” boundary layer situation where a high Reynolds number flow passes a body, e.g., a wing. In such cases inviscid-viscid interaction is taken into account in order to resolve the singularity in a technique called the inverse method . In our situation, however, there is no inviscid flow outside the layer. In Secs. 2.6–2.7 we propose a way to resolve the trouble in the following manner. We first include an additional degree of freedom in the velocity profile to make it non-self-similar, just like in the Karman-Pohlhausen method  for the usual boundary layer theory. To describe the evolution, in $`r`$, of this free parameter, we couple the layer thickness to the pressure by assuming hydrostatic pressure. This serves as an alternative to the inverse method in the absence of a potential external flow. The resulting model for a stationary solution is two coupled ordinary differential equations, and reproduces the type I flow with a separation bubble — the one shown in Fig. 3(a). Comparison with the experiment is made in Sec. 2.7. It is possible to approximate analytically the stationary solution found in the model. In Sec. 2.8 the analysis is presented separately for the regions before and after the jump (i.e. two “outer” solutions) and the “inner” solution inside the jump region. An interesting observation on the inner solution is that a formal parameter $`\mu `$ can be introduced so that Rayleigh’s shock condition is recovered in the limit $`\mu 0`$.

### 2.2 The full model

We write down the complete model to describe the circular hydraulic jump under the assumption that the flow is laminar and radially symmetric without any angular velocity component. We take the radial and vertical coordinates $`\stackrel{~}{r}`$ and $`\stackrel{~}{z}`$, and denote the velocity components by $`\stackrel{~}{u}`$ and $`\stackrel{~}{w}`$, respectively.333 We use tildes for the dimensional variables, dependent or independent. Dimensionless variables will be expressed by the same symbols but without tildes. In figures, however, we do not use tildes for simplicity. The governing equations are the continuity equation

$$\stackrel{~}{u}_{\stackrel{~}{r}}+\frac{\stackrel{~}{u}}{\stackrel{~}{r}}+\stackrel{~}{w}_{\stackrel{~}{z}}=0$$
(6)
and the Navier-Stokes equations:

$$\begin{array}{c}\stackrel{~}{u}_{\stackrel{~}{t}}+\stackrel{~}{u}\stackrel{~}{u}_{\stackrel{~}{r}}+\stackrel{~}{w}\stackrel{~}{u}_{\stackrel{~}{z}}=\frac{1}{\rho }\stackrel{~}{p}_{\stackrel{~}{r}}+\nu \left(\stackrel{~}{u}_{\stackrel{~}{r}\stackrel{~}{r}}+\frac{1}{\stackrel{~}{r}}\stackrel{~}{u}_{\stackrel{~}{r}}\frac{\stackrel{~}{u}}{\stackrel{~}{r}^2}+\stackrel{~}{u}_{\stackrel{~}{z}\stackrel{~}{z}}\right)\hfill \\ \stackrel{~}{w}_{\stackrel{~}{t}}+\stackrel{~}{u}\stackrel{~}{w}_{\stackrel{~}{r}}+\stackrel{~}{w}\stackrel{~}{w}_{\stackrel{~}{z}}=\frac{1}{\rho }\stackrel{~}{p}_{\stackrel{~}{z}}g+\nu \left(\stackrel{~}{w}_{\stackrel{~}{r}\stackrel{~}{r}}+\frac{1}{\stackrel{~}{r}}\stackrel{~}{w}_{\stackrel{~}{r}}+\stackrel{~}{w}_{\stackrel{~}{z}\stackrel{~}{z}}\right)\hfill \end{array}$$
(7)
where subscripts denote partial differentiations such as $`\stackrel{~}{u}_{\stackrel{~}{t}}=\stackrel{~}{u}/\stackrel{~}{t}`$. For the boundary conditions we impose no-slip on the bottom:

$$\stackrel{~}{u}(\stackrel{~}{z}=0)=\stackrel{~}{w}(\stackrel{~}{z}=0)=0.$$
(8)
The dynamic boundary conditions on the free surface $`\stackrel{~}{z}=\stackrel{~}{h}(\stackrel{~}{t},\stackrel{~}{r})`$ are

$$\begin{array}{c}\stackrel{~}{p}\frac{2\rho \nu }{1+\stackrel{~}{h}_{\stackrel{~}{r}}^2}\left[\stackrel{~}{h}_{\stackrel{~}{r}}^2\stackrel{~}{u}_{\stackrel{~}{r}}+\stackrel{~}{w}_{\stackrel{~}{z}}2\stackrel{~}{h}_{\stackrel{~}{r}}\left(\stackrel{~}{w}_{\stackrel{~}{r}}+\stackrel{~}{u}_{\stackrel{~}{z}}\right)\right]|_{\stackrel{~}{z}=\stackrel{~}{h}}=\sigma \stackrel{~}{k}\hfill \\ \nu \left[\left(\stackrel{~}{h}_{\stackrel{~}{r}}^21\right)\left(\stackrel{~}{w}_{\stackrel{~}{r}}+\stackrel{~}{u}_{\stackrel{~}{z}}\right)2\stackrel{~}{h}_{\stackrel{~}{r}}\left(\stackrel{~}{u}_{\stackrel{~}{r}}\stackrel{~}{w}_{\stackrel{~}{z}}\right)\right]|_{\stackrel{~}{z}=\stackrel{~}{h}}=0\hfill \end{array}$$
(9)
where $`\sigma `$ is the coefficient of surface tension and $`\stackrel{~}{k}`$ is the mean local curvature of the free surface. We also need to satisfy the kinematic boundary condition on the free surface:

$$\stackrel{~}{h}_{\stackrel{~}{t}}+\stackrel{~}{u}\stackrel{~}{h}_{\stackrel{~}{r}}=\stackrel{~}{w}\text{on }\stackrel{~}{z}=\stackrel{~}{h}(\stackrel{~}{t},\stackrel{~}{r})\text{.}$$
(10)
We are mostly interested in stationary solutions in this section. When the flow is stationary, we may integrate (6) over $`\stackrel{~}{z}`$ from 0 to $`\stackrel{~}{h}`$, and use (10) to obtain

$$\stackrel{~}{r}_0^{\stackrel{~}{h}(\stackrel{~}{r})}\stackrel{~}{u}(\stackrel{~}{r},\stackrel{~}{z})\text{d}\stackrel{~}{z}=q=\frac{Q}{2\pi }.$$
(11)
This quantity, the total mass flux $`Q`$ or the mass flux per angle $`q`$, is a constant, given as a parameter in the experiment.

### 2.3 Boundary layer approximation

Since it is a formidable task to treat the full model as it stands, some simplifications need to be made. As explained in Sec. 1, the Reynolds number for the flow of the circular hydraulic jump is too large to justify the lubrication approximation, but is not large enough to use the inviscid approximation. Fortunately, the flow is “thin,” i.e. runs predominantly horizontally along the plate. Truncation of the full model by the boundary layer approximation is quite natural in such a situation, and has indeed been used in previous literature . In the boundary layer approximation pressure, viscous, and inertial terms in (7) are all assumed to be of the same order, but there are only a few dominant terms in each group. For instance, a viscous term $`\nu \stackrel{~}{u}_{\stackrel{~}{r}\stackrel{~}{r}}`$ is assumed to be negligible compared to $`\nu \stackrel{~}{u}_{\stackrel{~}{z}\stackrel{~}{z}}`$. The dominant terms in the first equation in (7) are determined in the usual manner: $`\stackrel{~}{u}_{\stackrel{~}{t}}`$ (if time-dependent), inertia terms $`\stackrel{~}{u}\stackrel{~}{u}_{\stackrel{~}{r}}`$ and $`\stackrel{~}{w}\stackrel{~}{u}_{\stackrel{~}{z}}`$, the pressure term $`\stackrel{~}{p}_{\stackrel{~}{r}}/\rho `$, and the dominant viscous term $`\nu \stackrel{~}{u}_{\stackrel{~}{z}\stackrel{~}{z}}`$. Similarly, from the second equation in (7) we assume the dominant balance between $`\stackrel{~}{p}_{\stackrel{~}{z}}/\rho `$ and $`g`$. Here, unlike the usual boundary layer theory, we have taken into account the effect of gravity. This will couple the surface height $`h`$ to the pressure, and will later turn out to be crucial for removing the singularities of the boundary layer approximation.

If we denote the characteristic radius and height by $`r_{}`$ and $`z_{}`$, respectively, then the second dominant balance requires the characteristic pressure to be $`\rho gz_{}`$. Then, the first balance relation requires

$$\frac{u_{}}{t_{}}=\frac{u_{}^2}{r_{}}=\frac{u_{}w_{}}{z_{}}=\frac{\rho gz_{}}{\rho r_{}}=\frac{\nu u_{}}{z_{}^2}$$
(12)
where $`u_{}`$ and $`w_{}`$ are typical radial and vertical velocities, respectively, and $`t_{}`$ is the characteristic time scale. The mass flux relation (11) requires that

$$u_{}r_{}z_{}=q.$$
(13)
while the continuity equation (6) requires

$$\frac{u_{}}{r_{}}=\frac{w_{}}{z_{}}.$$
(14)
Solving (12), (13), and (14) uniquely determines the characteristic scales:

$$\begin{array}{c}r_{}=\left(q^5\nu ^3g^1\right)^{1/8}2.7\text{[cm]},\hfill \\ z_{}=\left(q\nu g^1\right)^{1/4}1.5\text{[mm]},\hfill \\ u_{}=\left(q\nu g^3\right)^{1/8}12\text{[cm/s]},\hfill \\ w_{}=\left(q^1\nu ^3g\right)^{1/4}6.7\text{[mm/s]},\hfill \\ t_{}=\left(q\nu ^1g^1\right)^{1/2}0.22\text{[s]}\hfill \end{array}$$
(15)
where the estimated values correspond to a typical set of parameters used in the experiments: $`\nu 0.1`$ cm2/s (for mixture of ethylene-glycol and water) and $`Q30`$ cm3/s, i.e., $`q5`$ cm3/s. The values for $`r_{}`$ and $`z_{}`$ correspond well to a typical jump radius and fluid thickness in the experiments. Also, the predicted scaling can be experimentally tested by, for instance, measuring the dependence of the jump radius by changing the parameters such as $`q`$. In  evidence of the scaling and validity of the underlying assumption was given.

We now use the characteristic scales (15), together with the pressure scale $`p_{}=\rho u_{}^2`$, to non-dimensionalize the full equations. From (7), we obtain

$$\begin{array}{c}u_t+uu_r+wu_z=p_r+u_{zz}+ϵ^2\left(u_{rr}+\frac{1}{r}u_r\frac{u}{r^2}\right)\hfill \\ ϵ^2\left(w_t+uw_r+ww_z\right)=p_z1+ϵ^2w_{zz}+ϵ^4\left(w_{rr}+\frac{1}{r}w_r\right),\hfill \end{array}$$
(16)
where

$$ϵ=z_{}/r_{}=\left(q^3\nu ^5g^1\right)^{1/8}.$$
(17)
Since $`ϵ=z_{}/r_{}=0.05`$ for the typical parameter values above, the assumption that the flow is “thin” is well satisfied, and we shall drop the terms of order $`ϵ^2`$ and higher in the equations (16). We also focus on stationary solutions in the rest of the section, and thus we obtain the simplified equations of motion:

$$\begin{array}{c}uu_r+wu_z=p_r+u_{zz}\hfill \\ 0=p_z1.\hfill \end{array}$$
(18)
Correspondingly, within the error of $`O(ϵ^2)`$, the dynamic boundary conditions (9) are just

$$\begin{array}{c}p|_{z=h}=Wh_{rr}\hfill \\ u_z|_{z=h}=0.\hfill \end{array}$$
(19)
Here we have introduced the Weber number

$$W=\frac{\sigma z_{}}{\rho u_{}^2r_{}^2}=\frac{\sigma }{\rho gr_{}^2}=\frac{\mathrm{}^2}{2r_{}^2}=\sigma \rho ^1(q^5\nu ^3g^3)^{1/4}.$$
(20)
where $`\mathrm{}=(2\sigma /(g\rho ))^{1/2}`$ is the capillary length. For the parameter values above together with $`\sigma 70`$\[dyn/cm\] (maximum), we estimate that $`W0.01`$ and $`\mathrm{}3.8`$\[mm\]. Since $`W`$ is small, we neglect it in the study of stationary states.444 However, the term influences dispersion of short waves, so should be included in the stability analysis of stationary states, possibly together with the neglected terms of $`O(ϵ^2)`$ and higher in (18).,555 The Reynolds number, defined as $`R=u_{}z_{}/\nu =(q^3\nu ^5g)^{1/8}18`$. The Reynolds number at the nozzle outlet is much higher, but it becomes moderate near the jump. The second equation of (18) and the first condition of (19) with $`W`$ set to zero yield hydrostatic pressure:

$$p(r,z)=h(r)z.$$
(21)
Combining (18) and (21), we obtain the stationary boundary layer equations:

$$uu_r+wu_z=h^{}+u_{zz},$$
(22)
where the prime denotes the derivative with respect to $`r`$. This is supplemented by the dimensionless continuity equation:

$$u_r+\frac{u}{r}+w_z=0,$$
(23)
and mass flux condition:

$$r_0^{h(r)}u(r,z)𝑑z=1.$$
(24)
The boundary conditions have been reduced to:

$$\begin{array}{c}u(r,0)=w(r,0)=0\hfill \\ u_z|_{z=h(r)}=0.\hfill \end{array}$$
(25)
In addition to these conditions, boundary conditions in the radial direction also need to be specified. We do not elaborate on them, however, since the in- and outlet conditions arise naturally without the need for prescription when we obtain a simplified system.

The boundary layer equations (22)–(25) form a closed system and can be solved numerically, but pose a difficulty when separated regions exist. Suppose that there is a separation point at $`r=r_s`$ and $`z=0`$ on a flat plate where the skin friction $`u_z`$ vanishes. In its vicinity one finds  that generic solutions of (22) develop singularities of the Goldstein-type $`u\sqrt{r_sr}`$, $`w1/\sqrt{r_sr}`$. On the other hand, experiments  show separation and reversed flow just behind a jump, so it is necessary to overcome this difficulty, which is well-known in the “usual” boundary layers around a body immersed in a high Reynolds number external flow. No such singularities are observed in numerics of the full Navier-Stokes equations in that case, and thus the trouble is thought to be due to truncation of the terms involving higher derivatives in $`r`$, i.e. the terms in (16) of the order $`O(ϵ^2)`$ and higher. An attempt to include those terms leads to intractable equations, so the inverse method  is often used. In this method the feedback from the boundary layer into the external potential flow is taken into account, and the coupled system is iteratively solved to remove the singularity. Without such an external flow present for the circular hydraulic jump, Higuera  has still obtained the velocity and height profiles from the boundary layer equations. His method, called marginal separation, is to force the boundary layer equations through the point of separation by choosing a special non-divergent velocity profile at the point. The physical reasoning for the choice of such a particular profile is rather unclear. Since our aim is also to obtain a simple tractable model, we have chosen a different approach.

### 2.4 Averaged equations

Rather than solving the partial differential equation (22) itself, we shall be content with satisfying only the mass and momentum conservation laws, derived from averaging (22) over the transverse $`z`$-direction. To do this we make an ansatz for the radial velocity profile $`u`$. One might expect that the singularities at separation points do not contribute to the averages and do not cause any harm. Such an expectation is too naive as shown in the next section, since the model still shows singular behavior near the jump if the simplest velocity profile is assumed. Nevertheless, we show in Sec. 2.7 that the model becomes capable of going through the jump smoothly once enough flexibility is introduced in the assumed profile.

We first define the average velocity at $`r`$ by

$$v=\frac{1}{h}_0^hu(r,z)𝑑z.$$
(26)
The total mass flux condition (24) can be written as

$$rhv=1.$$
(27)
Next, for each fixed $`r`$, we integrate the radial momentum equation (22) over $`z`$ from 0 to $`h(r)`$, and use the continuity equation (23) with the surface boundary conditions (25). We obtain the averaged momentum equation

$$\frac{1}{rh}\frac{d}{dr}\left[r_0^hu^2\text{d}z\right]=h^{}\frac{1}{h}u_z|_{z=0}.$$
(28)
Using $`v`$ and

$$G=\frac{1}{h}_0^h\left(\frac{u}{v}\right)^2𝑑z,$$
(29)
we obtain

$$v(Gv)^{}=h^{}\frac{1}{h}u_z|_{z=0}.$$
(30)
Equations (27) and (30) are the total mass and momentum equations.

### 2.5 Similarity profile for $`u`$

The simplest assumption for the radial velocity profile is a self-similar ansatz:

$$u(r,z)/v(r)=f(\eta )$$
(31)
where $`\eta =z/h(r)`$ takes values between 0 (bottom) and 1 (surface). Using (23), the ansatz can be rewritten in the alternate form: $`w(r,z)=\eta h^{}u(r,z)`$. It is also equivalent to the requirement that the local inclination of the streamlines at $`(r,z)`$ be proportional to $`\eta h^{}=zh^{}(r)/h(r)`$. Clearly, such an ansatz is too simple and “rigid” to describe a flow with separation. However, this is the assumption used in the previous literature, and we summarize its consequences. For more details, see .

The conditions (25) and (27) now imply

$$\begin{array}{c}f(0)=0,\hfill \\ f^{}(1)=0,\hfill \\ _0^1f(\eta )𝑑\eta =1.\hfill \end{array}$$
(32)
They are not sufficient to uniquely determine $`f`$, and we choose one that is physically reasonable. Thus, a parabolic profile $`f(\eta )=3\eta 3/2\eta ^2`$ is a simple candidate. Using this choice, $`G=6/5`$ is a constant from (29), and (30) becomes

$$\frac{6}{5}vv^{}=h^{}\frac{3v}{h^2}.$$
(33)
Other choices for $`f`$ lead to the same equation with different numerical coefficients. Since all such equations, corresponding to different choices of $`f(\eta )`$, can be further transformed to

$$vv^{}=h^{}\frac{v}{h^2}$$
(34)
by suitably including numerical coefficients in the characteristic scales (15), the choice of $`f`$ is not important in the study of qualitative behaviour and of parameter dependence.

Using (27), the equation reduces to a single ordinary differential equation for $`v(r)`$:

$$v^{}\left(v\frac{1}{v^2r}\right)=\frac{1}{vr^2}v^3r^2.$$
(35)
This Kurihara-Tani equation was derived and studied in , in its dimensional form, and in . The results can be summarized as follows. To find a solution corresponding to a hydraulic jump, the velocity $`v`$ should be large for small $`r`$, and decrease smoothly as $`r`$ increases. However, the model does not have such a solution. The coefficient of $`v^{}`$ on the left hand side generically vanishes at some $`r`$ where $`v^{}`$ diverges. If (35) is solved in a parametric form on the $`(r,v)`$-plane, all solutions spiral around and into the fixed point $`(r,v)=(1,1)`$, that is a stable focus in the plane. Therefore, one must still connect solutions in the interior and the exterior by means of, e.g., a Rayleigh shock across which mass and momentum flux are conserved. When this is carried out, one finds that the shock occurs very close to $`r=1`$ in the dimensionless coordinates, implying that the radius of the jump in the dimensional coordinates scales roughly as $`r_{}`$ in (15), i.e.:

$$R_j\left(q^5\nu ^3g^1\right)^{1/8}.$$
(36)
This scaling relation (36) was compared to experiments  by changing $`q`$ for several different $`\nu `$. The radius of the jump indeed scaled with the mass flux $`q`$, but the exponent observed in the experiment was about 3/4 rather than 5/8 suggested by (36). To explain the discrepancy, $`R_j`$ was calculated more accurately . It was first proven that there is no solution for $`v(r)`$ to the Kurihara-Tani equation that extends to $`r=\mathrm{}`$. All solutions were found to diverge at some $`r=r_{\text{end}}`$ (constant) like $`h\left\{\mathrm{log}(r_{\text{end}}/r)\right\}^{1/4}`$. By identifying this singularity as the end of the plate where the water runs off, one may always find the solution of (35) diverging at the end of the plate of a given radius $`r=r_{\text{end}}`$. By following the solution to smaller $`r`$, the solution before the jump and the position of the shock are uniquely determined assuming a connection via a Rayleigh shock. The shock location constructed in this way showed a good agreement  with the experiment.

### 2.6 Profile with a shape parameter

An ansatz more flexible than (31) must be used for resolving the flow pattern in the vicinity of the jump. We shall allow the function $`f`$ in (31) to depend also on $`r`$. The simplest modification we can make is to assume $`f=f(\eta ,\lambda (r))`$ so that the velocity profile is characterized by a single “shape parameter” $`\lambda (r)`$. The approach follows the ideas developed by von Karman and Pohlhausen  for the usual boundary layer flow around a body. There, separation of the boundary layer can occur when the pressure gradient, imposed by the external inviscid flow, becomes adverse. In our case, there is no external flow, but there is a pressure gradient, along the bottom $`z=0`$, that is proportional to $`h^{}(r)`$ due to the hydrostatic pressure (21). Thus, the possibility arises that the flow separates on $`z=0`$ near the jump where $`h^{}`$ is large and pressure is increasing in $`r`$, as in the usual boundary layer flow.

As an improvement over the parabolic profile, we approximate the velocity profile by the cubic:

$$u(r,z)/v(r)=a\eta +b\eta ^2+c\eta ^3,$$
(37)
where $`a`$, $`b`$, $`c`$ are now functions of $`r`$. Due to the boundary condition (25) and mass flux condition (27), the coefficients $`a`$, $`b`$, and $`c`$ can be expressed in terms of one parameter $`\lambda `$ as, for example:

$$a=\lambda +3,b=(5\lambda +3)/2,c=4\lambda /3.$$
(38)
The separation condition

$$u_z|_{z=0}=0$$
(39)
is now equivalent to $`a=0`$, or $`\lambda =3`$. The $`u`$-profile is parabolic when $`c=0`$, or $`\lambda =0`$.

Now that we have two unknowns $`h(r)`$ and $`\lambda (r)`$, two equations are necessary. We use the averaged momentum equation (30) as the first equation. Note that $`G`$ is now not a constant, but depends on the shape parameter $`\lambda `$. From (29), we obtain

$$G(\lambda )=\frac{6}{5}\frac{\lambda }{15}+\frac{\lambda ^2}{105}.$$
(40)
Following the Karman-Pohlhausen choice, we choose the second equation to be the momentum equation (22) evaluated at $`z=0`$:

$$h^{}=u_{zz}|_{z=0}.$$
(41)
This connects the pressure gradient on $`z=0`$ with $`\lambda `$. Using (38) and (40), the two equations (30) and (41) can be written as

$$\begin{array}{c}v\frac{d}{dr}\left\{G(\lambda )v\right\}=h^{}\frac{v}{h^2}(\lambda +3)\hfill \\ h^{}=\frac{v}{h^2}(5\lambda +3)\hfill \end{array}$$
(42)
which can be simplified to

$$\begin{array}{c}(G(\lambda )v)^{}=\frac{4\lambda }{h^2}\hfill \\ h^{}=v\frac{5\lambda +3}{h^2}.\hfill \end{array}$$
(43)
Finally, eliminating $`v`$ using (27), we obtain a nonautonomous system of two ordinary differential equations for $`h(r)`$ and $`\lambda (r)`$:

$$\begin{array}{c}h^{}=\frac{5\lambda +3}{rh^3}\hfill \\ \frac{dG}{d\lambda }\lambda ^{}=\frac{4r\lambda }{h}+G(\lambda )\frac{h^4(5\lambda +3)}{rh^4}.\hfill \end{array}$$
(44)
This is the model for the stationary circular hydraulic jump. It does become singular, but only on the lines $`h=0`$ and $`\lambda =7/2`$ which does not cause any trouble in describing a flow with a separated zone ($`\lambda <3`$). We show in the next section that the highly simplified model indeed contains solutions which describe the observed circular hydraulic jumps. A similar approach using momentum and energy conservation was used in , but they did not succeed in finding continuous solutions through the jump.

### 2.7 Numerical solution of the integrated model

The model (44) can be solved as a boundary value problem by specifying two boundary conditions for different values of $`r`$. Thus we impose

$$(r_1,h_1(r_1))\text{and}(r_2,h_2(r_2)),r_1<r_2$$
(45)
where the values are taken from the measured surface height data. There is no fitting parameter once they are chosen, and the function $`h(r)`$ and the shape parameter $`\lambda (r)`$ are determined. In particular, we do not need to specify the shape parameter as a part of the boundary conditions. This is an advantage of the simplified model since one no longer needs to specify the velocity profile at the inlet and/or outlet boundaries, which is not easy to do. In fact, we see that specifying both $`h`$ and $`\lambda `$ at one $`r`$, either inside the jump or outside, and solving (44) as an initial value problem is unstable. The system is extremely sensitive to the initial condition if one integrates (44) in the direction of increasing $`r`$ from a small $`r`$ or in the direction of decreasing $`r`$ from a large $`r`$. Therefore, we choose $`r_1`$ and $`r_2`$ near 1, typically $`r_1`$ around 0.4-0.8 and $`r_2`$ around 1.2-1.6. Then, a straightforward shooting method from either boundary is sufficient to obtain a solution. After this is achieved, the solution is extended to $`r<r_1`$ and to $`r>r_2`$ by integrating (44) backward from $`r_1`$ and forward from $`r_2`$, respectively. Integrations in these directions are stable.

Figure 4(a) shows two solutions of such a boundary value problem. They correspond to the two type I solutions in Fig. 2, reproduced here as dot-dashed curves. From each curve the boundary data are taken at $`\stackrel{~}{r}_1=11.8`$\[mm\] (corresponding to dimensionless value $`r_1=0.42`$) and $`\stackrel{~}{r}_2=30.0`$\[mm\] (to $`r_2=1.07`$). The computed solutions $`h(r)`$ corresponding to the data are shown in solid curves. Each curve shows a gradual decrease for small $`\stackrel{~}{r}`$ as $`\stackrel{~}{r}`$ increases, reaches a minimum at some $`\stackrel{~}{r}15`$\[mm\], and then undergoes a sharp jump at $`\stackrel{~}{r}22`$\- 23\[mm\], and a slow decay after the jump. The location of the jump is about 10% off in each case, and the slope behind the jump is noticeably different. However, the qualitative behavior is well captured by the simple model. Figure 4(b) shows the shape parameter $`\lambda `$. The velocity profile changes suddenly almost simultaneously with the rapid increase of the surface height, and a region where $`\lambda <3`$, corresponding to separation, is observed in each case.666 If the downstream height is further reduced, however, the shape parameter $`\lambda `$ does not reach $`\lambda =3`$, and there is no separated region. Thus, our model predicts that a (weaker) jump without an eddy is possible. The flow near the bottom still decelerates just after the jump. The parameter $`\lambda (r)`$ recovers and appears to converge to $`\lambda =0`$ (the parabolic profile) as $`r`$ becomes large.

The flow structure is more directly shown in Fig. 5, where the $`u`$-velocity profiles are computed from $`\lambda `$ at equidistant locations in $`r`$. Since magnitudes of the velocity vary a lot between small and large $`r`$, the profiles are scaled by the average velocity, so that the profiles of $`u(r,z)/v(r)`$ are shown. The stream function $`\psi `$ is computed from the definition

$$u=\psi _z/r,w=\psi _r/r.$$
(46)
The dimensionless stream function varies from $`\psi =0`$ on $`z=0`$ to $`\psi =1`$ on $`z=h`$. Inside the separated region $`\psi <0`$. The contours at $`\psi =0.1,0,0.1\mathrm{},1`$ are shown in the figure. That is, a region between two neighboring contour curves carries 10% of the mass flux.

The surface velocity $`U`$ predicted from the model is shown in Fig. 6. The parameters are the ones used in Fig. 5. The model again misses the location of the jump by about 20%, so measurements and the curve from the model are offset, but qualitative features are well reproduced. The velocity outside the jump is small and decays like $`U1/r`$, as can be seen from the log-log plot in the inset. This is consistent with an almost constant $`h`$ and a nearly parabolic velocity profile, which we analytically demonstrate in the next section. On the other hand, the surface velocity decreases almost linearly before the jump. This region is harder to explain intuitively, but an analytical approximation is also obtained in the next section. At the jump a rapid, cusp-like drop in the velocity is noticed.

Finally, we discuss the dependence of the solutions on the external height $`h_{ext}`$. Both in experiments and in the model the height inside the jump is little affected by the change in the external boundary condition $`h_2(r_2)`$. The numerical solutions as well as the measured surface profiles in Fig. 4(a,b) apparently overlap in the interior to the jump. Of course, the two solutions must correspond to different trajectories of the model (43) and cannot collapse exactly onto a single curve. However, the closeness of the solution curves in the interior to the jump is the cause of the difficulty of solving the initial value problem starting from a small $`r`$.

If the external height is further increased, a transition from type I to II is observed in the experiment, as illustrated in Fig. 2 and Fig. 3. Unfortunately, no such transition is reproduced in the model when $`h_2`$ is increased. Instead, one finds a computed solution of the model similar to the ones in Fig. 4 even for a much larger $`h_2`$. A physical mechanism to “break” the wave into a type II flow appears to be missing. In fact, a solution with a roller is prohibited by the model (43). The surface velocity on a roller is negative (inward). According to (38), the velocity at the surface is

$$U=v(a+b+c)=v\frac{9\lambda }{3},$$
(47)
where $`v>0`$ is the average velocity. Thus, $`U<0`$ iff $`\lambda >9`$. However, since we start with $`\lambda 0`$ and the line $`\lambda =7/2`$ makes (43) singular, a solution with a roller is not possible. It seems likely that this behavior can be traced back to the assumed pressure distribution (21) which does not provide any pressure gradient along the surface $`z=h`$. In a recent simulation of the circular hydraulic jump by Yokoi et al. pressure buildup just behind the jump is observed and claimed to be crucial in breaking the jump. The non-hydrostatic pressure arises partly due to the surface tension in (19.1), but also due to the truncated viscous terms in (18) and (19). We do not know at present how best to extend our model to include the type II flows.

### 2.8 Asymptotic analysis of the averaged system

In this section we approximate the solutions of (43) analytically using formal perturbation expansions. We obtain explicit expressions for two “outer” regions: the region before the jump and the one after the jump. Moreover, we derive a single ordinary differential equation for the “inner” region near the jump. Analysis in the inner region connects a previous model using a Rayleigh shock with our model.

#### 2.8.1 Outer solution 1 (before the jump)

First, we analyse the region before the jump where thickness of the fluid as well as the radius are small, compared to the exterior region. We denote the typical thickness, in the dimensionless coordinates, as $`\theta `$, and treat it as a formal small parameter. We rescale the variables into $`H`$, $`R`$, and $`V`$ as

$$\begin{array}{c}h=\theta H,\hfill \\ r=\theta ^\alpha R,\hfill \\ v=\theta ^{1\alpha }V,\hfill \end{array}$$
(48)
and require consistent balance of the terms in (44) or, equivalently, (43). The rescaling for $`v`$ in the third equation of (48), is chosen to ensure mass conservation (27) for all $`\theta `$. In terms of the new variables, (43) can be written as

$$\begin{array}{c}\theta ^{2\alpha 1}\frac{d}{dR}\left(G(\lambda )V\right)=\theta ^2\frac{4\lambda }{H^2},\hfill \\ \theta ^{1\alpha }H^{}=\theta ^{\alpha 3}V\frac{5\lambda +3}{H^2}.\hfill \end{array}$$
(49)
From the first equation the only consistent choice is to take $`\alpha =1/2`$. Then, in order to balance the power of $`\theta `$ on both sides of the second equation, we need

$$\lambda =3/5+\theta ^4\lambda _1+\mathrm{}.$$
(50)
The form is also motivated by Fig. 4 in which $`\lambda `$ stays close to the value $`0.6`$ before the jump.

To find $`H(R)`$ and the correction $`\lambda _1`$, substitute (50) into the first equation of (49). To the lowest order in $`\theta `$ we obtain

$$G(0.6)\left(\frac{1}{H^2R}\frac{dH}{dR}+\frac{1}{HR^2}\right)=\frac{12}{5H^2},$$
(51)
where $`G(0.6)=1088/8751.243`$. Solving this equation yields

$$H=\frac{C_1}{R}+\frac{4}{5G(0.6)}R^2,$$
(52)
where $`C_1`$ is an arbitrary integration constant. The functional form agrees with Watson’s self-similar solutions . We also compare the lowest order term of $`\theta `$ in the second equation of (49), and find that

$$\lambda _1=\frac{RH^3}{5}\frac{dH}{dR}.$$

By substituting $`H`$ in (52) we obtain an approximate expression for $`\lambda `$:

$$\lambda =\frac{3}{5}+\theta ^4\left[\frac{H^4}{5}\frac{12}{25G(0.6)}R^2H^3\right].$$
(53)
We test the approximations (52) and (53) in Fig. (4). The dashed curves are the theoretical curves of $`H(R)`$ and $`\lambda (R)`$, shown in the dimensional coordinates. They match the numerical solutions and the measurements well before the jump. Here, the formal parameter $`\theta `$ is taken as unity, and the one free parameter $`C_1`$ was fitted to be $`0.25`$.

#### 2.8.2 Outer solution 2 (after the jump)

Let us now consider the behavior of (43) for large $`r`$. We again introduce a formal small parameter $`\theta `$, but we now rescale $`r=\theta ^1R`$. If we moreover assume that the height is of order $`1`$, i.e., $`h=H`$, then the rescaling of the velocity is necessarily $`v=\theta V`$ due to (24). Using these new variables, Eqs. (43) become:

$$\begin{array}{c}\theta ^2\frac{d}{dR}\left(G(\lambda )V\right)=\frac{4\lambda }{H^2}\hfill \\ \frac{dH}{dR}=V\frac{5\lambda +3}{H^2}.\hfill \end{array}$$
(54)
In order to balance the terms in the first equation we choose

$$\lambda =\theta ^2\lambda _1+\mathrm{}.$$
(55)
This is again consistent with the bottom panel of Fig. 4 where $`\lambda `$ apparently tends to $`0`$, corresponding to the parabolic profile. Then, the terms of order unity in the second equation are

$$\frac{dH}{dR}=\frac{3}{RH^3}$$
(56)
whose solution is

$$H=\left(12\mathrm{log}\frac{R_{\text{end}}}{R}\right)^{1/4}$$
(57)
where $`R_{\text{end}}`$ is an integration constant representing the radius where the height goes to $`0`$. Thus, (43), as well as the simpler Kurihara-Tani model (33), becomes singular when $`r\mathrm{}`$. This seems to be a general property of models based on the boundary layer equations . The absence of regular solutions for the system (22)-(25) when $`r\mathrm{}`$ was proved in . We have attributed this lack of asymptotic solutions to the influence of the finite size of the plate. Indeed, a solution with vanishing height such as (57) reminds one very much of a flow running off the edge of a circular plate.

The height $`H(R)`$, given by equation (57), is a very slowly varying function of $`R`$. There is a long regime $`1RR_{\text{end}}`$ where the height appears almost constant. In this intermediate regime the leading order of (54.1) becomes

$$G(0)\frac{d}{dR}\left(\frac{1}{RH}\right)=\frac{4\lambda _1}{H^2}$$
(58)
where $`G(0)=6/5`$. Therefore,

$$\lambda =\theta ^2\lambda _1=\theta ^2\frac{G(0)H^2}{4}\left(\frac{1}{RH^2}\frac{dH}{dR}+\frac{1}{R^2H}\right)\frac{3}{10r^2}\left(\frac{3}{H^3}H\right).$$
(59)
We conclude that $`\lambda (R)1/R^20`$ which explains the observed approach to the parabolic velocity profile for large $`r`$.

#### 2.8.3 Inner solution near the jump: conservation of momentum

Finally, we analyze the region around the hydraulic jump. Recall that in the Kurihara-Tani theory (33) the jump was obtained by fitting a Rayleigh shock. In this section, we show that our model (44) is a natural generalization of the equation.

To do this we return to (42), and introduce a formal parameter $`\mu `$ in the left-hand side of the second equation.

$$\begin{array}{c}v\frac{d}{dr}\left\{G(\lambda )v\right\}=h^{}\frac{v}{h^2}(\lambda +3)\hfill \\ \mu h^{}=\frac{1}{rh^3}(5\lambda +3).\hfill \end{array}$$
(60)
where $`v=1/(rh)`$. The first equation describes the balance of inertia, hydrostatic pressure, and viscous forces. The value $`\mu =1`$ corresponds to (42).

Setting $`\mu =0`$ gives $`\lambda =0.6`$. Then the first equation becomes the Kurihara-Tani equation (33), except that the coefficient $`6/5=1.2`$ is changed to $`G(0.6)1.243`$ here, since the profile is not parabolic. (As discussed before, the velocity profile is not so important in their model as long as it is self-similar.) Since our model corresponds to $`\mu =1`$, the parameter $`\mu `$ interpolates between the two models, but the correspondence of the two is not obvious because the limit $`\mu =0`$ is a singular limit. We treat $`\mu `$ as a formal small parameter, and carry out a singular perturbation analysis to investigate the connection as well as to obtain an approximation in the jump region.

In Kurihara-Tani model a shock is needed to extend the solution from small to large values of $`r`$. Suppose the shock is situated at $`r=r_0`$. Consider a small region of size $`\mu `$ around $`r=r_0`$, and rescale the coordinate as $`r=r_0+\mu X`$. Then, in the inner coordinate $`X`$, Eq. (60) becomes

$$\begin{array}{c}\frac{1}{r_0h}\frac{d}{dX}\left\{\frac{G(\lambda )}{r_0h}\right\}=\frac{dh}{dX}+O(\mu )\hfill \\ \frac{dh}{dX}=\frac{5\lambda +3}{r_0h^3}+O(\mu ).\hfill \end{array}$$
(61)
We see that $`\lambda =0.6`$ with $`h`$ an arbitrary constant are the only possible fixed points of (61). Thus the solutions must satisfy $`\lambda 0.6`$ for $`X\pm \mathrm{}`$. This correctly matches the external solution before the jump, but not after the jump, where $`\lambda 0`$.777 Note that the singularity of the outer solution after the jump (57)-(59) for $`r0`$ does not allow correct matching for $`X+\mathrm{}`$ when $`\mu 0`$. Nevertheless, our method reproduces the structure of the separation zone quite well.  The first equation can be integrated once, giving the momentum conservation.

$$\frac{G(\lambda )}{r_0^2h}+\frac{h^2}{2}=C_3$$
(62)
with an integration constant $`C_3`$. Now we solve the second equation of (61) for $`\lambda `$, and substitute it into this equation. Using (40) in the form $`G(\lambda )=\frac{1}{105}\left(\lambda \frac{7}{2}\right)^2+\frac{13}{12}`$, we obtain an ordinary differential equation for $`h`$ only:

$$\frac{1}{105}\left(\frac{r_0h^3}{5}\frac{dh}{dX}+\frac{41}{10}\right)^2+\frac{13}{12}+\frac{r_0^2h^3}{2}=C_3r_0^2h.$$
(63)
We look for a solution $`h(X)`$ with $`hh_1`$ as $`X\mathrm{}`$ and $`hh_2`$ as $`X+\mathrm{}`$ where $`h_1`$ and $`h_2`$ are constants. Then, Eq. (63) with the first boundary condition determines the constant $`C_3`$ in terms of $`r_0`$ and $`h_1`$. Eliminating $`C_3`$ we obtain

$$\frac{1}{105}\left[\left(\frac{r_0h^3}{5}\frac{dh}{dX}+\frac{41}{10}\right)^2h_1\left(\frac{41}{10}\right)^2h\right]+\frac{13}{12}(h_1h)\frac{r_0^2}{2}h_1h(h_1^2h^2)=0.$$
(64)
Plugging the second boundary condition into this equation yields a relation between $`h_1`$ and $`h_2`$, given $`r_0`$.

$$h_1h_2^2+h_1^2h_22h_c^3=0.$$
(65)
where

$$h_c=(G(0.6)/r_0^2)^{1/3}$$
(66)
is the critical height for the circular hydraulic jump.888 In dimensional variables, the critical height is $`\stackrel{~}{h}_c=(G(0.6)q^2/g\stackrel{~}{r}_0^2)^{1/3}`$. This is identical to the critical height (4) that appeared in the Rayleigh shock, apart from the numerical factor and the influence of $`\stackrel{~}{r}_0`$ reflecting the radial geometry. The viscosity $`\nu `$ only enters in the coefficient of $`dh/dX`$ in the dimensional version of (64), thus does not affect $`\stackrel{~}{h}_c`$. Solving this equation, we obtain an equation analogous to the shock condition (3):

$$\frac{h_2}{h_1}=\frac{1}{2}\left(1+\sqrt{1+8(h_c/h_1)^3}\right)=\frac{2}{1+\sqrt{1+8(h_c/h_2)^3}}.$$
(67)
It is easy to see that $`h_c`$ is always between $`h_1`$ and $`h_2`$, i.e., $`h_1<h_c<h_2`$ or $`h_2<h_c<h_1`$. The Froude number in this case could naturally be defined as $`F(X)^2=(h_c/h(X))^3`$.999 However, it is not clear whether $`F`$ defined in this way can be a measure of super- and subcriticality since the governing equations are not the shallow water equations and therefore propagation of disturbances do not obey the well-known velocity $`\sqrt{gh}`$.

When $`h_1`$ is close to $`h_c`$, the final height $`h_2`$ is close to $`h_c`$ as well. Then, the Froude number is close to unity for all $`X`$, and the jump is weak, i.e. $`h_ch_1=\delta 1`$. Then, we see from the balance of the terms in (64) that $`h=h_c+\delta Y(\delta x)`$. The leading balance reduces to

$$Y^{}=\gamma (1Y^2)$$

with

$$\gamma =\frac{196875}{1312}\left(\frac{7}{17}\right)^{2/3}r_0^{5/3}83.1r_0^{5/3}.$$
(68)
Thus, in the weak jump limit, the height is given by

$$h(x)=h_c+\delta \mathrm{tanh}(\delta \gamma x).$$
(69)
It is interesting to note that we can connect from $`h_1`$ at $`X=\mathrm{}`$ to $`h_2`$ at $`X=+\mathrm{}`$ if $`h_1<h_2`$, but not if $`h_1>h_2`$, just like in the Rayleigh shock. This requirement comes from the equation (64) self-consistently rather than making a hypothesis on the energy loss like we did in (5). To see this, consider the stability of the fixed points $`h_1`$ and $`h_2`$ with respect to the governing equation (64) for $`h`$.101010Of course, this stability analysis is to study existence of stationary solutions, and not to study the stability of such solutions in the time-dependent theory. Linearizing (64) around the uniform solutions $`h_i`$ (where $`i=1,2`$), we obtain an equation for the perturbation $`\delta h_i`$ in the height:

$$\frac{d}{dX}(\delta h_i)=K_i\delta h_i$$

where

$$K_i=\frac{2625}{41}\frac{r_0\{2h_c^3+h_1(h_1^23h_i^2)\}}{2h_i^3h_1}.$$
(70)
If $`h_1<h_c<h_2`$, then $`K_1>0>K_2`$, showing that the fixed point $`h=h_1`$ is unstable and $`h=h_2`$ stable. A trajectory departing from $`h_1`$ at $`X=\mathrm{}`$ and arriving at $`h_2`$ at $`X=+\mathrm{}`$ is not prohibited, and we can indeed find such a trajectory shown in Fig. 7. In contrast, if $`h_1>h_c>h_2`$, then the stability of the fixed points is reversed, and there is no trajectory going from $`h_1`$ to $`h_2`$.

When $`h_1<h_c<h_2`$ so that such a trajectory exists, the departure from $`h_1`$ is generally rapid, giving an impression of a “sharp corner” at the beginning of the jump, and the arrival at $`h_2`$ is much smoother just as shown in Fig. 7. This is because the magnitude of the stability coefficient $`K_1`$ is large compared to that of $`K_2`$. The feature is most pronounced when $`h_1`$ is small (so, $`h_2`$ is large). It vanishes as $`(h_2h_1)0`$ when $`K_1`$ and $`K_2`$ both tend to zero.

In Fig. 7 we compare solutions of (64) with the two solutions of the full numerical solution of (43) shown in Fig. 4. The jump region is enlarged. Solutions of (64), shown as solid curves, are computed by fitting the values for $`h_1`$ and $`h_2`$, and solving the equation using $`r_0`$ obtained from (65) and (66). We chose an initial condition to be somewhere inside the jump, and integrated (64) forward and backward from it. Since (64) has a translational invariance with respect to $`X`$, the initial condition fixes the location of the jump without affecting the shapes of $`h`$ or $`\lambda `$. The analysis assuming $`\mu 0`$ performs surprisingly well against the numerical solution for $`\mu =1`$. The size of the jump region is now of order $`\mu `$, i.e., unity, and the internal structure is non-trivial. The single ordinary equation (64) is capable of describing the eddy formation in this region.

## 3 Flow down an inclined plane

### 3.1 Introduction to the problem

The properties of waves running down an inclined plane is a subject of great theoretical and practical importance, and has attracted the attention of many researchers. Starting with the pioneering work of Kapitsa & Kapitsa , some of the major contributions to this field are found in  . The physical picture is the following. A fixed flux of fluid is constantly poured onto the inclined plane from above. The fluid forms a stream moving downwards under the action of gravity – an idealized model of a river. If the influx of fluid upstream is suddenly increased, it causes the height upstream to increase, and the extra mass of fluid to propagate downstream. In a river, this may be caused by the melting of snow at regions neighbouring the river’s source, or by sudden rain. A river bore, on the other hand, is introduced at the mouth of the river by a tidal wave, for instance, and moves upstream. In both cases, a solitary wave can be formed, moving at a constant velocity $`c`$ without changing its shape.

We are particularly interested in kink-like solitary wave solutions going from one constant height $`h_1`$ to another $`h_2`$. One can identify such a solution with a heteroclinic orbit, connecting two stationary states . The speed $`c`$ depends on how much the fluid level is increased, i.e., the heights $`h_1`$ and $`h_2`$. Alternatively, we can consider $`c`$ as a parameter, and study the existence of the stationary solution $`h\text{const.}`$ depending on $`c`$. It is rather straightforward to see that two solutions with $`hh_1`$ and $`hh_2`$ exist if $`c`$ is sufficiently large. However, even if $`c`$ is in that regime, it is hard to judge whether there exists a smooth solution connecting the two states. Based on the method of averaging in Sec. 2, we develop a simple model which helps us to derive criteria for their existence and to compute the wave form. The model also enables us to ask whether they appear as “Rayleigh shocks” in the sense that the flow is supercritical in front of the kink structure and subcritical behind it. As we shall elaborate, the distinction between super- and subcritical flows is a concept inherent in inviscid shallow water theory, and is not at all obvious for a viscous flow since now the waves will show dispersion as well as damping. Indeed, we find that the wave velocities corresponding to the largest wave lengths will always propagate both forward and backwards, as in a subcritical flow. Nevertheless, if we focus on wavelengths of the order of the depth of the fluid layer, a clear distinction can be made.

There is another kind of flow in the linear geometry in which a sudden thickening of height is observed. This solution is not only relevant for, e.g., the flow of water exiting from a sluice but is also a direct analog of the circular hydraulic jump. The flow streams rapidly in a region immediately after the sluice, and then abruptly slows down at a certain downstream position. It is stationary (i.e. $`c=0`$) with a constant discharge, and is not obtained as a state connecting two “equilibrium” heights. In fact, the rapid flow before the jump cannot be extended arbitrarily far upstream. We shall show that our models provide physically reasonable solutions in this case, too.

In Secs. 3.2 and 3.3 we write down the complete system for the inclined plane problem, non-dimensionalize it, simplify it using the boundary layer approximation, and average over the thickness in two ways. These steps are in parallel with those in Sec. 2, but we go through them briefly not only for completeness but also since the geometry and the characteristic scales are different. To seek stationary and traveling wave solutions, we write the equations in a coordinate frame moving at a constant speed in Sec. 3.4. Traveling waves are studied in detail in Sec. 3.5, and the stationary jumps in Sec. 3.6.

### 3.2 The governing equations

We consider a viscous, incompressible, two-dimensional flow. The coordinate system is $`\stackrel{~}{x}`$ in the downstream direction parallel to the inclined plane, and $`\stackrel{~}{y}`$ in the perpendicular direction above the plate. Denote the velocities in these directions by $`\stackrel{~}{u}(\stackrel{~}{x},\stackrel{~}{y},\stackrel{~}{t})`$ and $`\stackrel{~}{w}(\stackrel{~}{x},\stackrel{~}{y},\stackrel{~}{t})`$, respectively, the pressure by $`\stackrel{~}{p}(\stackrel{~}{x},\stackrel{~}{y},\stackrel{~}{t})`$, and the height by $`\stackrel{~}{h}(\stackrel{~}{x},\stackrel{~}{t})`$. The governing equations for this problem are the continuity equation

$$\stackrel{~}{u}_{\stackrel{~}{x}}+\stackrel{~}{w}_{\stackrel{~}{y}}=0$$
(71)
and the Navier-Stokes equations

$$\begin{array}{c}\stackrel{~}{u}_{\stackrel{~}{t}}+\stackrel{~}{u}\stackrel{~}{u}_{\stackrel{~}{x}}+\stackrel{~}{w}\stackrel{~}{u}_{\stackrel{~}{y}}=\frac{1}{\rho }\stackrel{~}{p}_{\stackrel{~}{x}}+g\mathrm{sin}\alpha +\nu \left(\stackrel{~}{u}_{\stackrel{~}{x}\stackrel{~}{x}}+\stackrel{~}{u}_{\stackrel{~}{y}\stackrel{~}{y}}\right)\hfill \\ \stackrel{~}{w}_{\stackrel{~}{t}}+\stackrel{~}{u}\stackrel{~}{w}_{\stackrel{~}{x}}+\stackrel{~}{w}\stackrel{~}{w}_{\stackrel{~}{y}}=\frac{1}{\rho }\stackrel{~}{p}_{\stackrel{~}{y}}g\mathrm{cos}\alpha +\nu \left(\stackrel{~}{w}_{\stackrel{~}{x}\stackrel{~}{x}}+\stackrel{~}{w}_{\stackrel{~}{y}\stackrel{~}{y}}\right)\hfill \end{array}$$
(72)
Here, $`\alpha `$ is the angle of the inclined plane (between 0 and $`\pi /2`$) measured downward from the horizontal line, and the subscripts denote the partial derivatives as before. The boundary conditions are identical to those of the radial geometry, i.e., (8)–(10), by reading $`\stackrel{~}{r}`$ as $`\stackrel{~}{x}`$ and $`\stackrel{~}{z}`$ as $`\stackrel{~}{y}`$. The local mass flux is:

$$\stackrel{~}{q}(\stackrel{~}{x},\stackrel{~}{t})=_0^{\stackrel{~}{h}(\stackrel{~}{x},\stackrel{~}{t})}\stackrel{~}{u}𝑑\stackrel{~}{y}.$$

Integrating the continuity equation (71) in $`\stackrel{~}{y}`$ over the thickness and using the boundary conditions, we obtain the flux conservation equation:

$$\stackrel{~}{h}_{\stackrel{~}{t}}+\stackrel{~}{q}_{\stackrel{~}{x}}=0.$$
(73)
The equations above form a complete system apart from the inlet and outlet conditions. They possess a trivial stationary solution (Nusselt solution) with a constant $`\stackrel{~}{h}`$ and the parabolic velocity profile:

$$\stackrel{~}{u}(\stackrel{~}{x},\stackrel{~}{y},\stackrel{~}{t})\frac{g\mathrm{sin}\alpha }{\nu }\left(\eta \frac{\eta ^2}{2}\right),$$
(74)
where $`\eta =\stackrel{~}{y}/\stackrel{~}{h}`$. Given this equilibrium flow, the local flux $`\stackrel{~}{q}`$ is also uniform and steady, and is a function of $`\stackrel{~}{h}`$:

$$\stackrel{~}{q}=_0^{\stackrel{~}{h}}\stackrel{~}{u}𝑑\stackrel{~}{y}=\frac{g\stackrel{~}{h}^3\mathrm{sin}\alpha }{3\nu }.$$
(75)
In a non-equilibrium flow we assume that the inclined plane is infinitely long, and the flow sufficiently far downstream approaches this equilibrium flow. We then treat the flow rate $`\stackrel{~}{q}`$ for $`\stackrel{~}{x}\mathrm{}`$ as the characteristic mass flux $`q_{}`$. The corresponding height $`\stackrel{~}{h}`$ using (75) is used as the length scale $`h_{}`$, and $`v_{}=q_{}/h_{}`$ becomes the characteristic velocity. We non-dimensionalize the governing equations by these scales. The continuity equation is unchanged in form:

$$u_x+w_y=0,$$
(76)
and the Navier-Stokes equations become

$$\begin{array}{c}u_t+uu_x+wu_y=p_x+\frac{3}{R}+\frac{1}{R}(u_{xx}+u_{yy})\hfill \\ w_t+uw_x+ww_y=p_y\frac{3}{R\mathrm{tan}\alpha }+\frac{1}{R}(w_{xx}+w_{yy})\hfill \end{array}$$
(77)
where the pressure is normalized to $`\rho u_{}^2`$, and the Reynolds number is

$$R=\frac{v_{}h_{}}{\nu }=\frac{q_{}}{\nu }=\frac{gh_{}^{}{}_{}{}^{3}\mathrm{sin}\alpha }{3\nu ^2}.$$
(78)
The dimensionless mass flux is $`q(x,t)=hv`$ in terms of the average velocity

$$v(x,t)=\frac{1}{h}_0^hu𝑑y$$
(79)
whereby (75) becomes

$$q=hv=h^3$$
(80)
in an equilibrium flow of height $`h`$.

### 3.3 Boundary layer equations and averaged models

Since the flow on the inclined plane is expected to be predominantly in the $`x`$-direction, the boundary layer approximation should be applicable  as long as separation does not occur. In a similar manner as the radial case, the dominant terms of (77) are:

$$\begin{array}{c}u_t+uu_x+wu_y=p_x+\frac{3}{R}+\frac{1}{R}u_{yy}\hfill \\ 0=p_y\frac{3}{R\mathrm{tan}\alpha }.\hfill \end{array}$$
(81)
The dynamic boundary conditions on $`z=h`$ reduce, as before, to:

$$\begin{array}{c}p|_{y=h}=Wh_{xx}\hfill \\ u_y|_{y=h}=0\hfill \end{array}$$
(82)
with the Weber number in this case being

$$W=\frac{\sigma }{\rho h_{}v_{}^{}{}_{}{}^{2}}=\frac{9\sigma }{\rho gh_{}^{}{}_{}{}^{2}\mathrm{sin}^2\alpha }.$$
(83)
From (81.2) and (82), the pressure is hydrostatic with contribution from the surface tension:

$$p(x,y,t)=\frac{3}{R\mathrm{tan}\alpha }\left(h(x,t)y\right)+Wh_{xx}$$
(84)
so, (81.1) becomes

$$u_t+uu_x+wu_y=\frac{3}{R}\frac{3}{R\mathrm{tan}\alpha }h_x+\frac{1}{R}u_{yy}+Wh_{xxx}.$$
(85)
The mass conservation (73) is non- dimensionalized to

$$h_t+(hv)_x=0.$$
(86)
Now, we make an ansatz for the $`u`$-profile, and average over the thickness in order to obtain two simplified models. First, we use the self-similar velocity profile:

$$u(x,y,t)/v(x,t)=f(\eta )$$
(87)
where $`\eta =y/h(x,t)`$ and the function $`f(\eta )`$ satisfies

$$\begin{array}{c}f(0)=0\hfill \\ f^{}(1)=0\hfill \\ _0^1f(\eta )𝑑\eta =1.\hfill \end{array}$$
(88)
Plug this ansatz into (85), multiply it by $`h`$, and average over $`y`$ to obtain

$$(hv)_t+G(hv^2)_x=\frac{3h}{R}\frac{3}{R\mathrm{tan}\alpha }hh_x\frac{3v}{Rh}+Whh_{xxx}$$
(89)
together with the mass conservation (86). Here,

$$G=\frac{1}{h}_0^h(u/v)^2𝑑y=_0^1f^2(\eta )𝑑\eta $$

is a constant for a given profile in this model. We shall use $`G=6/5`$ for concreteness, corresponding to the parabolic profile $`f=3(\eta \eta ^2/2)`$. Equation (89) is the Cartesian analogue of the Kurihara-Tani equation (33), with time-dependent and surface tension terms.

Next, we assume a variable one-parameter profile for $`u`$. As before, we use a third-order polynomial

$$u(x,y,t)=v(x,t)(a\eta +b\eta ^2+c\eta ^3)$$
(90)
with $`a=\lambda +3`$, $`b=(5\lambda +3)/2`$, and $`c=4\lambda /3`$ chosen to satisfy the conditions (88) for $`f`$. The shape parameter $`\lambda (x,t)`$ is the single variable characterizing the velocity profile. To describe the evolution of $`\lambda (x,t)`$ and $`h(x,t)`$ we choose the same set of equations as in the circular hydraulic jump. The first equation is the mass flux equation (86). In addition, we use the momentum equation (85) multiplied by $`h`$ and averaged in $`y`$, and also (85) evaluated at $`y=0`$:

$$\begin{array}{c}(hv)_t+(hv^2G(\lambda ))_x=\frac{3h}{R}\frac{3}{R\mathrm{tan}\alpha }hh_x\frac{v}{Rh}(\lambda +3)+Whh_{xxx}\hfill \\ 0=\frac{3}{R}\frac{3}{R\mathrm{tan}\alpha }h_x\frac{v}{Rh^2}(5\lambda +3)+Wh_{xxx}\hfill \end{array}$$
(91)
where $`G(\lambda )`$ is given by (40) as before. This system can be cast into the more compact form:

$$\begin{array}{c}(hv)_t+(hv^2G(\lambda ))_x=\frac{4v\lambda }{Rh}\hfill \\ h_x\mathrm{cot}\alpha =1\frac{v}{3h^2}(5\lambda +3)+\frac{WR}{3}h_{xxx}.\hfill \end{array}$$
(92)
In the following we call (89) with (86) the “similarity model”111111 The similarity model is the “Shkadov model” considered in  when $`W0`$. and (92) with (86) the “one-parameter model”. Both models inherit the trivial uniform solution from the complete Navier-Stokes model: $`h=v=q1`$, and $`\lambda 0`$ (parabolic profile) for the one-parameter model.

### 3.4 Stationary solutions in a moving coordinate frame

Here, we are concerned with either stationary solutions or traveling waves whose surface profiles may show abrupt changes. Both types of solutions can be sought as stationary solutions in a moving coordinate system with a suitable constant velocity $`c`$, including the possibility $`c=0`$. Thus, we use the traveling wave coordinate $`\xi =xct`$, and rewrite the models within this frame.

Using the chain rule, the mass conservation (86) used in both models becomes

$$ch_\xi +(hv)_\xi =0$$

which can be integrated to

$$ch+hvQ\text{(const.)}$$
(93)
where $`Q`$ is the mass flux, viewed in the moving frame.121212 Note that the flux $`q(x,t)`$ in the laboratory frame is, in general, not a constant. The discharge at the inlet, e.g., at $`x=\mathrm{}`$ must be varied in time accordingly.  The flow must approach the uniform equilibrium flow $`h=1`$ in the $`\xi \mathrm{}`$ limit. Suppose it also approaches another equilibrium flow $`h=h_2`$ in the $`\xi \mathrm{}`$ limit. Then, using (80), the condition becomes

$$ch_2+h_2^3=Q=c+1.$$
(94)
Of course, $`h_2=1`$ is a solution of this equation. In this case we might still be able to find a non-trivial solution of a pulse-like solitary wave form. Such solutions have previously been studied well , and we do not further seek this type of solutions. For a solution of (94) other than $`h_2=1`$, we need

$$c=h_2^2+h_2+1.$$
(95)
The solution that can be positive is

$$h_2=\frac{1+\sqrt{4c3}}{2}$$

which is positive if and only if $`c>1`$.

When $`c>1`$ two different equilibrium solutions exist, and we hope to find a kink-like solution which connects the two limiting flows. However $`c>1`$ is only the necessary condition for its existence. Sufficiency for the existence depends on the models and the parameters: $`R`$, $`\alpha `$, and $`c`$. In Sec. 3.5 we shall clarify the parameter regime for finding such solutions. It turns out that the velocity profiles in this type of solutions do not deviate much from parabolic even in the one-parameter model. In this sense they correspond to somewhat “mild” jumps in terms of the flow structure.

In Sec. 3.6 we find another family of solutions which approaches $`h=1`$ as $`\xi \mathrm{}`$ when $`c<1`$. These solutions do not start from an equilibrium state at $`\xi =\mathrm{}`$. Instead, they are only valid for $`\xi `$ larger than some value $`\xi _0`$. In the similarity model they are not interesting since they approach $`h=1`$ smoothly. However, within the one-parameter model, an abrupt change is developed in both the surface and velocity profiles, sometimes with separation. We interpret this solution, when $`c=0`$, as the analogue of the circular hydraulic jump in the Cartesian geometry.

The presence of surface tension makes the order of the equations higher and makes it more difficult to compute the solutions even when they exist. We assume that $`W`$ is small and negligible, and set $`W=0`$ in this section. Under this assumption we convert the averaged models into the moving coordinate frame at velocity $`c`$. Equation (89) in the similarity model becomes:

$$c(hv)_\xi +\frac{6}{5}(hv^2)_\xi +\frac{3}{R\mathrm{tan}\alpha }hh_\xi =\frac{3v}{Rh}+\frac{3h}{R}.$$
(96)
Using the condition (93), $`v`$ can be eliminated. We obtain a first order differential equation for $`h`$:

$$\frac{dh}{d\xi }=\frac{15}{R}\frac{(h1)(h^2+h+1c)}{c^2h^26(1c)^2+15h^3/(R\mathrm{tan}\alpha )}.$$
(97)
Similarly, (92) in the one-parameter model is converted to:

$$\begin{array}{c}c(hv)_\xi +(hv^2G(\lambda ))_\xi =\frac{4v\lambda }{Rh}\hfill \\ h_\xi \mathrm{cot}\alpha =1\frac{v}{3h^2}(5\lambda +3)\hfill \end{array}$$
(98)
to be solved with (93). One variable, for instance $`v`$, can be eliminated so that the system becomes two-dimensional for $`h`$ and $`\lambda `$.

In the following sections we treat these averaged models as “dynamical systems”, and view $`\xi `$ as a time-like variable. Fixed points of these systems correspond to the uniform, equilibrium solutions of the original time-dependent equations. Note that stability in terms of the variable $`\xi `$ is not equivalent to temporal stability of the original time-dependent equations.

### 3.5 Traveling wave solutions

Due to the relationship (95) which is a one-to-one map between $`c`$ and $`h_2`$ in the range $`c>1`$, we may treat $`h_2`$ or $`c`$ as the primary parameter interchangeably. Using $`h_2`$ as a parameter corresponds physically to varying the height and discharge upstream and then observing the corresponding change in the wave velocity. The condition $`c>1`$ is equivalent to $`h_2>0`$, and $`h_2>1`$ if $`c>3`$. The two regimes $`h_2>1`$ and $`h_2<1`$ are qualitatively different. For $`h_2>1`$ the discharge at $`\xi \mathrm{}`$ is increased, and a forward-facing front travels downstream. As we shall see in this section, this state exists for small enough $`R`$. In contrast, $`h_2<1`$ corresponds to a backward-facing front which is found to exist for large enough $`R`$ but seems to us very likely unstable. Thus, we concentrate on the case $`h_2>1`$ in the following.131313 If we used the geometric mean of the up- and downstream heights $`\sqrt{\stackrel{~}{h}_1\stackrel{~}{h}_2}`$ as the characteristic length, we would obtain equations whose symmetric appearance makes it easy to study the forward- and backward-facing fronts simultaneously. However, we have chosen to scale by the downstream height $`\stackrel{~}{h_1}`$ in order to treat the traveling waves as well as the stationary jumps.

#### 3.5.1 The similarity model

Since (97) is a first order autonomous ordinary differential equation, the necessary condition for the existence of a heteroclinic orbit starting from $`h_2(>1)`$ and arriving at $`h=1`$ is that the fixed point $`h=1`$ is stable and $`h_2`$ is unstable. By linearization, the fixed point $`h=1`$ is found to be stable if

$$c^26(1c)^2+15/(R\mathrm{tan}\alpha )>0$$
(99)
or,

$$R\mathrm{tan}\alpha <\frac{15}{6(1c)^2c^2}=\frac{15}{5h_2^4+10h_2^3+3h_2^22h_21}f_1(h_2)$$
(100)
where the denominator is positive for $`c>3`$. Similarly, $`h_2`$ is found to be unstable if

$$R\mathrm{tan}\alpha <\frac{15h_2}{h_2^42h_2^3+3h_2^2+10h_2+5}f_2(h_2).$$
(101)
The denominator of $`f_2`$ vanishes only at $`h_2=h_2^{\text{max}}2.13`$ for the region $`h_2>1`$. If $`h_2>h_2^{\text{max}}`$, then $`f_2<0`$ and (101) cannot be satisfied. We discard this region of $`h_2`$. For $`1<h_2<h_2^{\text{max}}`$ one finds that $`f_2(h_2)>1>f_1(h_2)`$. Thus, the necessary condition for the existence is simply (100). Once the necessary condition is fulfilled, sufficiency is guaranteed. To see this, we only need to ensure that the denominator on the right hand side of (97) does not vanish in the region $`1<h<h_2`$. Suppose it vanished at $`h_s`$, then we would have

$$c^2h_s^26(1c)^2+15h_s^3/(R\mathrm{tan}\alpha )=0.$$
(102)
Comparison with (99) gives us

$$c^2(1h_s^2)+15(1h_s^3)/(R\mathrm{tan}\alpha )>0.$$

It is clear that $`h_s>1`$ is impossible. Thus, $`h_s<1`$, and there is no vanishing denominator in $`1<h<h_2`$. In Fig. 8(a) we show computed solutions of (97) for three different Reynolds numbers. The parameters $`\alpha `$ and $`h_2`$ are fixed, such that (100) becomes $`R<6.95`$. Within this range, a larger $`R`$ makes the propagating front sharper.

#### 3.5.2 The one-parameter model

We can eliminate $`v`$ from (93) and (98), and think of trajectories on the phase portrait for $`(h,\lambda )`$. We look for a heteroclinic orbit starting from a fixed point $`(h_2,0)`$ and arriving at $`(1,0)`$ as $`\xi \mathrm{}`$. It is necessary for its existence that the point $`(h_2,0)`$ has at least one unstable direction and $`(1,0)`$ has at least one stable direction. Linearizing around the equilibrium point as $`h=h_e+\delta h`$ and $`\lambda =0+\delta \lambda `$, where $`h_e=1`$ or $`h_2`$, we obtain:

$$\left(\begin{array}{c}\delta h_\xi \\ \delta \lambda _\xi \end{array}\right)=J\left(\begin{array}{c}\delta h\\ \delta \lambda \end{array}\right).$$

It is straightforward to calculate the $`2\times 2`$ Jacobian matrix $`J`$, and show that

$$detJ=\frac{60(c3h_e^2)\mathrm{tan}\alpha }{Rh_e^7}.$$
(103)
For the point $`(h_2,0)`$ we have $`c3h_e^2=1+h_22h_2^2<0`$ when $`h_2>1`$. This means that $`detJ<0`$ for $`h_2>1`$, and the fixed point is always a saddle, having exactly one unstable direction.

For the point $`(1,0)`$ we have $`detJ>0`$ since $`c3h_e^2=h_2^2+h_22>0`$ when $`h_2>1`$. Thus, we must also compute the trace of $`J`$ for $`h_e=1`$ which can be shown to be

$$\text{tr}J=\frac{60}{R}+(3361c+25c^2)\mathrm{tan}\alpha .$$

For the stability of $`(1,0)`$ we need tr$`J<0`$. Since $`3361c+25c^2>0`$ for $`c>3`$, this condition becomes

$$R\mathrm{tan}\alpha <\frac{60}{3361c+25c^2}=\frac{60}{311h_2+14h_2^2+50h_2^3+25h_2^4}f_s(h_2).$$
(104)
When this is satisfied, the fixed point is locally attracting, and a trajectory may reach it from any direction. Indeed, we find numerically that the condition (104) also seems to be sufficient. For any $`R`$ and $`\alpha `$ we have tried in the range (104), a heteroclinic solution was found. Computed solutions for three different values of $`R`$ are shown in Fig. 8(b) and (c). The parameters $`\alpha `$ and $`h_2`$ are identical to the ones used for the similarity model in Fig. 8(a). The condition (104) yields $`R<5.59`$. The height profiles in (b) are essentially identical to the ones in (a). This is because the shape parameter $`\lambda `$ shown in (c) does not deviate much from $`\lambda =0`$, the parabolic profile.

In Fig. 8(b) and (c), the solution is oscillatory around $`h=h_1`$ and $`\lambda =0`$ for $`R=5.5`$. This is a feature seen when $`R`$ becomes close to the critical value given by (104). It happens when the type of the fixed point $`(1,0)`$ changes from a stable node to a stable focus. The point is a focus when $`detJ>(\text{tr}J)^2/4`$, which is equivalent to $`f_+(h_2)<R\mathrm{tan}\alpha <f_{}(h_2)`$ where

$$f_\pm (h_2)=\frac{60}{79h_2+16h_2^2+50h_2^3+25h_2^4\pm 2\sqrt{5D}}$$
(105)
and

$$D=2+3h_29h_2^219h_2^3+3h_2^4+15h_2^5+5h_2^6.$$

It can be seen that $`f_+(h_2)<f_s(h_2)<f_{}(h_2)`$ for $`h_2>1`$. Therefore, a heteroclinic solution can be found and exhibits oscillations in a small region $`f_+(h_2)<R\mathrm{tan}\alpha <f_s(h_2)`$. In Fig. 8(b) and (c) this condition corresponds to $`4.81<R<5.59`$, so only the solution for $`R=5.5`$ shows oscillations.

### 3.6 Stationary jumps

If $`c<1`$, the two averaged models have only one fixed point $`h=1`$. Therefore, one might imagine that it is too limited to show any jump-like structures. Nevertheless, we look for trajectories that approach to the fixed point as $`\xi \mathrm{}`$. Even though $`c=0`$ is the physically most interesting case, we treat the general case $`0c<1`$. Since there is no $`h_2`$, we use $`c`$ as the prime parameter in this section.

#### 3.6.1 The similarity model

The sole fixed point $`h=1`$ must be stable to be the limiting point of a trajectory as $`\xi \mathrm{}`$. For $`0c<1`$, the condition is similar to (99) but with reversed inequality

$$c^26(1c)^2+15/(R\mathrm{tan}\alpha )<0.$$
(106)
The singular height $`h_s`$ of the governing equation is still given by (102), and, using a similar argument as before, it is easy to see that $`0h_s<1`$ is impossible when $`c<1`$. Thus, there is a trajectory which approaches $`h=1`$ from below if (106) holds. When $`1>c>(6\sqrt{6})/50.71`$, $`c^26(1c)^2>0`$ and (106) cannot be satisfied. When $`c<(6\sqrt{6})/5`$, the condition is equivalent to

$$R\mathrm{tan}\alpha >\frac{15}{6(1c)^2c^2},$$
(107)
which is satisfied in a range of $`R\mathrm{tan}\alpha `$ since the denominator of the right hand side is positive.

Computed solutions for $`R=50`$, 70, and 100 are shown in Fig. 9 as dashed curves using $`\alpha =3`$\[deg\] and $`c=0`$. The condition (106) becomes $`R>47.7`$, and is satisfied for all three. Each solution simply approaches $`h=1`$ smoothly, clearly reflecting the first order nature of the model (97). As $`\xi `$ decreases, the height vanishes at a finite $`\xi `$ and an inlet must be placed before this happens. If $`h`$ is very small, (97) simplifies to $`dh/d\xi =5/\{2R(1c)\}`$. The solution is

$$h(\xi )=\frac{2.5}{R(1c)}(\xi \xi _0)$$
(108)
for some $`\xi =\xi _0`$ where $`h=0`$. There is no abrupt change in the solutions that resembles a stationary shock structure. If we use $`R`$ smaller than the critical value, then there is no solution converging to $`h=1`$. Therefore, we view the similarity model as inadequate for describing stationary jumps.

#### 3.6.2 The one-parameter model

The sole fixed point of this model when $`c<1`$ is $`(h,\lambda )=(1,0)`$. The Jacobian and its determinant is still given by (103), but now $`c3h_e^2=c3<0`$ and, thus, $`detJ<0`$. Therefore, the fixed point is always a saddle in this range of $`c`$, and there is one direction convergent to the fixed point as $`\xi \mathrm{}`$. It is easy to compute the corresponding trajectory by integrating backward in $`\xi `$ from the vicinity of the fixed point. This solution seems to exist for all values of $`R`$, $`\alpha `$ and $`c<1`$. We are interested in solutions which approach $`h=1`$ from below, and tend to $`h=0`$ at some $`\xi =\xi _0`$ as $`\xi `$ decreases. (To be physical, an inlet condition must be specified at some $`\xi >\xi _0`$.) We can analyze the solutions asymptotically near $`\xi _0`$ by assuming that $`hA(\xi \xi _0)`$ as $`\xi \xi _0+0`$. Then, using (93) and $`Q=1c`$ in (94), we obtain $`v(1c)/\{A(\xi \xi _0)\}`$. Substituting these into (98.2) yields

$$\lambda 0.6+\frac{3A^3}{5(1c)}(1A\mathrm{cot}\alpha )(\xi \xi _0)^3.$$

Finally, comparing coefficients of the dominant terms in (98.1) determines $`A`$ as

$$A=\frac{12}{5RG(0.6)(1c)}\frac{1.93}{R(1c)}.$$

We observe two qualitatively different types depending on the parameter values. If $`\lambda `$ increases at the point $`\xi =\xi _0`$, then the solution reaches the parabolic profile $`\lambda =0`$ monotonically. This occurs when $`R`$ is large, and three computed solutions are shown in Fig. 9 as solid curves. The height profile is qualitatively identical to the ones from the similarity model shown in dashed curves. They do not show any jump structure.

On the other hand, if $`\lambda `$ decreases at $`\xi _0`$, then the trajectory makes an excursion to smaller $`\lambda `$, sometimes into the separation zone $`\lambda <3`$, before recovering toward $`\lambda =0`$. The condition to obtain the second type is $`A\mathrm{cot}\alpha >1`$, or,

$$R\mathrm{tan}\alpha <\frac{12}{5G(0.6)}(1c)1.94(1c)$$
(109)
with $`G(\lambda )`$ given by (40). Two solutions satisfying this condition are shown in Fig. 10(a) and (b). Both the height profile and the shape parameter vary in a similar manner to the one we obtained in the circular hydraulic jump. The phase portrait in (c) demonstrates how rapid and large the excursion can become for small $`R`$. This type of solution could be realized, for instance, as a stationary flow ($`c=0`$) exiting a sluice gate placed at some $`\xi >\xi _0`$.141414 A full-scale channel flow such as a river certainly requires a turbulence modelling, but we have been able to construct a miniature experimental model in which the flow remains laminar. However, our preliminary observation is that a pair of edge waves are created from the ends of the gate, which makes the flow three-dimensional.

## 4 Linear stability of equilibrium states

It is quite difficult to carry out linear stability analysis around the stationary solutions and traveling wave solutions found so far. They have non-uniform profiles obtained only numerically and some of the solutions have singular points beyond which they cannot be continued. Moreover, the inlet boundary condition can strongly affect the stability properties of the solutions. We shall therefore focus on the linear geometry, and only study stability of the equilibrium flow $`h`$ const. The results are, however, expected to be applicable to the equilibrium flow sufficiently far downstream of the jump in the stationary solutions and to flows sufficiently up- and downstream of the moving front in case of the traveling wave solutions. Since the dispersion relation scales with the chosen characteristic length, as described in Sec. 4.4, we only need to consider the flow $`h1`$. Both the similarity model (89) and the one-parameter model (92) are considered, including the surface tension term which is expected to be relevant  for stability. One of our aims is, of course, to judge when infinitesimal disturbances grow and whey they decay, but their propagation velocities are also of our great interest. By comparing the velocities to a reference velocity, which is zero for the stationary jump and $`c(>3)`$ for the traveling wave, we are able to classify different parts of the solutions as either super- or subcritical.

### 4.1 Dispersion relations

The first step is to linearize the models around the fixed point $`h=v=1`$ and, for the one-parameter model, $`\lambda =0`$. We assume infinitesimal disturbances $`\delta h`$, $`\delta v`$, and $`\delta \lambda `$, and decompose them into Fourier modes:

$$\delta h,\delta v,\delta \lambda e^{i(kx\omega t)}.$$
(110)
Plugging them into the linearized equations for the similarity model (86) and (89), we obtain:

$$\omega ^2+\omega \left(\frac{3i}{R}\frac{12}{5}k\right)+\left(\frac{9i}{R}k+\frac{6}{5}k^2\frac{3}{R\mathrm{tan}\alpha }k^2Wk^4\right)=0.$$
(111)
Solving the equation the dispersion relation is found to be

$$\omega _\pm =\frac{3i}{2R}+\frac{6}{5}k\pm \sqrt{D_0}$$
(112)
where the discriminant is

$$D_0=\frac{9}{4R^2}+\frac{27i}{5R}k+3k^2\left(\frac{2}{25}+\frac{1}{R\mathrm{tan}\alpha }\right)+Wk^4.$$
(113)
Similarly, from the one-parameter model (86) and (92), we obtain the dispersion relation:

$$\omega _\pm =\frac{6i}{5R}+\frac{61}{50}k\pm \frac{3}{5}\sqrt{D_1}$$
(114)
where

$$D_1=\frac{4}{R^2}+\frac{178i}{15R}k+\left(\frac{421}{900}+\frac{20}{3R\mathrm{tan}\alpha }\right)k^2+\frac{i}{9\mathrm{tan}\alpha }k^3+\frac{20W}{9}k^4+\frac{iRW}{27}k^5.$$
(115)
Note that this model also has only two dispersion relations, $`\omega _+(k)`$ and $`\omega _{}(k)`$ because the second equation of (92) does not include time-derivatives.

### 4.2 Long wave limit

We first study the long wave limit $`k0`$ by taking only the lowest order terms in $`k`$. For the similarity model, the dispersion relation (112) becomes

$$\begin{array}{c}\omega _+=3k+ik^2(R\mathrm{cot}\alpha )+O(k^3)\hfill \\ \omega _{}=\frac{3i}{R}\frac{3}{5}kik^2(R\mathrm{cot}\alpha )+O(k^3).\hfill \end{array}$$
(116)
As $`k0`$, the group velocities $`d\omega _+/dk3`$ and $`d\omega _{}/dk(3/5)`$. Therefore, waves corresponding to $`\omega _{}`$ propagate upstream, and the flow is subcritical irrespective of $`R`$. By studying the dominant imaginary components of $`\omega _\pm `$, we also find that the reverse propagating branch $`\omega _{}`$ is always stable, i.e. the disturbances decay, for small enough $`k`$ whereas the forward propagating branch $`\omega _+`$ is stable only for small enough Reynolds number satisfying

$$R\mathrm{tan}\alpha <1.$$
(117)
The limiting dispersion is identical in the one-parameter model apart from numerical coefficients. For small $`k`$, (114) becomes

$$\begin{array}{c}\omega _+=3k+ik^2(\frac{5}{4}R\mathrm{cot}\alpha )+O(k^3)\hfill \\ \omega _{}=\frac{12i}{5R}\frac{14}{25}kik^2(\frac{5}{4}R\mathrm{cot}\alpha )+O(k^3).\hfill \end{array}$$
(118)
Thus, the flow is always subcritical since the long waves in the $`\omega _{}`$ branch propagate upstream with velocity $`14/25`$. Again, this branch is stable for any $`R`$ while the $`\omega _+`$ branch is stable only for small Reynolds numbers:

$$R\mathrm{tan}\alpha <4/5.$$
(119)
### 4.3 Intermediate range of $`k`$

It is quite unexpected that the flow is subcritical for any $`R`$. One would intuitively expect that disturbances cannot propagate upstream for sufficiently rapid flows. An explanation can be made by a more careful study of the dispersion relations (112) and (114), or, in particular, the discriminants $`D_0`$ and $`D_1`$.

We first consider the similarity model. If the $`O(k^2)`$ term dominates in $`D_0`$, then the corresponding group velocities become

$$c_\pm =\frac{d\omega _\pm }{dk}\frac{6}{5}\pm \sqrt{\frac{6}{25}+\frac{3}{R\mathrm{tan}\alpha }}.$$
(120)
Both $`c_+`$ and $`c_{}`$ become positive for

$$R\mathrm{tan}\alpha >5/2.$$
(121)
We attempt to estimate such a range of $`k`$. For brevity we assume $`R\mathrm{tan}\alpha 25/2`$ so that the coefficient of $`k^2`$ in $`D_0`$ can be approximated by $`3/(R\mathrm{tan}\alpha )`$. If the magnitude of the $`O(k^2)`$ dominates in $`D_0`$, then we must have

$$\frac{3k^2}{R\mathrm{tan}\alpha }\frac{9}{4R^2},\frac{27k}{5R},Wk^4,$$

that is,

$$\mathrm{max}[\sqrt{\frac{3\mathrm{tan}\alpha }{4R}},\frac{9\mathrm{tan}\alpha }{5}]k\sqrt{\frac{3}{RW\mathrm{tan}\alpha }}.$$
(122)
Using $`R=30`$, $`\alpha =5`$\[deg\], and $`W=0.01`$, for instance, the condition (121) and (122) gives a window $`0.16k10.7`$ in which we can hope that the $`O(k^2)`$ term dominates.

Rather than attempting a more accurate estimate of the zone, we demonstrate that such an interval can be in fact quite long, by plotting the real part of $`\omega _\pm (k)`$ for (112) in Fig. 11(a). Three different values of $`R`$ are used while $`\alpha `$ and $`W`$ are fixed. The $`\omega _+`$ branch, shown as dashed curves, has a positive slope for any $`k`$. Both phase and group velocities of this branch are positive. On the other hand, the $`\omega _{}`$ branch, shown as solid curves, qualitatively changes with $`R`$. For $`R=25`$ its slope appears to be negative for all $`k`$, indicating a subcritical flow. However, for a larger $`R`$ there is an interval of $`k`$ in which the slope becomes positive. In the limit $`k0`$, the branch still has a negative slope in accordance with the analysis of the long wave limit in the previous section. However, the subcritical region near $`k=0`$ can be very small. One sees in Fig. 11(a) that the curve has a positive slope already when $`k>0.05`$ and $`R=35`$. The slope continues to be positive until $`k=2`$, corresponding to a wavelength of half the thickness of the equilibrium flow. Since the system length is finite in practice, the subcritical flow in the $`k0`$ limit cannot be achieved, and the flow becomes essentially supercritical for all the wave numbers observed. This defines the super- and subcritical flows within our viscous model, and confirms the intuitive picture of having a supercritical flow when the flow is sufficiently rapid.

The situation is qualitatively identical in the one-parameter model. We obtain

$$R\mathrm{tan}\alpha >20/11$$
(123)
and

$$\mathrm{max}[\sqrt{\frac{3\mathrm{tan}\alpha }{5R}},\frac{50\mathrm{tan}\alpha }{89}]k\mathrm{min}[\frac{60\mathrm{tan}\alpha }{R},\sqrt{\frac{3}{RW\mathrm{tan}\alpha }},\left(\frac{180}{R^2W\mathrm{tan}\alpha }\right)]$$
(124)
as the corresponding equations to (121) and (122), respectively. Again using $`R=30`$, $`\alpha =5`$\[deg\], and $`W=0.01`$, the interval becomes $`0.05k0.18`$. The upper limit comes from the $`O(k^3)`$ term in $`D_1`$, and is estimated to be rather small since we have only compared the magnitudes. In fact, when we plot the real part of the dispersion relation (114) in Fig. 11(b), we find that the $`\omega _{}`$ branch has a positive group velocity for a much longer range of $`k`$. The supercritical flow near the $`k=0`$ limit is very small once again if $`R`$ becomes as large as $`R=25`$.

### 4.4 Super- and subcriticality for moving fronts

The intermediate-$`k`$ behavior enables us to decide whether a given equilibrium flow is “inherently” super- or subcritical. This distinction is made based on wave velocities with respect to the laboratory frame. A more classical distinction of the two types arises in the context of the shock theory, as reviewed in Sec. 2.1. In this case velocities are measured with respect to a moving front; we call the flow “supercritical” if the group velocity of all the waves is less than the front velocity $`c`$, and “subcritical” if there is a wave component whose group velocity is larger than $`c`$. Here, we briefly note that the averaged equations can describe this traditional classification, too.

Take a moving front such as the one shown in Fig. 8. We concentrate on the long wave limit $`k0`$. For $`\xi \mathrm{}`$ the flow approaches an equilibrium flow with $`h=1`$. Linear waves propagate forward and backward with the group velocities $`d\omega _+/dk=3`$ and $`d\omega _{}/dk=3/5`$ according to the dispersion relation for the similarity model (116). This is a subcritical situation in the laboratory frame, but, since the front velocity is $`c=1+h_2+h_2^2>3`$, both these waves propagate into the front. Therefore, the flow is supercritical with respect to the front.

To derive the dispersion relation of the equilibrium flow with height $`h_2`$ for $`\xi \mathrm{}`$, consider rescaling the height by $`h_2`$. That is, we use this height as the characteristic length so that a wave number $`k`$ must be multiplied by $`h_2`$. Since the flow rate is $`q_2=h_2^3`$ from (80), the velocity has to be scaled by $`q_2/h_2=h_2^2`$. Thus, the group velocities for this flow in the laboratory frame are $`d\omega _+/dk=3h_2^2`$ and $`d\omega _{}/dk=(3/5)h_2^2`$. It is easy to show that $`3h_2^2>c=1+h_2+h_2^2`$ for $`h_2>1`$. Thus, one wave component propagates into the front while the other moves away from it so that the flow behind the front is subcritical.

Therefore, the moving front has a supercritical flow on the shallower side and a subcritical flow on the deeper side, and can be regarded as a classical shock. Using the one-parameter model instead of the similarity model is qualitatively identical.

### 4.5 Short wave limit

We now come back to the stationary equilibrium flow, and study the dispersion relation in the short wave range. Since the derivation of the averaged equations relies on the assumption of predominantly horizontal flow, it is not our aim to accurately resolve wave components when $`k`$ is large. We only hope that the short waves decay so that they do not interfere with meaningful dynamics when we simulate the time-dependent model. Unfortunately, the one-parameter model performs poorly in this respect compared to the similarity model.

The dispersion relation of the similarity model (112) can be approximated in the large $`k`$ limit as

$$\begin{array}{c}\text{Re }\omega _\pm =\pm \sqrt{W}k^2+O(k)\hfill \\ \text{Im }\omega _\pm =\frac{3}{2R}+O(1/k).\hfill \end{array}$$
(125)
Thus, short waves in (89) are damped out if $`W>0`$.

If we neglect the surface tension and set $`W=0`$, the dispersion relation for large $`k`$ is

$$\omega _\pm =c_\pm k\frac{3i}{2R}\frac{c_\pm 3}{c_\pm 6/5}+O(k^1)$$
(126)
where $`c_\pm `$ is the velocity of the corresponding wave given by

$$c_\pm =\frac{6}{5}\pm \sqrt{\frac{6}{25}+\frac{3}{R\mathrm{tan}\alpha }}.$$
(127)
Since $`c_{}<6/5`$ from (127), the branch $`\omega _{}`$ is always stable, as can be seen from (126). On the other hand, since $`c_+>6/5`$, the condition for the stability of the branch $`\omega _+`$ is $`c_+<3`$, which is equivalent to

$$R\mathrm{tan}\alpha <1.$$
(128)
For a large $`R`$ the equilibrium state is no longer stable, but this is reasonable in the absence of surface tension.

Now, we turn into the dispersion relation of the one-parameter model (114). For large $`k`$, it behaves as

$$\omega _\pm \pm k^{5/2}\sqrt{iW/75}\text{if }W>0\text{,}$$
(129)
and as

$$\omega _\pm \pm k^{3/2}\sqrt{i/(25\mathrm{tan}\alpha )}\text{if }W=0\text{.}$$
(130)
In either case one of the branches has an unstable component as $`k\mathrm{}`$, irrespective of $`R`$ or $`\alpha `$. We have been unable to find a natural modification to the one-parameter model which prevents this unphysical behavior. Its cause may well be that the evolution of short waves is not well represented by the boundary layer approximation we started with. In fact, in the boundary layer equations (81) the higher order derivatives of $`x`$ that are thought to be crucial for stability of the high-$`k`$ modes are neglected. In this view the similarity model (89) provides surprisingly reasonable behaviour for large $`k`$, even starting from (81).

## 5 Conclusions

In this article we have presented a simple but fairly quantitative method of reducing flows with strongly deformed free surfaces to a manageable system of equations. By assuming a “flexible” velocity profile whose shape parameter is another dependent variable, flows with an internal eddy can be described. In the radial geometry our results compare well with experiments and we have obtained analytic expressions for the circular hydraulic jump.

We have also studied the flow down an inclined plane. The reduced equations possess not only the traveling wave solutions (heteroclinic orbits) studied previously but also stationary jump solutions. We have found that the stationary solutions show a stronger change in the velocity profile than the traveling waves.

Finally, we have classified different parts of the flows into super- and subcritical by studying the dispersion relation around the equilibrium flow. This classification is standard for inviscid shallow water flow and in shock theory, but is is not obvious in the context of viscous flow. Indeed, for sufficiently long waves the averaged equations show that supercritical flow is not possible. However, waves with intermediate lengths can make the flow essentially supercritical.

The only but serious defect of our reduced model which we have been unable to overcome is its short wavelength behavior. As it stands now, some artificial dissipation term to stabilize the short waves is necessary before time-dependent simulations are attempted. To our dismay a more natural treatment of this problem has so far eluded us.

## Acknowledgements

The core part of this work was carried out while authors were at the Center for Chaos & Turbulence Studies (CATS) at the Niels Bohr Institute to which SW and VP are grateful for hospitality and an inspiring environment. SW thanks the Institute for Mathematics & its Applications (IMA) of the University of Minnesota for providing him with a place and atmosphere to continue the work. Research supported in part under Grant-in-Aid for Scientific Research of JSPS. VP acknowledges the hospitality of University of Chicago and support through the NSF grant No. DMR 9415604 and MRSEC, NSF Grant No. DMR 9808595.