A Nonlinear Formulation of Radiation Stress and Applications to Cnoidal Shoaling
aa r X i v : . [ phy s i c s . f l u - dyn ] F e b A Nonlinear Formulation of Radiation Stress andApplications to Cnoidal Shoaling
Martin O. Paulsen · Henrik Kalisch g Abstract
In this article we provide formulations of energy flux and radiation stress consistent with thescaling regime of the Korteweg-de Vries (KdV) equation. These quantities can be used to describethe shoaling of cnoidal waves approaching a gently sloping beach. The transformation of thesewaves along the slope can be described using the shoaling equations, a set of three nonlinearequations in three unknowns: the wave height H , the set-down ¯ η and the elliptic parameter m .We define a numerical algorithm for the efficient solution of the shoaling equations, and we verifyour shoaling formulation by comparing with experimental data from two sets of experiments aswell as shoaling curves obtained in previous works. Keywords:
Surface waves · KdV equation · Conservation laws · Radiation stress · Set-down
In the present work, the development of surface waves across a gently sloping bottom is in view. Ourmain goal is the prediction of the wave height and the set-down of a periodic wave as it enters anarea of shallower depth. This problem is known as the shoaling problem and has a long history, withcontributions from a large number of authors going back to Green and Boussinesq.In order to reduce the problem to the essential factors, we make a number of simplifying as-sumptions. First, the bottom slope is small enough to allow the waves to continuously adjust to thechanging depth without altering their basic shape or breaking up. Second, reflections are assumed tobe negligible so that the energy flux of a wave is conserved as it transforms on the slope. Third, theperiod T of a wave is assumed to be constant as it progresses. This last point is sometimes termedthe conservation of waves as the number of waves in a group will remain constant during the shoalingprocess.Rather than following a wave in space and time as it propagates over the slope, we estimate waveproperties as functions of the depth using conservation of period T and energy flux, and a prescribedchange in momentum due to the effect of the radiation stress. These conservation equations are dueto the assumptions stated above, and they form the basis for the description of the shoaling processas first envisioned by Rayleigh [34], and subsequently used by a number of authors. The crucial factorwhich guarantees that the shoaling assumptions are valid is that the relative change in water depthover a wavelength is smaller than the wave steepness. This point is explained in more detail in [32, 43].In the context of linear wave theory, the process outlined above is the classical approach to waveshoaling and can be found in textbooks on coastal engineering such as [12, 40]. This approach dependson the assumption that the waves can be described by a sinusoidal function of the form H cos( kx − ωt ),where H is the wave height, k is the wavenumber, ω is the circular frequency, t is the time and x is a1patial coordinate along which the bottom is sloping. Recall that the wavelength is λ = πk and thewave period is T = πω , and these are related by the dispersion relation ω = gk tanh ( kh ) at localdepth h . In the linear shoaling theory, the conservation of energy flux is central. The energy flux isgenerally formulated as the average over one period of the energy density E times the group speed C g . The energy density is given in terms of the fluid density ρ , the gravitational acceleration g andthe wave height H as E = ρgH , and the group speed is given by C g = dωdk .If the bathymetry is sloping only gently upwards, then it may be assumed that the waves adjustadiabatically to the changing conditions, and that reflections and distortions are negligible. Thewave-height transformation of a shoaling wave is obtained by imposing conservation of the energyflux EC g across the shoaling region. If a wave measurement H , T at some point offshore is given,then the wave height at a point closer to shore may be computed from the conservation of energy fluxgiven in the form HC g = H C g , (1)where the group speed C g at local depth h can be found from the conservation of the wave period T . Similar considerations using the linear definition of the radiation stress yield a formula for theset-down [28].The authors of [44] used cnoidal functions in connection with the linear formulation of the energyflux defined above to obtain a hybrid theory of shoaling. Unfortunately, this approach led to adiscontinuity in wave height at the matching point between the linear shoaling equation (1) and thenonlinear theory based on the cnoidal functions. The problem was remedied to some degree in [45] byrequiring continuity in wave height effected through the use of conversion tables. This approach ledto good agreement with the experimental data, but was cumbersome to implement, and as alreadynoted in [45], the continuity of the wave height at the matching point in the shoaling curve led to adiscontinuity in energy flux.The problem can be resolved by defining a shoaling equation based on the full water-wave problemusing the streamfunction method [12] for example. Such an approach has been documented in [46].While accurate, the streamfunction approach depends on discretizing the steady Euler equationsleading to large systems of equations which need to be solved numerically. If the approach is tobe used in connection with ocean-wave statistics, the approach based on the steady Euler equationsmight be too computationally demanding. Moreover, the agreement of the combined linear-cnoidalapproach of [45] already gave good agreement with experiments while being much less expensive.In present work we will show that the discrepancy between conservation of wave height andconservation of energy flux can be resolved in the context of the KdV theory by using a definitionof the energy flux consistent with the KdV equation [4, 15], and without using the fully nonlinearapproach of [46]. Recall that the KdV equation is given in the form η t + c η x + 32 c h ηη x + c h η xxx = 0 , (2)where η ( x, t ) is the deflection of the free surface from rest at a point x and a time t , g is the gravitationalacceleration, h is the local water depth, and c = √ gh is the long-wave speed.Shoaling in the context of steady solutions of the KdV equation is then described as follows.Given a steady wavetrain of wave height H , wavelength λ and period T , the variation of the waveas it propagates from depth h A to h B is obtained by imposing conservation of T , conservation of theenergy flux q E , and balance of forces using the bottom forcing by the pressure force, and the radiationstress.In order to execute this plan, one needs to have in hand formulations for the energy flux q E andthe radiation stress S xx in the context of the KdV equation. These expressions can be developed2igure 1: A periodic wave of wave height H and wavelength λ propagating at the free surface of a fluid ofundisturbed depth h . The free surface is given by z = η ( x, t ). using ideas first put forward in [3, 4]. Using the methods in these papers, it can be seen that we have q E = c (cid:16) h η + 54 h + h ηη xx (cid:17) and S xx = ρg (cid:16) η + 32 η + h η xx (cid:17) . Note that the energy flux is defined in a pointwise sense while the radiation stress is defined in anaverage sense. However, for the shoaling problem q E will also be averaged over one wave period.The plan of the paper is as follows. In the next section, the formulation of momentum and energybalance in the context of the KdV equation will be recalled, leading to the expression for flow force q I and energy flux q E . In Section 3, a formulation of the radiation stress consistent with the KdVequation will be found. Section 4 contains the formulation of the nonlinear shoaling equations, detailson the numerical implementation and a comparison with wave tank experiments with particular focuson the set-down. Section 6 details the comparison with other nonlinear shoaling theories withoutset-down, and in particular with the shoaling curves obtained by Svendsen and Buhr Hansen [45]. InSection 7, we provide a comparison with the experimental results of [45]. In the Conclusion we putour work into context and mention possibilities for further work. We start this section with a brief description of the water-wave problem of an inviscid, incompressibleand homogeneous fluid with a free surface. Due to the assumption that the relative change in waterdepth over one wavelength is less than the wave steepness, the problem can be locally formulated witha constant undisturbed depth h (see Figure 1). If the flow is assumed to be irrotational, the problemcan be formulated in terms of a velocity potential φ ( x, z, t ) in addition to the surface deflection η ( x, t ).The velocity field ( u ( x, z, t ) , v ( x, z, t )) is then given as the gradient of φ , and the pressure P can befound using the Bernoulli equation. In terms of φ and η , the problem is described by the Laplaceequation φ xx + φ zz = 0 for − h < z < η ( x, t ) (3)3n the fluid domain, and the boundary conditions η t + φ x η x − φ z = 0 on z = η ( x, t ) ,φ t + 12 ( φ x + φ z ) + gη = 0 on z = η ( x, t ) ,φ z = 0 on z = − h . It is well known that this problem is difficult to treat both mathematically and numerically.In particular, it is not known whether solutions exist on relevant time scales [22], and numericaldiscretization of the full problem on the field scale remains completely out of reach. Thus in practicalsituations, an asymptotic approximation of the Euler equations is usually required. The traditionalasymptotic theories are the linearized (Airy) theory and the hyperbolic (Saint-Venant) shallow-watertheory such as outlined in [42]. The work of Boussinesq [8] and Korteweg and de Vries [18] led toanother asymptotic regime, the so-called Boussinesq scaling which can be used for long waves of smallto moderate amplitude. In the present work, we consider waves in the Boussinesq regime, and since weconsider waves propagating in the shoreward direction only, we restrict attention to the Korteweg-deVries (KdV) equation. Indeed, the KdV equation describes unidirectional waves in the presence ofweakly nonlinear and dispersive effects in the case when these corrections are approximately balanced.In order to make this statement precise, it is useful to define the relative amplitude α = a/h andthe shallowness parameter β = h /λ , where a = H/ λ a representative wavelength of the wavefield to be studied. Using an asymptotic expansion and adimensional argument it can be shown that the KdV equation is a good model for long waves of smallamplitude if terms of order O ( α , αβ, β ) are neglected.In order to derive the KdV equation, one defines the non-dimensional variables x = λ ˜ x , z = h (˜ z − t = λc ˜ t , and φ = gaλc ˜ φ , then assumes that the velocity potential takes the form˜ φ = ˜ f − β ˜ z f ˜ x ˜ x + O ( β ) , (4)where ˜ f is the velocity potential evaluated at the bed. Following the description in [48], the free-surface boundary conditions can then be used to obtain the KdV equation (2). In the process, itcomes to light that the horizontal velocity at the bottom for a right-moving wave is given by˜ u = ˜ η − α ˜ η + 13 β ˜ η ˜ x ˜ x ˜ x + O ( α + β ) . (5)Moreover, as noted for example in [4], combining (5) and the derivative of (4) with respect to x one finds that the horizontal component of the velocity field is given by˜ φ ˜ x (˜ x, ˜ z, ˜ t ) = ˜ η − α ˜ η + β (cid:16) − ˜ z (cid:17) ˜ η ˜ x ˜ x + O ( α , αβ, β ) (6)in the KdV approximation. This relation will be needed later in the derivation of the the momentumand energy balance, and the formulation of the radiation stress. In addition, we need the verticalcomponent of the velocity field and the pressure. The vertical velocity is given by˜ φ ˜ z (˜ x, ˜ z, ˜ t ) = − β ˜ z ˜ η ˜ x + O ( αβ, β ) . In order to express the pressure in terms of the surface excursion η , assuming unit density for themoment, one first defines the dynamic pressure in the usual way: P ′ = P − P atm + gz, P ′ = ag ˜ P ′ . P ′ = ˜ η − β (˜ z − η ˜ x ˜ x + O ( αβ, β ) . (7)Approximations of the momentum and energy balance laws in the KdV scaling regime will now bederived. In the context of the full Euler equations with surface boundary conditions, the momentumbalance is expressed by ∂∂t Z η − h φ x dz + ∂∂x Z η − h (cid:8) φ x + P (cid:9) dz = 0 . (8)Using the previously defined scaling, the momentum balance is α ∂∂ ˜ t Z α ˜ η ˜ φ ˜ x h d ˜ z + ∂∂ ˜ x Z α ˜ η (cid:8) α ˜ φ x + α ˜ P ′ + (1 − ˜ z ) (cid:9) h d ˜ z = 0 . Substituting the two expressions ˜ φ x , ˜ P ′ and evaluating the integral one finds (cid:16) α ˜ η + 34 α ˜ η + 16 αβ ˜ η ˜ x ˜ x (cid:17) ˜ t + (cid:16)
12 + α ˜ η + 32 α ˜ η + 13 β ˜ η ˜ x ˜ x (cid:17) ˜ x = O ( α , αβ, β ) . From this relation we identify the non-dimensional momentum density˜ I = α ˜ η + 34 α ˜ η + 16 αβ ˜ η ˜ x ˜ x , and the non-dimensional momentum flux˜ q I = 12 + α ˜ η + 32 α ˜ η + 13 β ˜ η ˜ x ˜ x . Returning the expression to its dimensional forms through the scaling I = c h ˜ I and q I = c h ˜ q I yields I = c (cid:16) η + 34 h η + h η xx (cid:17) , and q I = c (cid:16) h η + 32 h η + h η xx (cid:17) . (9)Note that the expression for q I combines the momentum flux and the pressure force. FollowingBenjamin and Lighthill [6], we will use the term flow force for this quantity. Similarly, we can givethe energy balance in the full Euler equations as ∂∂t Z η − h n |∇ φ | + g ( z + h ) o dz + ∂∂x Z η − h n |∇ φ | + g ( z + h ) + P o φ x dz = 0 , and following the same procedure as above will lead to the expressions E = c (cid:16) h η + 14 h η + h ηη xx + h η x (cid:17) , and q E = c (cid:16) h η + 54 h η + h ηη xx (cid:17) , (10)for the energy density and energy flux respectively.5ote that q E is of second order, while the energy flux in the linear approximation is of first order.Indeed, it will be instructive to compare (10) to the well known expression for the energy flux in thelinear theory. Denoting the amplitude of a linear wave by a = H/ C g asbefore, it can shown (see for example [21]) that the energy flux averaged over one period is q linear E = gη C g = 12 ga C g . On the other hand, for cnoidal waves of very large wavelength and very small amplitude (10) can beapproximated by q E ∼ c h η = c gη . Averaging over one period and using the fact that the cnoidal wave is similar to a sinusoidal functionfor very small amplitudes, we obtain q E ∼ ga c . Finally, the two expressions are seen to be approximately equal if it is recognized that c is the groupvelocity for long linear waves in shallow water.Note that in previous works on cnoidal shoaling such as [44, 45], the linear formulation of theenergy flux was used. In addition to being obviously inconsistent, this approach also had practicaldrawbacks. Indeed as already mentioned, the method put forward in [45, 44] led to a discontinuityin the wave height at the matching point between linear and nonlinear theory. As will be shownmomentarily, keeping the nonlinear terms in the expression for q E resolves this problem, and indeedis essential in order to obtain the correct shoaling behavior for larger wave amplitudes.In the next section, we will use the definition of the flow force q I , and ideas outlined above toderive an expression for the radiation stress in the context of the KdV equation. Incorporating theradiation stress into the shoaling equations will then enable us to identify the set-down in additionto the wave height of a shoaling wave. In linear water-wave theory, the radiation stress is well understood as a second-order response to aperiodic wave train which is akin to the energy flux. Wave radiation stress is an important ingredientin the study of wave-current interactions, and also plays a major role in the understanding of wavesetup in the context of breaking waves. The reader may consult [27] for the mathematical developmentof an expression for the radiation stress in the context of linear waves. A discussion of the radiationstress based on a more heuristic physical understanding is given by Longuet-Higgins and Stewart [28].In the present work, we aim to expand the definition of the radiation stress to the nonlinear case, inparticular in the context of the KdV equation. Using the momentum balance equation derived in theprevious section, a nonlinear version of the radiation stress will be formulated. With this expressionin hand, some consequences of the radiation stress on wave shoaling will be investigated.Analyzing the definition given in [28], one can simply think of radiation stress as the total flowforce of a progressive wave averaged over one period minus the hydrostatic pressure force at rest.As previously discussed we consider a wave propagating solely in the x − direction, and neglect alltransverse effects. Then the definition of the principal component of the radiation stress is given by S xx = Z η − h ( ρu + P ) dz + Z − h ρgzdz. The first term expresses the total flux of momentum across a plane integrated from the bottom to thefree surface and with unit width, and as usual, the overbar denotes that the average over one wave6igure 2:
Schematic of the forces acting on a control volume given by a differential interval of length dx andreaching through the whole fluid column. period with respect to time is taken. The second term expresses the flow force in the absence of anywave motion.In the context of the KdV equation, the flow force is given by equation (9) when rescaled for afluid with unit density. Consequently, the x − component of the radiation stress is obtained in theKdV theory as S xx = q I + Z − h ρgzdz, = ρgh (cid:16) h h η + 32 h η + h η xx (cid:17) − ρgh , = ρg (cid:16) h η + 32 η + h η xx (cid:17) . (11)This formulation represents an improvement over the formulation given in [43] which only recordedthe middle term ρg η in the final expression for the radiation stress. A preliminary version of thisquantity was also presented in [19]. In order to apply the radiation stress to the formulation of theshoaling problem, let us consider a wave encountering a gently sloping beach. As indicated in Figure2, the momentum flux is reduced in the onshore direction due to an opposing horizontal force exertedby the bed, and opposing the fluid pressure. When accounting for the set-down of the mean surfacelevel, the flow force on the offshore side of the boundary of the control volume is given in terms ofradiation stress by q I = S xx + Z η − h ρg ( η − z ) dz = S xx + 12 ρg ( η + h ) . This relation can be found for example in [12, 28]. As already mentioned above, the momentum fluxis reduced by the reaction force exerted by the sea bed on the fluid due to the weight of the fluid.This force is equal in magnitude to the horizontal component of the pressure force at the bottom ρg ( h + η ) dz , where dz denotes the vertical variation in the sea bed over the horizontal distance dx .Thus the horizontal component of the pressure force exerted on the bottom can be written as ρg ( h + η ) dz = ρg ( h + η ) h x dx. (12)As indicated in Figure 2, the difference in flow force between the shoreward face and the offshore faceof the control volume is given by ddx h S xx + 12 ρg ( h + η ) i dx. (13)7vidently, momentum balance requires (12) and (13) to be equal, and simplifying yields the relation dS xx dx + ρg ( η + h ) dηdx = 0 . (14)In order to develop the shoaling equation, we integrate over the control interval [ x, x + dx ] to get Z x + dxx dS xx dx + Z x + dxx ρgη dηdx + Z x + dxx ρgh dηdx = 0 . The first two integrals are straightforward to compute. The third integral is evaluated by usingintegration by parts and then approximated using the trapezoidal rule. Defining the difference of anyquantity F as ∆ F = F | x + dx − F | x we get, the relation∆ S xx = − ρg η − ρg h + h x dx )∆ η. (15)Recall that we are interested in the averaged behavior of the wave height for a wave-train approachingthe beach. We therefore find it useful to define the changes in radiation stress averaged over a periodby ∆ S xx = − T Z T n ρg η + ρg h + h )∆ η o dt, (16)where S xx is defined by (11). With this identity in hand, it will be possible to include the developmentof the set-down η in the shoaling problem. The formulation of the complete shoaling equation will bethe subject of the next section. An explicit set of equations describing the shoaling of long waves on a gently sloping beach will nowbe given. The idea is to use the well known cnoidal wave solution for periodic waves in the KdVequation together with the assumption that the bottom slope is gentle enough so that reflections canbe neglected, and the waves are able to adjust adiabatically to the new local depth. The wave profile η will distort, but the underlying cnoidal shape and periodicity are preserved.A solution of the KdV equation for periodic waves of constant form was first discovered by Ko-rteweg and de Vries in 1895 [18]. Assuming a traveling-wave solution of constant shape, one maymake the ansatz η ( x, t ) = f ( x − ct ), and reduce the KdV equation to an ordinary differential equationof the form − h f ′ ) = F ( f ) = 14 c h f + (cid:16) c − c (cid:17) f + Af + B, (17)with A and B being constants of integration. This differential equation has a number of solutions,but for practical applications only periodic real-valued and bounded solutions correspond to realisticwave profiles. Since F ( f ) is a third-order polynomial it has three roots. Periodic solutions occur onlyin the case where the roots are distinct, so that they can be labeled f < f < f , and this conventionallows us to write F ( f ) = ( f − f )( f − f )( f − f ) . Graph of the function − F ( f ). The roots are labelled by f < f < f . The red lines denote thevalues of f such that − F ≥ Requiring the solution to be real, it can be seen from (17) that − F ( f ) must be positive. Examiningthe phase plane plot in Figure 3, it is clear that periodic solutions will oscillate between f and f .Since f < f by assumption, f will denote the wavecrest while f will be the trough. Consequently,the wave height H is given by the difference H = f − f . As is well known, the periodic solutions of(17) are given in terms of the Jacobian elliptic function cn . The solution of the KdV equation is then η ( x, t ) = f − ( f − f )cn (cid:16) K ( m ) (cid:0) tT − xλ (cid:1) , m (cid:17) , (18)where m is the elliptic modulus defined by m = f − f f − f , and K ( m ) is the complete elliptic integral offirst kind [24]. The wavespeed c and the wavelength λ are then given by c = c (cid:16) f + f + f h (cid:17) , λ = K ( m ) s h f − f ) . (19)In order to use these formulas in practice, it is convenient to take the wave height H , the mean surfacelevel η and the elliptic parameter m as given parameters. The roots of F can then by computedexplicitly and are given as follows. f = η − HE ( m ) mK ( m ) ,f = f + Hm ,f = f − H. (20)The shoaling problem can be formulated as follows. Suppose we know that a periodic wave train ofwave height H , wavelength λ and mean-surface level η is given at a local depth h A . Assuming that thewaves can be approximately described by the the KdV equation, a unique cnoidal wave solution canbe found, and the frequency ν , the energy flux q E and the radiation stress S xx be found. The shoalingproblem now consists of finding an appropriate cnoidal solution at a smaller depth h B . Relying onthe conservation of frequency and energy flux, and a prescribed change in the radiation stress, theshoaling equations can be written as ν B ( f , f , f , h B ) = ν A = c A λ A , (21) q BE ( f , f , f , h B ) = q AE = 1 T Z T n h A η + 54( h A ) η + h A ηη xx o dt, (22) S Bxx ( f , f , f , h B ) = S Axx − T Z T n ρg η + ρg h B + h A )∆ η o dt. (23)9he first equation is conservation of frequency ν , while the second equation expresses conservationof energy flux integrated over a period T [4]. Finally, the third equation is (16), and is indeed aformulation of conservation of momentum expressed in terms of radiation stress. In previous work[20], a similar system was found, and then solved for the parameters f , f and f . However, themethod used in [20] had several drawbacks. First, due to a lack of a definition of radiation stress,the equation for momentum conservation was replaced by conservation of the mean fluid level. Thisapproach was therefore not able to predict the wave set-down. Moreover, solving for the parameters f , f and f was not optimal from a numerical point of view, and the algorithm terminated longbefore the highest wave was reached. In the present work, we use the explicit form of the parametersto solve directly for H , m and η . This approach is more expedient from a numerical point of view, andallows us to go much higher up on the shoaling curve. In terms of the wave height H , the conservationof frequency turns into the following cubic equation for H : − g ( E ( m ) mK ( m ) − m + 1) K ( m ) h m H + 3 g ( η h + 1)( E ( m ) mK ( m ) − m + 1)16 K ( m ) h m H + − g ( η h + 1) K ( m ) h m H + 1 T = 0 . (24)Similarly we can manipulate the equation describing the changes in radiation stress to find an expres-sion for the set-down. Indeed, (23) reduces to the quadratic equation − η + A η + B = 0 , (25)with A = (cid:16) H − h B h A − H cn ( m ) − Hm + 3 HE ( m ) mK ( m ) (cid:17) and B = S Axx ( m ) − η A ( m )2 − H cn ( m ) − ( h B ) η xxB ( m )3 + 32 h A η A ( m ) − h B η A ( m ) −
32 ( H − Hm + HE ( m ) mK ( m ) ) + 3 H cn ( m )( H − Hm + HE ( mmK ( m ) ) . Note that we had to integrate various powers and derivatives of cn ( ξ ; m ). These formulas arebased on calculations that can be found in [24] and [1], and are given in explicit form in the appendix.Having this representation in hand, we can in principle write (24) as H = F ( m, η ) and take (25) to be η = G ( m, H ) in explicit terms as two coupled equations. To solve these two equations simultaneously,we may iterate between them at current local depth h B as follows; first initialize the procedure with η ( m ) and then find H by H ( m ) i +1 = F ( m, η i ( m )) , solving the cubic equation. This can in turn be used to update the set-down by η ( m ) i +1 = G ( m, H i +1 ( m )) . Repeating this process, we continue to approximate H and η until a stopping criterion has beenreached. Using this reduction allows us to solve the nonlinear system (21),(22),(23), by the followingprocedure. From the value of m at depth h A , we increment up, at each step iterating the aboveequations for H and η , and then checking whether (22) is satisfied to a specified tolerance. If that isthe case, we consider the system solved at depth h B , and move on to the next step with depth h C ,thereby moving up the slope. 10 Implementation of the shoaling equations
Having explained the nonlinear shoaling equations, and the strategy for finding approximate solu-tions, we will now implement the equations and produce shoaling curves for various deep-water data.The method used here proceeds in three stages, such as originally developed in [44, 45]. Suppose anincoming wave of wavelength λ is given at a starting depth h . While this depth may not be deep water, it is assumed that h > . λ , so that we label this starting value the deep-water wavelength.First, the linear shoaling equation is used up to the point h/λ = 0 .
1. At this point the cnoidaltheory is valid [39], and we propose a matching technique to obtain the fundamental parameters ofthe nonlinear wave. Lastly, we use the nonlinear KdV theory to follow the shoaling curve into thenonlinear region.
Stage 1. Linear theory:
The essential ingredient in this step is the linear dispersion relation ω − gk tanh ( kh ) = 0 , (26)and the conservation of energy. Since the circular frequency ω is conserved, the nonlinear equation ω − gk B tanh ( k B h B ) = 0 needs to be solved numerically for the wave number k B using a Newtonsolver. The wave height is then found from the conservation of energy as H B = H A s C Ag C Bg , (27)where C g = dωdk is the group velocity depending on the wave number and depth. The derivation ofequation (27) is described in detail by for example in [12]. The set-down can be solved similarly usingthe procedure outlined in [28] using the radiation stress. Stage 2. Matching linear and nonlinear theory:
Due to increasing waveheight, the nonlinearalgorithm must be utilized on the final stretch of the shoaling curve. The transition must be madein a region where both theories are valid, and some general conditions on the transition depth aregiven in [45]. In order to initialize the nonlinear solver at the transition depth, three parameters needto be specified. We choose to match the wave height H and the set-down η with the linear theory.In addition one of the parameters λ , c , ν , q E can be be matched. To do this we choose one of theseparameters and then keeping H and η constant, we find an elliptic modulus m ∈ (0 ,
1) which leads toequality in the chosen parameter. In principle we can only match one parameter, and it is not clearwhich one will be the most convenient. To investigate this issue, we define the root problems λ lin − λ nonlin ( m ) = 0 , c lin − c nonlin ( m ) = 0 ,ν lin − ν nonlin ( m ) = 0 , q linE − q nonlinE ( m ) = 0 . Here the linear parameter is a fixed constant while the nonlinear quantity is a function of m given that H and η are given. Plotting the root problems we observe that the wavelength is most sensitive tovarying m , and therefore is the best candidate for a matching parameter. Figure 4 shows a particularexample. Stage 3. Cnoidal shoaling:
The final step in the shoaling procedure is solving for wave height andset-down using the scheme defined in Section 6. First, define H and η as functions of m as given byformula (24) and (25). Then use a nonlinear solver to find m from equation (22). With m in hand,one can determine the wave height H ( m ) at the current local depth. Repeating this procedure willthen determine changes of wave height and set-down at consecutive points, allowing us to move upthe slope. 11igure 4: Root problems defined by the parameters λ, c, ν, q E as functions of m . Having defined the shoaling equations and the numerical approach, we may plot the developmentof the wave height and set-down for a practical shoaling problem. A comparison of our numericalresults with the wave tank data obtained by Bowen et al. [9] and the classical linear theory is shownin Figure 5. In this figure, the red circles represent the wave height of the shoaling wave according tothe experimental data of [9], and the red asterics represent the experimentally determined set-downand set-up, i.e. the deviation of the mean depth from the still water line (S.W.L.). The dashed bluecurves show the prediction of the linear theory as already presented in [9]. Finally, the solid curvesrepresent the shoaling model put forward in the present article. The waves under consideration hadan initial wave length of λ = 202cm and an initial wave height H = 6 . nonlinearmid-elevation rise was acknowledged in [25] has been found previously with higher-order methods.For example, a combination of third-order hyperbolic waves near the shore and Stokes waves furtheroffshore was used in [17], and Cokelet’s extension of Stokes’ approximation of periodic waves was usedin [41]. In both works, a similar rise of the nonlinear set-down was observed. In the present work,the computed set-down matches the experimental data even well beyond the breaking point, but wehave to acknowledge that the KdV model ceases to be valid for breaking waves. It was argued in [43] that the wave set-down is negligible in many practical cases. With this as-sumption, the shoaling problem can be simplified since only two equations need to be solved. Thisapproach has been taken in many works on wave shoaling.The particular aim of this section is to improve the method put forward in [44]. As already men-tioned, the use of a first-order approximation of the energy flux [44] gave rise to a discontinuity in waveheight at the matching point between the linear and their nonlinear theory due to the inconsistency ofthe approximation of the energy flux. A partial fix was presented in [45] where a continuous shoalingcurve is obtained by use of conversion tables but at the cost of not conserving the energy flux. Morerecently, the shoaling theory was extended in [20] using the nonlinear definition of the energy flux12igure 5:
Profile of the mean water level η and the wave height H as functions of the depth, compared todata points taken from [9]. In this case, we have wave period T = 1 .
14s incident wave height H = 6 . H b = 8 . β ) = 0 . developed in [3, 4]. However, due to the specific way the problem was formulated, it was not possibleto reach all the way up the shoaling curve due to numerical instabilities.In the present section, it will be shown that with the definition of the shoaling equation fromSection 4, and the implementation explained in Section 5, we can find a continuous shoaling curvethat reaches all the way up to the highest wave and agrees well with the more accurate methods of[36] using Cokelet’s theory and that of [46] using the fully nonlinear steady Euler equations.The algorithm is the same as in Section 5, except that in equation (25) we set η = 0 at each step.We then use the same three stages as explained in Section 5 to determine the development of thewave height in the shoaling region. First formulate H = H ( m ) according to formula (24) given inSection 6. Then use the assumption of zero set-down to find the roots given by (20). Having the rootsas functions of m we may use energy conservation to define a nonlinear equation as done in equation(22). Solving for m we are free to determine the wave height at a specified depth h .We note that the present implementation is able to determine the wave height further into theshoaling region as compared to the curves presented by [20]. This is due to the simplicity of theimplementation, solving two nonlinear equations rather than three. We also see in Figure 6 that ourtheory is in fairly good agreement with the shoaling curves presented in [44] which are in reasonableagreement with higher-order theories. Moreover, it is clearly visible in Figure 6 that our curves donot feature a discontinuity in wave height.Finally, comparisons of shoaling curves computed using the method at hand with data provided bywave tank experiments are presented. The experimental data are taken from [45] where considerablepains were taken to stay within the correct scaling regime. The waves were generated with a piston-type wavemaker, then propagated over a flat bottom with still water depth of 36 cm before shoaling13igure 6: Shoaling curves based on present theory in blue compared to the theory presented by Svendsenand Brink-Kjær [44] in red. Deep water values H /λ = { . , . , . , . , . } . on a 1 : 35 beach. We initialize the code with wave height and wavelength at depth h = 36cm at thetoe of the slope. The code then computes wave height H with the linear model up to h/λ > . H /λ = 0 . In this paper we have shown how to derive expressions for the energy flux and radiation stress of awavetrain in the context of the well known KdV equation. The derivations are based on a generalapproach developed in [3, 4], and yield approximations of the same order as the KdV equation itself.Compared to previous definitions of such quantities as for example in [43, 44], these expressions featureadditional terms which are required by the theoretical underpinnings of the method used here.The energy flux and radiation stress are used to develop shoaling equations for long waves of smallto moderate amplitude, and since such waves can be described by cnoidal functions, the shoalingequations can be posed as a 3 × Comparisons of shoaling curves from six experiments with a range of initial wave steepness.
Acknowledgments
The authors thank Volker Roeber for helpful discussions. This work has received funding from the EuropeanUnion’s
Horizon 2020 research and innovation programme under grant agreement No. 763959.
A Integrals of cnoidal functions
In order to define the two shoaling models we need to handle integrals of various powers of derivativesof η . In particular we need to determine η , η , η , η xx and ηη xx in order to evaluate (25) in terms of m and the current depth h . First note that we can write the time-averaged η given by (18) in themore convenient form η = 1 T Z T η (cid:16) K ( m )( tT − xλ ) (cid:17) dt = f + H Z cn (2 K ( m ) ξ ; m ) dξ, (for more details see [43]). Similarly considerations apply to the powers of η which involve integralsof the form Z cn (2 Kξ ; m ) dξ = 14 mK (cid:16) E − (1 − m ) K (cid:17) , Z cn (2 Kξ ; m ) dξ = 13 m (cid:16) m − m + 2(4 m − EK (cid:17) , Z cn (2 Kξ ; m ) dξ = 15 m (cid:16) m − Z cn (2 Kξ ; m ) dξ + 3(1 − m ) Z cn (2 Kξ ; m ) dξ (cid:17) . η xx and ηη xx also include such terms due to theidentity η xx (2 Kξ ; m ) = 4 K H (2 − m + (8 m − (2 Kξ ; m ) − m cn (2 Kξ ; m )) , and similar relations which can be found in [1]. Combining the results above we deduce that η = f + H Z cn (2 Kξ ; m ) dξ,η = f + 2 Hf Z cn (2 Kξ ; m ) dξ + H Z cn (2 Kξ ; m ) dξ,η = H Z cn (2 Kξ ; m ) dξ + 3 H f Z cn (2 Kξ ; m ) dξ + 3 Hf Z cn (2 Kξ ; m ) dξ + f ,η xx = 3 H m ( − m Z cn (2 Kξ ; m ) dξ + 4 m Z cn (2 Kξ ; m ) dξ − Z cn (2 Kξ ; m ) dξ − m + 1) ,ηη xx = 3 H f Z cn (2 Kξ ; m ) dξ − H f Z cn (2 Kξ ; m ) dξ + 32 Hf f − Hf f + 32 H f Z cn (2 Kξ ; m ) dξ − H f Z cn (2 Kξ ; m ) dξ + 6 H f m Z cn (2 Kξ ; m ) dξ − H f m Z cn (2 Kξ ; m ) dξ − Hf f m + 32 3 Hf f m − H f m Z cn (2 Kξ ; m ) dξ + 32 H f m Z cn (2 Kξ ; m ) dξ − Hf f Z cn (2 Kξ ; m ) dξ + 3 Hf f Z cn (2 Kξ ; m ) dξ − H f m Z cn (2 Kξ ; m ) dξ + 92 H f m Z cn (2 Kξ ; m ) dξ + 6 Hf f m Z cn (2 Kξ ; m ) dξ − Hf f m Z cn (2 Kξ ; m ) dξ − Hf f m Z cn (2 Kξ ; m ) dξ + 92 Hf f m Z cn (2 Kξ ; m ) dξ. eferences [1] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, andMathematical Tables, National Bureau of Standards Applied Mathematics series 55, (1965).[2] A. Ali and H. Kalisch, Energy balance for undular bores , Compt. Rend. Mecanique (2010), 67–70.[3] A. Ali and H. Kalisch,
Mechanical balance laws for Boussinesq models of surface water waves , J. Non-linear Sci. (2012), 371–398.[4] A. Ali and H. Kalisch, On the formulation of mass, momentum and energy conservation in the KdVequation , Acta Applicandae Mathematicae (2014), 113–131.[5] J.A. Battjes and H.W. Groenendijk,
Wave height distribution on shallow foreshores , Coastal Engineering (2000), 161–182. no. 3[6] B.T. Benjamin and J. Lighthill, On cnoidal waves and bores , Proc. R. Soc. A (1954), 448–460.[7] M. Bjørkav˚ag and H. Kalisch,
Wave breaking in Boussinesq models for undular bores , Phys. Lett. A (2011), 1570–1578.[8] J. Boussinesq,
Th´eorie des ondes et des remous qui se propagent le long d’un canal rectangulaire horizon-tal, en communiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surfaceau fond , J. Math. Pures Appl. (1872), 55–108.[9] A.J. Bowen, D.L. Inman and V.P. Simmons, Wave ‘set-down’ and set-up , J. Geophysical Research (1968), 2569–2577.[10] M. Brocchini, A reasoned overview on Boussinesq-type models: the interplay between physics, mathe-matics and numerics , Proc. R. Soc. A (2013), 1–27.[11] J. Buhr Hansen and I.A. Svendsen,
Laboratory generation of waves of constant form , Proc. Coast. Eng.Conf., 14th, Cope, 1974.[12] R.G. Dean and R.A. Dalrymple, Water Wave Mechanics for Scientists and Engineers, (World Scientific,Singapore, 1984).[13] S. Grilli, R. Subramanya, I.A. Svendsen and J. Veeramony,
Shoaling of solitary waves on plane beaches ,J. Wtrwy., Port, Coast., and Oc. Engrg. (1994), 609–628.[14] L.H. Holthuijsen, Waves in Oceanic and Coastal Waters. (Cambridge University Press, Cambridge,2010).[15] S. Israwi and H. Kalisch,
Approximate conservation laws in the KdV equation , Phys. Lett. A (2019),854–858.[16] S. Israwi, H. Kalisch,
A mathematical justification of the momentum density function associated to theKdV equation , in press.[17] I.D. James, 1974.
Non-linear waves in the nearshore region: shoaling and set-up , Estuarine CoastalMarine Science (1974), 207–234.[18] D.J. Korteweg and G. de Vries, On the change of form of long waves advancing in a rectangular channeland on a new type of long stationary wave , Philos. Mag. (1895), 422–443.[19] Z. Khorsand, Flow Properties of Fully Nonlinear Model Equations for Surface Waves , The Universityof Bergen, (2017).[20] Z. Khorsand and H. Kalisch,
On the shoaling of solitary waves in the KdV equation , Coastal EngineeringProceedings, (2014).[21] P. Kundu and I. Cohen, Fluid Mechanics (Academic Press, San Diego, 2004).[22] D. Lannes, The Water Waves Problem. Mathematical Surveys and Monographs, vol. (Amer. Math.Soc., Providence, 2013).[23] D. Lannes and P. Bonneton, Derivation of asymptotic two-dimensional time-dependent equations forsurface water wave propagation , Phys. Fluids (2009), 016601.[24] D.F. Lawden, Elliptic Functions and Applications (Springer, New York, 2013).[25] B. Le M´ehaut´e, An Introduction to Hydrodynamics and Water Waves (Springer, New York, 1976).[26] J. Lighthill, Waves in Fluids (Cambridge University Press, Cambridge, 1978).[27] M.S. Longuet-Higgins and R.W. Stewart, Radiation stress and mass transport in gravity waves, withapplication to ‘surf beats’ , J. Fluid Mech. (1962), 481–504.
28] M.S. Longuet-Higgins and R.W. Stewart,
Radiation stresses in water waves; a physical discussion, withapplications , Deep Sea Research (1964), 529–562.[29] P.A. Madsen D.R. Fuhrman and B. Wang, A Boussinesq-type method for fully nonlinear waves inter-acting with rapidly varying bathymetry , Coast. Eng. (2006), 487–504.[30] P.A. Madsen and H.A. Sch¨affer, Higher-order Boussinesq-type equations for surface gravity waves:derivation and analysis , Phil. Trans. Roy. Soc. Lond. Ser. A (1998), 3123–3184.[31] O. Nwogu,
Alternative form of Boussinesq equations for nearshore wave propagation , J. Waterway, Port,Coastal and Ocean Engineering (1993), 618–638.[32] L.A. Ostrovsky and E.N. Pelinovsky,
Wave transformation on the surface of a fluid of variable depth ,Atmospheric and Oceanic Physics (1970), 552–555.[33] D.H. Peregrine, Long waves on a beach , J. Fluid Mech. (1967), 815–827.[34] Lord Rayleigh, Hydrodynamical notes , Phil. Mag. (1911), 177–195.[35] V. Roeber, K.F. Cheung and M.H. Kobayashi, Shock-capturing Boussinesq-type model for nearshorewave processes , Coastal Engineering (2010), 407–433.[36] T. Sakai and J.A. Battjes, Wave shoaling calculated from Cokelet’s theory , Coastal Engineering, (1980),65–84.[37] T. Saville, Experimental determination of wave set-up , Proc. 2nd Tech. Conf. on Hurricanes, MiamiBeach, FL, (1961).[38] F. Serre,
Contribution ´a l’ ´etude des ´ecoulements permanents et variables dans les canaux , La HouilleBlanche, (1953), 374–388.[39] O. Skovgaard, Sinusoidal and cnoidal gravity waves formulae and tables , Institute of Hydrodynamicsand Hydraulic Engineering, Technical University of Denmark, (1974).[40] R.M. Sorensen, Basic Wave Mechanics for Coastal and Ocean Engineers (John Wiley & Sons, New York,1993).[41] M. Stiassnie and D.H. Peregrine,
Shoaling of finite-amplitude surface waves on water of slowly-varyingdepth , J. Fluid Mech. (1980), 783–805.[42] J.J. Stoker, Water Waves: The Mathematical Theory with Applications, (Interscience Publishers, NewYork, 1957).[43] I.A. Svendsen, Introduction to Nearshore Hydrodynamic, (World Scientific, 2006).[44] I.A. Svendsen and O. Brink-Kjær, Shoaling of cnoidal waves , Proc. 13th Intern. Conf. Coastal Engrg.,Vancouver (1972), 365–383.[45] I.A. Svendsen and J. Buhr Hansen,
The wave height variation for regular waves in shoaling water ,Coastal Engineering, (1977), 261–284.[46] R.H. Swift and J.C. Dixon, Transformation of regular waves , Proceedings of the Institution of CivilEngineers (1987), 359–380.[47] G. Wei, J.T. Kirby, S.T. Grilli and R. Subramanya, A fully nonlinear Boussinesq model for surfacewaves. Part 1. Highly nonlinear unsteady waves , J. Fluid Mech. (1995), 71–92.[48] G. B. Whitham, Linear and Nonlinear Waves (Wiley, New York, 1974).[49] J.M. Witting,
A unified model for the evolution of nonlinear water waves , J. Comp. Phys. (1984),203–236.(1984),203–236.