A Nonlinear Formulation of Radiation Stress and Applications to Cnoidal Shoaling

In this article we provide formulations of energy flux and radiation stress consistent with the scaling regime of the Korteweg-de Vries (KdV) equation. These quantities can be used to describe the shoaling of cnoidal waves approaching a gently sloping beach. The transformation of these waves along the slope can be described using the shoaling equations, a set of three nonlinear equations 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 verify our shoaling formulation by comparing with experimental data from two sets of experiments as well as shoaling curves obtained in previous works.


Introduction
In the present work, the development of surface waves across a gently sloping bottom is in view. Our main goal is the prediction of the wave height and the set-down of a periodic wave as it enters an area of shallower depth. This problem is known as the shoaling problem and has a long history, with contributions 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 assumptions. First, the bottom slope is small enough to allow the waves to continuously adjust to the changing depth without altering their basic shape or breaking up. Second, reflections are assumed to be negligible so that the energy flux of a wave is conserved as it transforms on the slope. Third, the period T of a wave is assumed to be constant as it progresses. This last point is sometimes termed the conservation of waves as the number of waves in a group will remain constant during the shoaling process.
Rather than following a wave in space and time as it propagates over the slope, we estimate wave properties as functions of the depth using conservation of period T and energy flux, and a prescribed change in momentum due to the effect of the radiation stress. These conservation equations are due to the assumptions stated above, and they form the basis for the description of the shoaling process as first envisioned by Rayleigh [34], and subsequently used by a number of authors. The crucial factor which guarantees that the shoaling assumptions are valid is that the relative change in water depth over 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 wave shoaling and can be found in textbooks on coastal engineering such as [12,40]. This approach depends on the assumption that the waves can be described by a sinusoidal function of the form H 2 cos(kx−ωt), where H is the wave height, k is the wavenumber, ω is the circular frequency, t is the time and x is a spatial coordinate along which the bottom is sloping. Recall that the wavelength is λ = 2π k and the wave period is T = 2π ω , and these are related by the dispersion relation ω 2 = gk tanh (kh) at local depth h. In the linear shoaling theory, the conservation of energy flux is central. The energy flux is generally 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 and the wave height H as E = 1 8 ρgH 2 , 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 adjust adiabatically to the changing conditions, and that reflections and distortions are negligible. The wave-height transformation of a shoaling wave is obtained by imposing conservation of the energy flux EC g across the shoaling region. If a wave measurement H 0 , T 0 at some point offshore is given, then the wave height at a point closer to shore may be computed from the conservation of energy flux given in the form 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 the set-down [28]. The authors of [44] used cnoidal functions in connection with the linear formulation of the energy flux defined above to obtain a hybrid theory of shoaling. Unfortunately, this approach led to a discontinuity in wave height at the matching point between the linear shoaling equation (1) and the nonlinear theory based on the cnoidal functions. The problem was remedied to some degree in [45] by requiring continuity in wave height effected through the use of conversion tables. This approach led to good agreement with the experimental data, but was cumbersome to implement, and as already noted in [45], the continuity of the wave height at the matching point in the shoaling curve led to a discontinuity in energy flux.
The problem can be resolved by defining a shoaling equation based on the full water-wave problem using 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 equations leading to large systems of equations which need to be solved numerically. If the approach is to be used in connection with ocean-wave statistics, the approach based on the steady Euler equations might be too computationally demanding. Moreover, the agreement of the combined linear-cnoidal approach 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 and conservation of energy flux can be resolved in the context of the KdV theory by using a definition of the energy flux consistent with the KdV equation [4,15], and without using the fully nonlinear approach of [46]. Recall that the KdV equation is given in the form where η(x, t) is the deflection of the free surface from rest at a point x and a time t, g is the gravitational acceleration, h 0 is the local water depth, and c 0 = √ gh 0 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 wave as it propagates from depth h A to h B is obtained by imposing conservation of T , conservation of the energy flux q E , and balance of forces using the bottom forcing by the pressure force, and the radiation stress.
In order to execute this plan, one needs to have in hand formulations for the energy flux q E and the radiation stress S xx in the context of the KdV equation. These expressions can be developed using ideas first put forward in [3,4]. Using the methods in these papers, it can be seen that we have Note that the energy flux is defined in a pointwise sense while the radiation stress is defined in an average 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 energy balance 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 KdV equation will be found. Section 4 contains the formulation of the nonlinear shoaling equations, details on the numerical implementation and a comparison with wave tank experiments with particular focus on the set-down. Section 6 details the comparison with other nonlinear shoaling theories without set-down, and in particular with the shoaling curves obtained by Svendsen and Buhr Hansen [45]. In Section 7, we provide a comparison with the experimental results of [45]. In the Conclusion we put our work into context and mention possibilities for further work.

Momentum and energy balance in the KdV approximation
We start this section with a brief description of the water-wave problem of an inviscid, incompressible and homogeneous fluid with a free surface. Due to the assumption that the relative change in water depth over one wavelength is less than the wave steepness, the problem can be locally formulated with a constant undisturbed depth h 0 (see Figure 1). If the flow is assumed to be irrotational, the problem can 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 be found using the Bernoulli equation. In terms of φ and η, the problem is described by the Laplace equation in the fluid domain, and the boundary conditions 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 numerical discretization of the full problem on the field scale remains completely out of reach. Thus in practical situations, an asymptotic approximation of the Euler equations is usually required. The traditional asymptotic theories are the linearized (Airy) theory and the hyperbolic (Saint-Venant) shallow-water theory such as outlined in [42]. The work of Boussinesq [8] and Korteweg and de Vries [18] led to another asymptotic regime, the so-called Boussinesq scaling which can be used for long waves of small to moderate amplitude. In the present work, we consider waves in the Boussinesq regime, and since we consider waves propagating in the shoreward direction only, we restrict attention to the Korteweg-de Vries (KdV) equation. Indeed, the KdV equation describes unidirectional waves in the presence of weakly 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 0 and the shallowness parameter β = h 2 0 /λ 2 , where a = H/2 denotes a representative amplitude and λ a representative wavelength of the wavefield to be studied. Using an asymptotic expansion and a dimensional argument it can be shown that the KdV equation is a good model for long waves of small amplitude if terms of order O(α 2 , αβ, β 2 ) are neglected.
In order to derive the KdV equation, one defines the non-dimensional variables , and φ = gaλ c 0φ , then assumes that the velocity potential takes the form wheref is the velocity potential evaluated at the bed. Following the description in [48], the freesurface boundary conditions can then be used to obtain the KdV equation (2). In the process, it comes to light that the horizontal velocity at the bottom for a right-moving wave is given bỹ 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 bỹ in the KdV approximation. This relation will be needed later in the derivation of the the momentum and energy balance, and the formulation of the radiation stress. In addition, we need the vertical component of the velocity field and the pressure. The vertical velocity is given bỹ In order to express the pressure in terms of the surface excursion η, assuming unit density for the moment, one first defines the dynamic pressure in the usual way: Then using the Bernoulli equation and following the computations outlined in [4] leads to the nondimensional perturbation pressure in the KdV approximation: Approximations of the momentum and energy balance laws in the KdV scaling regime will now be derived. In the context of the full Euler equations with surface boundary conditions, the momentum balance is expressed by Using the previously defined scaling, the momentum balance is Substituting the two expressionsφ x ,P ′ and evaluating the integral one finds From this relation we identify the non-dimensional momentum densitỹ and the non-dimensional momentum flux Returning the expression to its dimensional forms through the scaling I = c 0 h 0Ĩ and q I = c 2 0 h 0qI yields Note that the expression for q I combines the momentum flux and the pressure force. Following Benjamin and Lighthill [6], we will use the term flow force for this quantity. Similarly, we can give the energy balance in the full Euler equations as and following the same procedure as above will lead to the expressions and for the energy density and energy flux respectively.
Note 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 the linear theory. Denoting the amplitude of a linear wave by a = H/2 and the group velocity by C g as before, it can shown (see for example [21]) that the energy flux averaged over one period is On the other hand, for cnoidal waves of very large wavelength and very small amplitude (10) can be approximated by Averaging over one period and using the fact that the cnoidal wave is similar to a sinusoidal function for very small amplitudes, we obtain Finally, the two expressions are seen to be approximately equal if it is recognized that c 0 is the group velocity for long linear waves in shallow water.
Note that in previous works on cnoidal shoaling such as [44,45], the linear formulation of the energy flux was used. In addition to being obviously inconsistent, this approach also had practical drawbacks. Indeed as already mentioned, the method put forward in [45,44] led to a discontinuity in the wave height at the matching point between linear and nonlinear theory. As will be shown momentarily, keeping the nonlinear terms in the expression for q E resolves this problem, and indeed is 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 to derive an expression for the radiation stress in the context of the KdV equation. Incorporating the radiation stress into the shoaling equations will then enable us to identify the set-down in addition to the wave height of a shoaling wave.

Radiation stress in the KdV approximation
In linear water-wave theory, the radiation stress is well understood as a second-order response to a periodic wave train which is akin to the energy flux. Wave radiation stress is an important ingredient in the study of wave-current interactions, and also plays a major role in the understanding of wave setup in the context of breaking waves. The reader may consult [27] for the mathematical development of an expression for the radiation stress in the context of linear waves. A discussion of the radiation stress 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, in particular in the context of the KdV equation. Using the momentum balance equation derived in the previous section, a nonlinear version of the radiation stress will be formulated. With this expression in 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 flow force 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 all transverse effects. Then the definition of the principal component of the radiation stress is given by The first term expresses the total flux of momentum across a plane integrated from the bottom to the free surface and with unit width, and as usual, the overbar denotes that the average over one wave period with respect to time is taken. The second term expresses the flow force in the absence of any wave motion.
In the context of the KdV equation, the flow force is given by equation (9) when rescaled for a fluid with unit density. Consequently, the x−component of the radiation stress is obtained in the KdV theory as This formulation represents an improvement over the formulation given in [43] which only recorded the middle term ρg 3 2 η 2 in the final expression for the radiation stress. A preliminary version of this quantity was also presented in [19]. In order to apply the radiation stress to the formulation of the shoaling problem, let us consider a wave encountering a gently sloping beach. As indicated in Figure  2, the momentum flux is reduced in the onshore direction due to an opposing horizontal force exerted by the bed, and opposing the fluid pressure. When accounting for the set-down of the mean surface level, the flow force on the offshore side of the boundary of the control volume is given in terms of radiation stress by This relation can be found for example in [12,28]. As already mentioned above, the momentum flux is 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 As indicated in Figure 2, the difference in flow force between the shoreward face and the offshore face of the control volume is given by Evidently, momentum balance requires (12) and (13) to be equal, and simplifying yields the relation In order to develop the shoaling equation, we integrate over the control interval [x, x + dx] to get The first two integrals are straightforward to compute. The third integral is evaluated by using integration by parts and then approximated using the trapezoidal rule. Defining the difference of any quantity F as ∆F = F | x+dx − F | x we get, the relation Recall that we are interested in the averaged behavior of the wave height for a wave-train approaching the beach. We therefore find it useful to define the changes in radiation stress averaged over a period by where S xx is defined by (11). With this identity in hand, it will be possible to include the development of the set-down η in the shoaling problem. The formulation of the complete shoaling equation will be the subject of the next section.

The nonlinear shoaling equations
An explicit set of equations describing the shoaling of long waves on a gently sloping beach will now be given. The idea is to use the well known cnoidal wave solution for periodic waves in the KdV equation together with the assumption that the bottom slope is gentle enough so that reflections can be 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 Korteweg and de Vries in 1895 [18]. Assuming a traveling-wave solution of constant shape, one may make the ansatz η(x, t) = f (x − ct), and reduce the KdV equation to an ordinary differential equation of the form 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 realistic wave profiles. Since F (f ) is a third-order polynomial it has three roots. Periodic solutions occur only in the case where the roots are distinct, so that they can be labeled f 3 < f 2 < f 1 , and this convention allows us to write Requiring the solution to be real, it can be seen from (17) that −F (f ) must be positive. Examining the phase plane plot in Figure 3, it is clear that periodic solutions will oscillate between f 2 and f 1 .
Since f 2 < f 1 by assumption, f 1 will denote the wavecrest while f 2 will be the trough. Consequently, the wave height H is given by the difference H = f 1 − f 2 . 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 where m is the elliptic modulus defined by m = f 1 −f 2 f 1 −f 3 , and K(m) is the complete elliptic integral of first kind [24]. The wavespeed c and the wavelength λ are then given by In order to use these formulas in practice, it is convenient to take the wave height H, the mean surface level η and the elliptic parameter m as given parameters. The roots of F can then by computed explicitly and are given as follows.
The shoaling problem can be formulated as follows. Suppose we know that a periodic wave train of wave height H, wavelength λ and mean-surface level η is given at a local depth h A . Assuming that the waves can be approximately described by the the KdV equation, a unique cnoidal wave solution can be found, and the frequency ν, the energy flux q E and the radiation stress S xx be found. The shoaling problem now consists of finding an appropriate cnoidal solution at a smaller depth h B . Relying on the conservation of frequency and energy flux, and a prescribed change in the radiation stress, the shoaling equations can be written as The first equation is conservation of frequency ν, while the second equation expresses conservation of energy flux integrated over a period T [4]. Finally, the third equation is (16), and is indeed a formulation 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 1 , f 2 and f 3 . However, the method 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. This approach was therefore not able to predict the wave set-down. Moreover, solving for the parameters f 1 , f 2 and f 3 was not optimal from a numerical point of view, and the algorithm terminated long before the highest wave was reached. In the present work, we use the explicit form of the parameters to solve directly for H, m and η. This approach is more expedient from a numerical point of view, and allows us to go much higher up on the shoaling curve. In terms of the wave height H, the conservation of frequency turns into the following cubic equation for H: Similarly we can manipulate the equation describing the changes in radiation stress to find an expression for the set-down. Indeed, (23) reduces to the quadratic equation with ).
Note that we had to integrate various powers and derivatives of cn 2 (ξ; m). These formulas are based 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 η 0 (m) and then find H by solving the cubic equation. This can in turn be used to update the set-down by Repeating this process, we continue to approximate H and η until a stopping criterion has been reached. Using this reduction allows us to solve the nonlinear system (21), (22), (23), by the following procedure. From the value of m at depth h A , we increment up, at each step iterating the above equations for H and η, and then checking whether (22) is satisfied to a specified tolerance. If that is the 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.

Implementation of the shoaling equations
Having explained the nonlinear shoaling equations, and the strategy for finding approximate solutions, 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 an incoming wave of wavelength λ 0 is given at a starting depth h 0 . While this depth may not be deep water, it is assumed that h 0 > 0.28λ 0 , so that we label this starting value the deep-water wavelength. First, the linear shoaling equation is used up to the point h/λ 0 = 0.1. At this point the cnoidal theory is valid [39], and we propose a matching technique to obtain the fundamental parameters of the nonlinear wave. Lastly, we use the nonlinear KdV theory to follow the shoaling curve into the nonlinear region.

Stage 1. Linear theory:
The essential ingredient in this step is the linear dispersion relation and the conservation of energy. Since the circular frequency ω is conserved, the nonlinear equation ω 2 − gk B tanh (k B h B ) = 0 needs to be solved numerically for the wave number k B using a Newton solver. The wave height is then found from the conservation of energy as where C g = dω dk is the group velocity depending on the wave number and depth. The derivation of equation (27) is described in detail by for example in [12]. The set-down can be solved similarly using the procedure outlined in [28] using the radiation stress.
Stage 2. Matching linear and nonlinear theory: Due to increasing waveheight, the nonlinear algorithm must be utilized on the final stretch of the shoaling curve. The transition must be made in a region where both theories are valid, and some general conditions on the transition depth are given in [45]. In order to initialize the nonlinear solver at the transition depth, three parameters need to 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 these parameters and then keeping H and η constant, we find an elliptic modulus m ∈ (0, 1) which leads to equality in the chosen parameter. In principle we can only match one parameter, and it is not clear which one will be the most convenient. To investigate this issue, we define the root problems 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 to varying m, and therefore is the best candidate for a matching parameter. Figure 4 shows a particular example.
Stage 3. Cnoidal shoaling: The final step in the shoaling procedure is solving for wave height and set-down using the scheme defined in Section 6. First, define H and η as functions of m as given by formula (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 will then determine changes of wave height and set-down at consecutive points, allowing us to move up the slope. Having defined the shoaling equations and the numerical approach, we may plot the development of the wave height and set-down for a practical shoaling problem. A comparison of our numerical results with the wave tank data obtained by Bowen et al. [9] and the classical linear theory is shown in Figure 5. In this figure, the red circles represent the wave height of the shoaling wave according to the experimental data of [9], and the red asterics represent the experimentally determined set-down and set-up, i.e. the deviation of the mean depth from the still water line (S.W.L.). The dashed blue curves show the prediction of the linear theory as already presented in [9]. Finally, the solid curves represent the shoaling model put forward in the present article. The waves under consideration had an initial wave length of λ 0 = 202cm and an initial wave height H 0 = 6.45cm before they reached the toe of the slope. The beach had a slope of 1 : 12.
The linear shoaling model agrees well with the experimental data up to moderate wave heights. On the other hand, the benefit of the nonlinear formulation can be clearly seen. Indeed, the nonlinear theory clearly yields a better fit of both wave height and set-down close to the breaking point of the wave. One interesting detail in the lower panel of Figure 5 is that the mean depth appears to rise in the nonlinear theory while it continues on downwards in the linear theory. This nonlinear mid-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 further offshore was used in [17], and Cokelet's extension of Stokes' approximation of periodic waves was used in [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 we have to acknowledge that the KdV model ceases to be valid for breaking waves.

Shoaling without set-down
It was argued in [43] that the wave set-down is negligible in many practical cases. With this assumption, the shoaling problem can be simplified since only two equations need to be solved. This approach 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 mentioned, the use of a first-order approximation of the energy flux [44] gave rise to a discontinuity in wave height at the matching point between the linear and their nonlinear theory due to the inconsistency of the approximation of the energy flux. A partial fix was presented in [45] where a continuous shoaling curve is obtained by use of conversion tables but at the cost of not conserving the energy flux. More recently, the shoaling theory was extended in [20] using the nonlinear definition of the energy flux Figure 5: Profile of the mean water level η and the wave height H as functions of the depth, compared to data points taken from [9]. In this case, we have wave period T = 1.14s incident wave height H 0 = 6.45cm and breaking wave height H b = 8.55cm. In the experiment, the beach slope is tan(β) = 0.082. The lower panel depicts the same scenario in a different aspect ratio with a zoom-in around the still water line (S.W.L.). developed in [3,4]. However, due to the specific way the problem was formulated, it was not possible to 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 from Section 4, and the implementation explained in Section 5, we can find a continuous shoaling curve that 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 the wave height in the shoaling region. First formulate H = H(m) according to formula (24) given in Section 6. Then use the assumption of zero set-down to find the roots given by (20). Having the roots as 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 the shoaling region as compared to the curves presented by [20]. This is due to the simplicity of the implementation, solving two nonlinear equations rather than three. We also see in Figure 6 that our theory is in fairly good agreement with the shoaling curves presented in [44] which are in reasonable agreement with higher-order theories. Moreover, it is clearly visible in Figure 6 that our curves do not feature a discontinuity in wave height.
Finally, comparisons of shoaling curves computed using the method at hand with data provided by wave tank experiments are presented. The experimental data are taken from [45] where considerable pains were taken to stay within the correct scaling regime. The waves were generated with a pistontype wavemaker, then propagated over a flat bottom with still water depth of 36 cm before shoaling We observe in Figure 7 that there is very good agreement between the numerical model (shown in blue) and the experimental data (indicated by red crosses). The only exception is the experiment with deep water steepness H 0 /λ 0 = 0.064, shown in the uppermost panel in Figure 7. This is a rather steep wave and it can be observed that even the linear theory fails to get good agreement. This error is then propagated and compounded by the nonlinear theory. On the other hand, as long as the deepwater steepness is moderate enough, the linear theory does an adequate job, and we observe fairly good agreement in the plots between experimental data and numerical simulations for both linear and nonlinear regimes. One should note that the higher-order theory by Cokelet was used numerically for the same data set in [36], but no better agreement was found than with the present theory. We have also done some preliminary computations with a full Boussinesq model [35], and similar agreement with the wavetank data was found after some experimentation with grid size and friction coefficients.

Conclusion
In this paper we have shown how to derive expressions for the energy flux and radiation stress of a wavetrain in the context of the well known KdV equation. The derivations are based on a general approach 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 feature additional 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 small to moderate amplitude, and since such waves can be described by cnoidal functions, the shoaling equations can be posed as a 3×3 system of equations for the three parameters of the cnoidal functions. The shoaling equations are then formulated in a numerically convenient way, and it is shown that they can be solved to yield shoaling curves up to the breaking point. In particular, our formulation resolves a problem encountered in earlier work [20] where the curves terminated prematurely due to numerical instability. The curves are compared with experimental data and higher-order theories, and are found to be quite accurate. In addition, since the shoaling equations are defined using approximations of energy flux and radiation stress which are consistent with the approximations made in the KdV Figure 7: Comparisons of shoaling curves from six experiments with a range of initial wave steepness. equation itself, the shoaling curves are continuous, and no ad-hoc fixes such as advocated for in [45] are necessary to get a continuous development of the wave height, Finally, since the radiation stress is included in the formulation of the shoaling equations, it is possible to compute the set-down as part of the shoaling problem. Comparisons with experimental data from [9] show that the set-down can be computed accurately. The nonlinear set-down defined here may also potentially be useful in the study of waves interacting with currents [43].
To put our study into context, note that due to the ready availability of computational resources, shoaling is nowadays usually computed by utilizing nearshore wave models which are used to forecast wave conditions in the coastal zone. These models are generally known as Boussinesq models, and probably the first such model was developed by Peregrine [33]. These models have become fairly sophisticated in recent years, being able to treat even fairly short waves due to extending the dispersion properties using ideas of Nwogu [31] and Witting [49]. Some of the models currently in use are described in [10,30,35]. In particular, as shown for example in [29], shoaling can be computed accurately with such models. Following ideas of [38], some models have been extended to be able to treat larger-amplitude waves, and it is now possible to simulate waves with fully nonlinear and highly dispersive models [23,47] or with the full Euler equations [13]. Nevertheless, for practical purposes, these models need to be initialized with data which are generally given by stochastic wave models, and in general the transition between stochastic and deterministic models is still poorly understood. The method of computing wave height and set-down put forward here may have some potential if coupled with stochastic input data (see [14]) since the computational complexity of our method is by far smaller than for any phase-resolving nearshore model and no tuning is necessary with the current model.

A Integrals of cnoidal functions
In order to define the two shoaling models we need to handle integrals of various powers of derivatives of η. In particular we need to determine η, η 2 , η 3 , η 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 the more convenient form (for more details see [43]). Similarly considerations apply to the powers of η which involve integrals of the form These expressions can be found in [24]. The terms η xx and ηη xx also include such terms due to the identity η xx (2Kξ; m) = 4K 2 H(2 − 2m + (8m − 4)cn 2 (2Kξ; m) − 6mcn 4 (2Kξ; m)), and similar relations which can be found in [1]. Combining the results above we deduce that η = f 2 + H