Wind-Induced Changes to Shoaling Surface Gravity Wave Shape

Unforced shoaling waves experience growth and changes to wave shape. Similarly, wind-forced waves on a flat-bottom likewise experience growth/decay and changes to wave shape. However, the combined effect of shoaling and wind-forcing on wave shape, particularly relevant in the near-shore environment, has not yet been investigated theoretically. Here, we consider small-amplitude, shallow-water solitary waves propagating up a gentle, planar bathymetry forced by a weak, Jeffreys-type wind-induced surface pressure. We derive a variable-coefficient Korteweg-de Vries-Burgers (vKdV-B) equation governing the wave profile's evolution and solve it numerically using a Runge-Kutta third-order finite difference solver. The simulations run until convective prebreaking -- a Froude number limit appropriate to the order of the vKdV-B equation. Offshore winds weakly enhance the ratio of prebreaking height to depth as well as prebreaking wave slope. Onshore winds have a strong impact on narrowing the wave peak, and wind also modulates the rear shelf formed behind the wave. Furthermore, wind strongly affects the width of the prebreaking zone, with larger effects for smaller beach slopes. After converting our pressure magnitudes to physically realistic wind speeds, we observe qualitative agreement with prior laboratory and numerical experiments on breakpoint location. Additionally, our numerical results have qualitatively similar temporal wave shape to shoaling and wind-forced laboratory observations, suggesting that the vKdV-B equation captures the essential aspects of wind-induced effects on shoaling wave shape. Finally, we isolate the wind's effect by comparing the wave profiles to the unforced case. This reveals that the numerical results are approximately a superposition of a solitary wave, a shoaling-induced shelf, and a wind-induced, bound, dispersive, and decaying tail.


I. INTRODUCTION
Wind coupled to surface gravity waves leads to wave growth and decay as well as changes to wave shape.However, many aspects of wind-wave coupling are not yet fully understood.Since the sheltering theory of wind-wave coupling by Jeffreys [1], a variety of mechanisms for wind-wave interactions have been put forward, often with a focus on calculating growth rates [2,3].Furthermore, these theories have been tested by many studies in the laboratory [4][5][6][7] and the field [8,9].Similarly, numerical studies modeled the airflow above waves using methods such as large eddy simulations [10][11][12] or modeled the combined air and water domain using Reynolds-averaged Navier-Stokes (RANS) solvers [13] or direct numerical simulations [14,15].While wave growth rates and airflow structure have received much attention, wind-induced changes to wave shape have been less studied.Unforced, weakly nonlinear waves on flat bottoms (e.g., Stokes, cnoidal, and solitary waves) are horizontally symmetric about the peak (i.e., zero asymmetry), but are not vertically symmetric (i.e., nonzero skewness, e.g., [16,17]).Wave shape (skewness and asymmetry) influences many physical phenomena, such as wave asymmetry sediment transport [18,19] and extreme waves [20,21].Laboratory experiments of wind blowing over periodic waves demonstrated that wave asymmetry increases with onshore wind speed in intermediate water [22] and deep water [23].Theoretical studies likewise showed that wind-induced surface pressure induces wave shape changes in both deep [24] and shallow [25] water.However, the influence of wind on wave shape has not yet been investigated theoretically for shoaling waves on a sloping bottom.
Additionally, the shoaling of unforced waves on bathymetry also induces wave growth and shape change.Field observations revealed the importance of nonlinearity in wave shoaling and its relation to skewness and asymmetry [26,27].Additionally, laboratory experiments of waves shoaling on planar slopes yield how the wave height and wave shape evolve with distance up the beach [28][29][30].Furthermore, numerical studies investigated wave shoaling all the way to wave breaking.A variety of methods were utilized, including pseudospectral models [31], the fully nonlinear potential flow boundary element method solvers [32,33], the large eddy simulation volume of fluid methods [33], and two-phase direct numerical simulations of both the air and water [34].Theoretical [35] and numerical [33] investigations of wave breaking showed that convective wave breaking depends on the surface water velocity u and the phase speed c and occurs when the Froude number Fr := u/c is approximately unity.The type of wave breaking (e.g., spilling, plunging, surging, etc.) is related to the beach slope β, initial wave height H 0 , and initial wave width L 0 through the Iribarren number Ir := β/ H 0 /L 0 [36,37].
Few studies looked at the combined effects of wind and shoaling of surface gravity waves.Experimental studies found that onshore wind increases the surf zone width [38] and decreases the wave height-to-water depth ratio at breaking [39], with offshore wind having the opposite effect.Additionally, numerical studies using two-phase RANS solvers of wind-forced solitary [40] and periodic [41] breaking waves demonstrated that increasingly onshore winds enhance the wave height at all points prior to breaking.Regarding wave shape, only Feddersen and Veron [23] and O'Dea et al. [42] investigated the combined influence of wind and shoaling.Feddersen and Veron [23] demonstrated experimentally that onshore winds enhance the shoaling-induced time-asymmetry.In field observations of random waves, cross-shore wind was weakly correlated to the overturn aspect ratio of strongly nonlinear, plunging waves with offshore (onshore) wind reducing (increasing) the aspect ratio [42].However, the wind variation was relatively weak and covariation with other parameters was not considered.Despite a growing literature of wave shape measurements and simulations, a theoretical description of wind-induced changes to wave shoaling (e.g., wave shape, breaking location, etc.) has not yet been developed.
This study will derive a simplified, theoretical model for wind-forced shoaling waves that takes the form of a variablecoefficient Korteweg-de Vries (KdV)-Burgers equation, which is a generalization of the standard KdV equation.When the bottom bathymetry is allowed to vary, the coefficients of the KdV equation are no longer constant and the system is described by a variable-coefficient KdV (vKdV) equation [43,44].The deformation of solitary-wave KdV solutions propagating without wind on a sloping-bottom vKdV system was studied both analytically [45] and numerically [31].Shoaling causes solitary wave initial conditions to deform and gain a rear "shelf" for small enough slopes [45].Alternatively, if the flat-bottomed KdV equation is augmented with a wind-induced surface pressure forcing, the constant-coefficient KdV-Burgers (KdV-B) equation results [24].Wind in the KdV-B equation induces a solitary-wave initial condition to continuously generate a bound, dispersive, and decaying tail with polarity depending on the wind direction [24], analogous to a KdV non-solitary-wave initial condition [46].The variable-coefficient KdV-Burgers equation combines both shoaling and wind forcing into one equation.Our development of a simplified, analytic model for the coupling of shoaling and wind-forcing highlights the relative importance of these phenomena and provides a concise framework for analyzing their competing effects.
In Sec.II, we apply a wind-induced pressure forcing over a sloping bathymetry to derive a vKdV-Burgers equation and determine a convective prebreaking condition.We then solve the resulting vKdV-Burgers equation numerically using a third-order Runge-Kutta solver and investigate the changes to wave shape and prebreaking location in Sec.III.Finally, we examine the relationship between pressure and wind speed, isolate the effect of wind from the effect of shoaling, and discuss how our findings relate to previous laboratory and numerical studies in Sec.IV.

A. Governing equations
We derive a vKdV-Burgers equation for wind-forced shoaling waves by utilizing the standard [27,47] simplifying assumptions of planar, two-dimensional waves propagating in the +x-direction on incompressible, irrotational, inviscid fluid without surface tension.Additionally, we choose the +z-direction to be vertically upwards with the z = 0 at the still water level and impose a bottom bathymetry at z = −h(x).The standard incompressibility, bottom boundary, kinematic boundary, and dynamic boundary conditions are We introduce the wave profile η(x, t), the velocity potential φ(x, z, t) derived from the water velocity u = ∇φ, the surface pressure p(x, t), the gravitational acceleration g, and the water density ρ w , which is much larger than the air density ρ a ≈ 1.225 × 10 −3 ρ w .Additionally, we remove the Bernoulli constant from the dynamic boundary condition Schematic showing the (periodic) simulation domain and relevant lengthscales.The blue line represents the water surface and wave profile η and the solid black line is the bottom bathymetry h(x).The solitary wave initial condition has a height H0 and effective half-width L0 and begins with its peak on the far left side, in the middle of flat region of depth h0.
The initial wave then propagates to the right with phase speed c up the beach with slope β until it reaches prebreaking (cf.Sec.II F).The positive/negative wind speed U corresponds to an onshore/offshore wind forcing.
by using the φ gauge freedom.Next, to examine the wind's effect on shoaling waves, we impose the analytically simple Jeffreys-type surface pressure p(x, t) forcing [1] p(x, t) = P ∂η(x, t) ∂x . ( The pressure constant P ∝ ρ a (U − c) 2 depends on the wave phase speed c and wind speed U (cf. Sec.IV A).For a wave propagating towards the shore, onshore winds yield P > 0 whereas offshore winds give P < 0. The Jeffreys-type forcing is likely most relevant for near-breaking waves [48] or strongly forced steep waves [49,50] (see discussion in Zdyrski and Feddersen [25]).Recent large eddy simulations [51], two-phase direct numerical simulations [52], and RANS [41] started investigating the coupling of periodic waves and wind.Some simulations [12,51] suggested a phase shift of approximately 135°to 155°between waves and the surface pressure, in contrast to the 90°s hift predicted by the Jeffreys-type forcing.However, these studies used periodic waves and are not directly applicable to the solitary waves we consider, where an angular phase shift is undefined.Xie [40] considered wind-wave coupling for shoaling solitary waves, but pressure distributions were not reported.Therefore, we will prioritize the analytical simplicity of the Jeffreys-type forcing for our analysis.

B. Model domain and model parameters
The model domain (Fig. 1) is similar to that of Knowles and Yeh [31] and consists of an initial flat section of length L f at a depth of h 0 and transitions smoothly at x = 0 into a planar beach region with constant slope β and characteristic beach width L b := h 0 /β.The bathymetry then smoothly transitions to a flat plateau of length 2L f at a depth of h = 0.1h 0 followed by a downward slope with slope −β.Finally, there is another flat section of length L f at a depth of h 0 before the domain wraps periodically, which simplifies the boundary conditions.
The initial condition will be a KdV solitary wave with height H 0 and width L 0 following Knowles and Yeh [31], and L 0 will be specified later.The solitary wave begins centered on the left boundary, in between the two flat, deep sections of length L f .From the defined dimensional quantities, we specify four nondimensional parameters Here, ε 0 is the nondimensional initial wave height, µ 0 is the square of the nondimensional initial inverse wave width, P 0 is the nondimensional pressure magnitude (normalized by the initial wave width), and γ 0 is ratio of the initial wave width to the beach width.Note that the wave width-to-beach width ratio γ 0 is related to the beach slope β as γ 0 = β/ √ µ 0 .Together, these four nondimensional parameters control the system's dynamics.

C. Nondimensionalization
We nondimensionalize our system's variables using the characteristic scales described in Sec.II B: the initial depth h 0 ; the initial wave's height H 0 ; the initial wave's horizontal lengthscale L 0 ; the gravitational acceleration g; and the pressure magnitude P .Using primes for nondimensional variables, we normalize as Zdyrski and Feddersen [25] did and define We later assume the nondimensional parameters ε 0 , µ 0 , γ 0 , and P 0 are small to leverage a perturbative analysis.For the constant slope β beach profile, the spatial derivative of the bathymetry is also small from the different nondimensionalizations of h and x).However, perturbation analyses are simplest when all nondimensional variables are O (1).Therefore, we leverage the two, horizontal lengthscales L 0 and L b (cf.Sec.II B) to define a nondimensional, stretched bathymetry h′ that depends on x/L b = γ 0 x ′ as h′ (γ Then, denoting derivatives with respect to γ 0 x ′ using an overdot, the derivative of h is ḣ′ : and the small slope becomes explicit as Now, the nondimensional equations take the form For the remainder of Sec.II, we remove the primes for clarity.

D. Boussinesq equations, multiple-scale expansion, and vKdV-Burgers equation
We follow the conventional Boussinesq equation derivation presented in, e.g., Mei et al. [53] or Ablowitz [54].The two modifications we include are the weakly sloping bottom, similar to the treatment in Johnson [43] and Mei et al. [53], and the inclusion of a pressure forcing like that of Zdyrski and Feddersen [25].However, the joint contributions of both pressure and shoaling are new.For the sake of brevity, we only detail the relevant differences here, but we still treat the derivation formally to ensure a proper ordering of small terms and obtain the parameters' validity ranges.First, we expand the velocity potential in a Taylor series about the bottom z = −h(x) as Substituting this expansion into the incompressibility Eq. ( 8) and bottom boundary condition (9) and assuming µ 0 ≪ 1 gives φ as a function of the velocity potential evaluated at the bottom ϕ := φ 0 .If we further assume that the bottom is very weakly sloping γ 0 ∼ µ 0 ≪ 1, this simplifies to Note that the assumption γ 0 ∼ µ 0 ≪ 1 implies a moderate slope and is used by several other authors [31,43,45].For reference, if µ 0 = γ 0 = 0.1, then this implies a physically realistic β = 0.03.
Substituting this φ expansion (13) into the kinematic boundary condition (10) and dynamic boundary condition (11) yields Boussinesq-type equations with a pressure-forcing term Note that replacing h with the total depth h total = h + ε 0 η shows that these are equivalent to the flat-bottomed Boussinesq equations with h total = 1+ε 0 η.In other words, any sloping-bottom terms ḣ only appear in the combination ∂ x h total = γ 0 ḣ + ε 0 ∂ x η.This is expected since the only sloping-bottom term µ 0 γ 0 ḣ∂ x φ in the governing equations ( 8), ( 9), (10), and (11) was dropped when we neglected terms of O µ 2 0 , γ 2 0 , γ 0 µ 0 .Since the bathymetry varies on the slow scale x/γ 0 , we expand our system in multiple spatial scales x n = γ n 0 x for n = 0, 1, 2, . .., so the derivatives become and the bathymetry is a function of the long spatial scale h = h(x 1 ).Then, we expand η and ϕ in asymptotic series of ε 0 Similar to Johnson [43], we replace x 0 and t with left-and right-moving coordinates translating with speed c(x 1 ) dependent on the stretched coordinate x 1 : Then, we replace the derivatives ∂ t and ∂ x0 with Now, we will assume that ε 0 ∼ P 0 ∼ µ 0 ≪ 1 and follow the standard multiple-scale technique [53,54].The order-one terms O ε 0 0 from Eqs. ( 14) and (15) yield wave equations for φ 0 and η 0 with right-moving solutions propagating with the slowly varying, linear shallow-water phase speed c(x 1 ) = h(x 1 ).Continuing to O(ε 0 ) of the asymptotic expansion gives Eliminating η 1 from these equations gives The left-hand operator ∂ 2 /∂ξ − ∂ξ + is the same as the O(1) differential operator (20).Therefore, the right-hand side must vanish to prevent φ 1 from developing secular terms.Thus, the right-hand side becomes the variable-coefficient Korteweg-de Vries-Burgers (vKdV-Burgers) equation Finally, multiplying Eq. ( 25) by ε 0 , adding the O(1) differential equation ∂ ξ− η 0 = 0 derived from Eq. ( 20), and transforming back to the original, nondimensional variables x and t yields This vKdV-Burgers equation represents an analytically simple system for studying the effects of wind-forcing and shoaling.In the absence of wind-forcing (P 0 = 0), it reduces to the variable-coefficient KdV equation [46], while it simplifies to the constant-coefficient KdV-Burgers equation in the case of flat bathymetry (c = 1) [25].The formal derivation revealed the depth-dependent c factor in the wind-forced P 0 term which was absent in the flat-bottomed case.The pressure term P 0 ∂ 2 x η 0 functions as a damping, positive viscosity for offshore P 0 < 0 wind, making Eq.( 26) a (forward) vKdV-Burgers equation.Conversely, onshore P 0 > 0 wind causes a growth-inducing, negative viscosity giving the backward vKdV-Burgers equation.Though the backward, constant-coefficient KdV-Burgers equation is ill posed in the sense of Hadamard [55], this is irrelevant here owing to the finite time the wave takes to reach the beach.

E. Initial conditions
Our initial condition will be the solitary-wave solutions of the unforced (P 0 = 0), flat-bottom KdV equation.These waves balance the KdV equation's nonlinearity η 0 ∂ x η 0 and dispersion ∂ 3 x η 0 terms, propagate without changing shape, and require that the height H 0 and width L 0 satisfy the constraint H 0 L 2 0 = constant.Therefore, we now fix the previously unspecified L 0 by choosing µ 0 = (3/4)ε 0 so L 0 acts like an effective half-width for the solitary wave initial condition [53] While the unforced KdV equation also possesses periodic solutions called cnoidal waves, we only consider solitary waves here.

F. Convective breaking criterion
The asymptotic assumptions used to derive the vKdV-Burgers equation ( 26) fail when the wave gets too large.Therefore, we require a condition to determine when the simulations should stop.Brun and Kalisch [35] defined a convective breaking condition for solitary waves on a flat-bottom depending on the surface water velocity u s (x, t) and the phase speed c using the local Froude number Fr := ε 0 u s (x, t)/c(t), with the ε 0 coming from nondimensionalization.Convective breaking occurs wherever max x (Fr) = 1, where max x represents the maximum over x.However, when the Froude number approaches the breaking value of unity, our weakly nonlinear asymptotic assumption used to derive the vKdV-Burgers equation is violated.Thus, we instead stop our simulations at the smaller prebreaking Froude number Fr pb := 1/3 and define the "prebreaking" time t pb as the first time such that max x (Fr) = Fr pb := 1/3.Likewise, we define x pb as the location where Fr = Fr pb , which will be very near the wave peak.
As the solitary wave propagates on a slope, the wave evolves over time and the phase speed c can be ambiguous.One option is to use the adiabatic approximation derived by the authors of [45] for unforced solitary waves on very gentle slopes with x peak the location of the wave peak.Alternatively, Derakhti et al. [33] used large eddy simulations to numerically investigate unforced solitary wave breaking on slopes ranging from β = 0.2-0.005for two different forms of c.They found wave breaking at max x (Fr) = 0.85 when using the speed of the numerically tracked wave peak c peak .However, they also found that the shallow-water approximation c shallow = h(x peak ) + ε 0 η(x peak ) [equivalent to c adi to O ε 2 0 ] was within 15% of c peak near breaking.Therefore, we will use Eq. ( 28) for c(t) owing to its simplicity and theoretical foundation.Finally, though these studies all considered unforced solitary waves, our results will show that c adi varies approximately 3% across pressure magnitudes P 0 for our simulations, so this approximation is valid.
We still need an expression for the wave velocity at the surface u s (x, t), which we now derive by modifying the example of Brun and Kalisch [35] to include sloping bathymetry and pressure forcing.Combining the vKdV-Burgers equation ( 26) and kinematic boundary condition (14) Here, ċ := ∂ x1 c(x 1 ) = O(1) in analogy to the previously defined ḣ.Assuming an ansatz we insert Eqs. ( 30) and ( 31) into Eq.( 29), drop terms of O ε 2 0 and solve for A, B, C, and D by using the independence of ε 0 , γ 0 , µ 0 , and P 0 : Note that A represents the nonlinear contribution, B the effect of shoaling, C the dispersive effect, and D the pressure forcing.Finally, the Taylor expansion of φ(x, z) Eq. ( 13) gives the fluid velocity at the surface u s (x, t) using u = ∂ x φ evaluated at z = ε 0 η: Therefore, the Froude number is calculated as with u s (x, t) given by Eq. (33).Equation (34) demonstrates that the choice of Fr pb = 1/3 keeps the problem in the correct asymptotic regime.Our simulated waves will reach prebreaking in depths of h(x pb ) = 0.59-0.72,so Eqs.( 33) and (34) show that Fr pb = 1/3 corresponds to a dimensional H(x pb )/h(x pb ) ≈ Fr pb / h(x pb ) = 0.39-0.43.
Considering an asymptotic regime around ε 0 of ε 0 , or 0.089-0.45,we see that H(x)/h(x) remains in the asymptotic regime.

G. Numerics
The vKdV-Burgers equation ( 26) lacks analytic, solitary-wave-type solutions, so we solve it numerically using a third-order explicit Runge-Kutta adaptive time stepper with the error controlled by a second-order Runge-Kutta method as implemented in SciPy [56].The domain's deep, flat regions have length L f = 20L 0 .Therefore, the spatial domain's entire, periodic width L x varies between 108 and 150, depending on the slope β.We discretize the spatial domain using a fourth-order finite difference method with N x = 2 2 ⌊L x /0.025/2 2 ⌉ points which yields a grid spacing of dx = L x /N x ≈ 0.025 and ensures that the grid can be refined by two steps.We employ adaptive time stepping to keep the relative error below 10 −3 and the absolute error below 10 −6 at each step.For all cases, the average time step is ∆t ≈ 2 × 10 −4 .The pressure is initially turned off until the solitary wave is one unit (i.e., a half-width L 0 ) away from the start of the beach slope.The pressure is linearly ramped up to its full value over the time the wave takes to cross a full-width 2L 0 .We included a hyperviscosity ν 4 ∂ 4 x η 0 with ν 4 = 1 × 10 −5 for numerical stability [57].Finally, we only considered solitary waves, as mentioned previously, because cnoidal wave trains require spin-up and damping since we only consider prebreaking waves.
We validated the solver against the unforced, flat-bottom analytical solution and had a normalized root-mean-square error of 1.6 × 10 −4 after nondimensional time t = 50 (the longest simulation time) as well as a normalized wave height change of To verify the numerical convergence on our nominal (fine) grid size, we also ran the simulation on a medium grid and a fine grid, each with a refinement ratio of 2. We then calculated the Romberg interpolation [58] and compared each grid's normalized root-mean-square error with respect to this interpolation to yield a grid refinement index (GCI) [59] for the unforced case with the steepest β = 0.025 and shallowest β = 0.01 slopes.Both slopes had a convergence order p > 2.1 which yielded grid convergence indices indices GCI nominal,medium < 3 × 10 −4 and GCI medium,coarse < 8 × 10 −4 much less than unity and implying that the results are grid-converged.Furthermore, the profiles for all times until prebreaking were nearly identical to the simulations of Knowles and Yeh [31] for an unforced solitary wave shoaling on a slope.Finally, the simulation reproduced the finding of Knowles and Yeh [31] that small waves (ε 0 ≪ 1) on weak slopes (γ 0 ≪ 1) yield Green's Law for the wave height H(x) := max t (η) ∝ h(x) 1/4 (with max t the maximum over time t), while moderate waves (ε 0 < 1) on very weak slopes (γ 0 ≪ 1) give Miles' adiabatic law H(x) ∝ h(x) −1 [60].
The vKdV-Burgers equation ( 25) is determined by two nondimensional parameter combinations: the pressure term P 0 /ε 0 and the shoaling term γ 0 /ε 0 .Recall that the dispersive term µ 0 /ε 0 is a redundancy which we fixed by specifying L 0 (cf.Sec.II E).We investigate this two-dimensional parameter space by choosing ε 0 = 0.2 and µ 0 = 0.15 and varying the beach slope β = 0.01-0.025and pressure |P | = 0 to 0.05 (cf.Sec.IV A for a discussion of the size of P ).This yields a total of 36 simulations (Table I).Note that Eq. ( 25) demonstrates changing ε 0 → λε 0 is equivalent to γ 0 → γ 0 /λ in the wave's comoving reference frame.Therefore, solutions for waves with different initial heights ε 0 can be generated from our solutions to the vKdV-Burgers equation in the laboratory frame Eq. ( 26) by scaling the height, boosting, and adjusting γ 0 .Note the asymptotic expansion assumed P 0 ∼ ε 0 , or P/(ρ w gL 0 ε 0 ) ∼ 1, but the pressure values we are using (Table I) are smaller than unity.Nevertheless, multiple-scale expansions are often accurate outside their parameters' validity ranges and this constraint would be satisfied asymptotically for smaller values of ε 0 .

H. Shape statistics
While previous studies investigated shoaling's effect on wave-breaking location, we will theoretically examine the combined influences of both wind and shoaling on wave location x pb at prebreaking t pb (Sec.II F).To estimate how x pb changes, we first calculate the shoreline x shore as the location where the bathymetry would intersect z = 0 if it had a constant slope β without our shallow plateau (Fig. 2).Then, we calculate the prebreaking zone width as L pz := x pb − x shore .For a given beach slope β, we will analyze the change in the prebreaking zone width relative to the unforced case ∆L pz := L pz − L pz P =0 normalized by the unforced prebreaking zone width L pz P =0 .This bulk statistic ∆L pz /( L pz P =0 ) determines the variance in prebreaking locations as a fractional change of the prebreaking zone width.
Additionally, since we wish to analyze the effect of wind and shoaling on wave shape, we will investigate four shape statistics that vary as the wave propagates.The first three are local shape parameters defined at each location x.First, we directly examine the maximum Froude number max t (Fr) expressed in Eq. (34).Second, we investigate the maximum height relative to the local water depth max t (η)/h(x) at each location x.Third, we consider the maximum slope max t (|∂η/∂x|).Both the relative height and maximum slope contribute to the convective breaking criterion max x (Fr) = Fr pb .Finally, we introduce a bulk shape parameter, the full width of the wave at half of the wave's maximum (FWHM) L W (t) normalized by the local water depth h(x).For our unforced KdV solitary wave initial condition (27), the FWHM divided by the initial depth is Schematic showing the definition of the prebreaking zone and shore locations.The blue line represents the water surface and wave profile η at prebreaking and the solid black line is the bottom bathymetry h(x).The bathymetry consists of a flat region of depth h0, a sloping region, and a shallow plateau.The shoreline x shore (black dot) is the location where the bathymetry would intersect the still water level if it had a constant slope (dashed line).The beach width L b is the distance from the toe of the beach slope to x shore and the prebreaking point x pb is the location on the wave where Fr = Fr pb , which will be very near the wave peak.The prebreaking zone width Lpz is the distance from x pb to x shore .
shape parameter defined at each point in time t with the local parameters defined at each point in space.Therefore, we define L W (x) = L W [t peak (x)] at the time t peak (x) when the wave peak passes location x.

III. RESULTS
Now, we use the results of the numerical simulations to investigate the effect of wind on solitary wave shoaling and shape.We will present shape statistics (Sec.II H) for the 20 different runs (Table I) to detail the wave-shape changes and prebreaking behavior across the parameter space.For the remainder of the paper we will return to dimensional variables.

A. Profiles of shoaling solitary waves with wind
First, we qualitatively investigate the effect of varying pressures P and bathymetric slopes β on solitary-wave shoaling by examining the wave profile η/h 0 , normalized by the initial depth h 0 , at three different times t (Fig. 3) corresponding to when the solitary wave first feels the slope (t = 0), the time of prebreaking (t = t pb ), and half-way between (t = t pb /2).Note that these t = 0 wave profiles (purple in Fig. 3) are nearly identical to the sech 2 (x/L 0 ) initial condition (27) since the waves have only propagated over a flat bottom [Figs.3(g) and 3(h)] and the pressure has not yet been turned on.Halfway to prebreaking (t = t pb /2, blue), the solitary wave has grown through shoaling with a steeper front face (+x side) and increased asymmetry for all P and β.At the time of prebreaking (t = t pb , green) the solitary wave has increased in height, steepened, and gained a substantial rear shelf for all P and β.These changes are likely even larger for waves propagating to full wave breaking (Fr ≈ 1).The generation of rear shelves by shoaling solitary waves as in Fig. 3 was first calculated by Miles [45] and resulted from the mass shed by the sech 2 wave as it narrowed.Onshore wind (P > 0) reinforces the shoaling-based wave growth and yields relatively narrow peak widths for both β [Figs.3(a) and 3(b)].In contrast, offshore wind (P < 0) reduces the wave shoaling, but results in wider peak widths [Figs.3(e) and 3(f)].For instance, the width at prebreaking on the milder slope β = 0.015 is L W /h(x) = 3.75 under onshore winds and L W /h(x) = 4.26 under offshore winds.These differences in wave-shoaling result in the offshore-forced (P < 0) solitary wave reaching prebreaking (x pb , ×'s in Fig. 3) farther onshore (shallower water) than the onshore-forced (P > 0) solitary wave.For the β = 0.015 case, the onshore-forced wave reaches prebreaking at x pb /h 0 = 20.8 while the offshore-forced wave reaches prebreaking farther offshore at x pb /h 0 = 26.1.Similarly, the larger beach slope [β = 0.025, Figs.3(b), 3(d), and 3(f)] causes waves to reach x pb in less horizontal distance, though they prebreak in shallow water than the milder beach slope (β = 0.015) waves.For the unforced case, the β = 0.015 wave prebreaks at x pb /h 0 = 23.3 while the steeper beach β = 0.025 causes prebreaking at x pb /h 0 = 15.2.At t = t pb , the rear shelf is wider and extends higher up the rear face for offshore winds [≈ 0.1h 0 in Fig. 3(e)] than for onshore winds [≈ 0.07h 0 in Fig. 3(a)].Again, we expect these differences to be The profile times shown depend on the Froude number (34) and therefore vary between the panels.The first profile (purple) occurs when the peak is located at x = −L0 where the pressure begins turning on and the time is defined so t = 0 here.The last profile (green) occurs when the convective prebreaking condition maxx(Fr) = Fr pb = 1/3 is met (cf.Sec.II F), and the middle profile (blue) occurs at a time halfway between the first and last profiles.Both columns have ε0 = 0.2 and µ0 = 0.  (34), and the ×'s on the last profiles (green) are the prebreaking locations x pb .We only display a subset of the full spatial domain.Note that the aspect ratio is chosen to highlight the wind's effect on the shoaling solitons.
even larger for fully breaking waves (Fr ≈ 1).As the control case, the unforced (P = 0) solitary wave has x pb located between the onshore and offshore wind cases with an intermediate rear shelf.Finally, the milder slope (β = 0.015) has a sharper, more pronounced rear shelf while the steeper slope (β = 0.025) has a more gently sloping rear shelf.
We next investigate the impact of onshore [Figs.4( (28) at each location multiplied by the prebreaking Froude number Fr pb = 1/3.The ×'s denote the locations with the highest Froude number and the ×'s on the last (green) profiles are the prebreaking locations x pb .The squares are the locations of the maximum slope magnitude |∂η/∂x| and the upside-down triangles represent the locations of the maximum wave velocity profile.We only display a subset of the full spatial domain.Note that the aspect ratio is chosen to highlight the wind's effect on the shoaling solitons.
forcing's sign.The wave slope [Figs.4(c) and 4(d)] highlights the shoaling-and wind-induced shape changes by accentuating the front-rear asymmetry.At t = 0 [purple in Figs.4(a) and 4(b)], the wave slope has odd-parity about the peak.However, as the wave propagates onshore, both the front and rear faces steepen, though the front face steepens more dramatically.The influence of the wind is most noticeable in three aspects: the offshore-forced wave [P = −0.05,Fig. 4(b)] is 10% smaller than the onshore forced wave [P = 0.05, Fig. 4(a)]; the offshore-forced rear-face wave slope [Fig.4(d)] is 16% smaller than the onshore-forced wave slope [Fig.4(c)], though the front-face slope is only 2% smaller; and the trailing shelf's slope extends further behind the offshore-forced wave [≈ 8h 0 , Fig. 4(d)] than the onshore-forced wave [≈ 5h 0 , Fig. 4(c)].The wave velocity profile u/ √ gh 0 [Eq.(33), Figs.4(e) and 4(f)] nearly mirrors the wave profile [Figs.4(a) and 4(b)], as is expected given that u ∝ η to leading order (33).Finally, the phase speed c adi [red, Eq. ( 28)] decreases as the wave shoals which enhances convective prebreaking, though c adi only varies 3% between the onshore and offshore winds.Note, in Figs.4(e) and 4(f), c adi is multiplied by Fr pb = 1/3 so that the intersection of the red curve with the wave velocity profile occurs at x pb , the location of prebreaking.We highlight that the prebreaking quantity u/c ≈ Fr = 1/3 is smaller than the unified breaking onset criteria u/c = 0.85 described by Derakhti et al. [33].

B. Shape statistics with shoaling and variations of prebreaking zone width with wind
Building on the previous qualitative descriptions of the wave shape, slope, and velocity, we now quantify the change in the shoaling wave's shape parameters for onshore and offshore P (Fig. 5).First, we consider the maximum Froude number max t (Fr) as a function of nondimensional position x/h 0 [Fig.5(a)].In the flat region (x < 0), the maximum Froude number is max t (Fr) = 0.1986, and it increases as the waves shoal to the prebreaking value max t (Fr) = Fr pb = 1/3 (light gray line).The wind has a significant impact on the location of prebreaking x pb , with onshore wind (red) causing the Froude number to increase faster and x pb to occur farther offshore than offshore wind (blue) does.This can also be seen in Figs.4(e) and 4(f), where the maximum velocities u/ √ gh 0 (upside-down triangles), which are proportional to max t (Fr), are growing faster for the onshore wind [Fig.4(e)] than the offshore wind [Fig.4(f)].Notably, at a fixed location x/h 0 , the max x (Fr) varies substantially (e.g., 0.28 to 0.32 at x/h 0 = 20).In addition, we consider the maximum height max t (η) at a fixed location and normalized by the local water depth h(x) [Fig.5(b)].For all pressures P , the solitary wave increases in height, but the onshore wind enhances this growth while the offshore wind partially suppresses the growth.Again, this is apparent in the evolution of the maximums η(x peak )/h 0 in Fig. 3, with the peak locations x peak closely approximated by the ×'s marking the location of maximum Fr.Since Fr ∝ η to leading order, the relative height at prebreaking is approximately 0.41 for all P [Fig.5(b)] with offshore-forced waves slightly larger (1%) than onshore-forced waves., the steepness is enhanced by onshore wind (P > 0), suppressed for offshore wind (P < 0), and approaches nearly the same prebreaking value of 0.15 for all wind speeds, being only 1% larger for onshore winds than offshore winds.The maximum slope at prebreaking is nearly constant because solitary waves have a fixed relationship between the wave height and wave width (and hence slope) as discussed in Sec.II E. And the height at prebreaking is approximately constant since Fr ∝ η to leading order and Fr = Fr pb is constant.However, this relationship is only approximate on a slope, with deviations due to nonlinearity, dispersion, shoaling, and wind forcing (32).Finally, we examine the FWHM L W , normalized by the local water depth h(x) [Fig.5(d)].While L W /h(x) ultimately decreases from its initial value of 4.56 for all pressure magnitudes, there is significant variation in the prebreaking value.For our parameters, L W /h(x) changes 215% more for onshore wind (P = 0.05) than offshore wind (P = −0.05)from start to prebreaking.Figures 4(a) and 4(b) show that the rear shelf does not rise to half the wave height, so the FWHM does not incorporate the shelf's width.Instead, the onshore-forced narrowing is occurring in the top region above the shelf.Hence, while the relative height and slope at prebreaking are largely similar for all the wind speeds, the FWHM at prebreaking is strongly affected by the wind speed indicating wind effects on shoaling shape.We expect the wind-induced changes to maximum wave slope and FWHM to be even more stronger approaching wave breaking (Fr ≈ 1).
We also investigate the change in the prebreaking zone width ∆L pz (Sec.II H) as a function of pressure P/(ρ w gL 0 ε 0 ) for four different values of the beach slope β (Fig. 6).First, ∆L pz is linearly related to the pressure magnitude, and the wind has a larger effect on ∆L pz for smaller beach slopes, with P/(ρ w gL 0 ε 0 ) = −0.05changing the prebreaking zone width by approximately 5% for the smallest slope β = 0.01.This is because the wind has more time to affect the wave before it reaches prebreaking.This wind-induced change in prebreaking location is visible in Fig. 3, where the prebreaking location x pb (×'s on green profiles) occurs closer to the shoreline (+x direction) for offshore winds P < 0 [Figs.3(e) and 3(f)] than for onshore winds P > 0 [Figs.3(a) and 3(b)].Additionally, we note that for the smallest slope β = 0.01, the fractional change in prebreaking zone width ∆L pz /( L pz P =0 ) is asymmetric with respect to pressure, with offshore P/(ρ w gL 0 ε 0 ) = −0.05yielding a 22% larger change than onshore P/(ρ w gL 0 ε 0 ) = 0.05 (Fig. 6).For fully breaking waves, the wind-induced changes to the breaking-zone width ∆L bz would likely be larger and the unforced surf zone width L bz P =0 would be smaller (as waves propagate closer to shore before fully breaking), making the fractional change ∆L bz /( L bz P =0 ) much larger.

C. Normalized prebreaking wave shape changes induced by wind and shoaling
As Fig. 5 quantified the shape statistics at prebreaking for all x, we now directly investigate the effect of pressure P and shoaling β on prebreaking wave shape by normalizing each prebreaking wave profile η by its maximum height max x (η) and aligning the prebreaking locations x pb /h 0 (Fig. 7).Each solution is dominated by the sech 2 wave centered near x − x pb = 0, which becomes taller and narrower as the wave shoals as required by energy conservation [45].Furthermore, while the sech 2 component is symmetric in time at a fixed location, it becomes slightly deformed when viewed at a fixed time as the front face moves slower than the rear face [46,61].We also observe a shelf behind the wave, which Miles [45] calculated by requiring that the right-moving mass-flux be conserved as the sech 2 narrows and sheds mass.While long-duration calculations of the Miles shelf reveal a nearly horizontal shelf extending far behind the wave [61,62], our shelf instead slopes gently downward, likely due to insufficient development time and distance.
In Fig. 7, we plot the prebreaking wave shape for fixed bottom slope β [Fig.7(a)] and fixed pressure magnitude P [Figs.7(b) and 7(c)].For a fixed slope [Fig.7(a)], the front wave faces at prebreaking are qualitatively very similar and match an unforced solitary wave of the same height.However, wind strongly affects the rear shelves as observed in Fig. 3.The offshore winds (blue) cause the shelf to be thicker and extend higher up the rear wave face than the offshore wind (reds) do, although the shelf intersects z = 0 at (x − x pb )/h 0 ≈ −10 for all wind speeds.Nevertheless, all of the changes in Fig. 7 would likely be enhanced for fully breaking waves as wind and shoaling effects have longer to act on the wave.
We also consider the wave shape at breaking for different values of the beach slope β with a fixed onshore [Fig.7(b)] or offshore [Fig.7(c)] wind.The rear half of the wave shows that bottom slope β affects the rear shelf differently than pressure P/(ρ w gL 0 ε 0 ) does.While the shelf intersected z = 0 at the same location for all wind speeds [Fig.7(a)], increasing β causes the intersection point (i.e., the base of the shelf) to move forward and closer to the peak.Finally, the offshore wind [Fig.7(c)] causes a noticeably larger shelf than the onshore wind [Fig.7(b)] for the weakest slope β = 0.01 (purple), with a similar pattern observed in Fig. 4(a) (β = 0.015) compared to Fig. 4(b) (β = 0.025).However, this difference is much smaller for the steeper (green) slopes, implying that stronger shoaling partially suppresses the wind-induced shape change because there is less time for pressure to act prior to prebreaking.

A. Wind speed
Our derivation in Sec.II coupled wind to the wave's motion through the use of a surface pressure Eq. ( 5).The resulting vKdV-Burgers equation ( 26) has a wind-induced term dependent on the pressure magnitude constant P/(ρ w gL 0 ε 0 ).We analyze the evolution and prebreaking of solitary waves for variable P (Sec.III).While the usage of P was the most natural since it is the physical coupling between wind and waves (in the absence of viscous tangential stress), measuring the surface pressure is challenging in field observations or laboratory experiments [7,9].Therefore, we also consider the evolution and prebreaking of the shoaling solitary waves as a function of the wind speed U .Zdyrski and Feddersen [25]  to our Eq.( 27)] with dimensional form with nondimensional height ε = H/h and width L = 2h/ √ 3ε in water of depth h.Zdyrski and Feddersen [25] then related wind speed U to the surface pressure magnitude using where U is measured at a height of half the solitary wave's width.This parametrization originated from observations by Donelan et al. [9] of nonseparated wind forcing for periodic, shallow-water waves and was adapted to solitary waves in shallow water by Zdyrski and Feddersen [25].The relationship is accurate for the wind and wave conditions of Donelan et al. [9], but should be interpreted qualitatively here.Note that the radicand differs by a factor of 2 from Zdyrski and Feddersen [25] owing to the different definitions of ε.Even though Eq. ( 36) was originally applied to flatbottomed KdV solitary waves (27), our assumption that γ = L/L b ≪ 1 implies that the bathymetry is approximately flat over the wave's width 2L.Therefore, we use Eq.(36) to translate between the pressure P/(ρ w gL 0 ε 0 ) and the TABLE II: Wind speeds as functions of pressure P/(ρ w gLε) and local depth h for solitary waves (35) with ε = 0.2.U onshore corresponds to P > 0 and U offshore to P < 0. The conversion from P/(ρ w gLε) to U is given in Eq. (36).wind speed U at any point on the sloping bathymetry by using the local ε and h and relating the initial pressure to the local pressure P/(ρ w gLε) = (ε 0 L 0 /εL)P/(ρ w gL 0 ε 0 ).
Table II shows the onshore (P > 0) and offshore (P < 0) wind speeds corresponding to the pressures used in our simulations for two representative depths h.It shows that the pressure magnitudes in our simulations correspond to physically reasonable wind speeds, with onshore U from 3.1 m s −1 to 13 m s −1 for water 1 m deep or 4.9 m s −1 to 20 m s −1 for water 2.5 m deep.Notice that unforced waves with P = 0 correspond to a wind speed matching the wave phase speed U = c, with c approximately the linear shallow-water phase speed c ≈ √ gh.In particular, this means that onshore P > 0 and offshore P < 0 winds with the same pressure magnitude |P | will have different wind-speed magnitudes |U |.Additionally, note that keeping P fixed implies that the wind speed U changes as the wave shoals.This is mostly due to the decrease in the phase speed c ∝ √ gh, with higher-order effects coming from the ε and L dependence of the radicand in Eq. (36).Finally, note that as the wave shoals and ε increases, the height at which the wind speed is measured z = L/2 = h/ √ 3ε decreases.
We now reexamine our results regarding the prebreaking zone width (Fig. 6) in terms of the wind speed U/ gh(x) using Eq.(36).In addition to changing the abscissa of the plot (Fig. 8), we also modify the definition of the change in the prebreaking zone width ∆L pz := L pz − L pz U=0 by comparing and normalizing each prebreaking zone width to the U = 0 case rather than the P = 0 case.This transformation changes the initially straight lines of Fig. 6 into approximate pairs of upward-and downward-facing ∆L pz curves shifted to the right by one unit (Fig. 8).Furthermore, we see that ∆L pz is now much flatter for onshore winds (U > 0) than for equal magnitude offshore winds (U < 0).This is due to the inflection point of the unforced case (P = 0) being shifted to the right at U/ √ gh = 1.

B.
Isolating the effect of wind For no wind (P = 0), solitary wave shoaling is well understood to generate a rear shelf [45].The variation in the rear shelf's thickness with P (Fig. 7) is reminiscent of the variability in the wind-generated bound, dispersive, and decaying tails of flat-bottom solitary waves [25].Additionally, Zdyrski and Feddersen [25] showed that flatbottom, wind-generated tails are analogous to the dispersive tails of KdV solutions with non-solitary-wave initial conditions [53].Both the rear shelf and wind-generated tail can be viewed as weak perturbations to the KdV equation by transforming the nondimensional vKdV-Burgers equation ( 25) into a constant-coefficient, perturbed KdV equation by defining ν := (3/2)(ε 0 /µ 0 )η 0 / h2 and τ := cdx 1 µ 0 /(6γ 0 ): The first term on the right-hand-side (RHS) is the shoaling term which leads to the rear shelf [45], and the second term is the wind-induced Burgers term [25].Our derivation assumed all terms in Eq. ( 37) were the same order, and indeed ∂ τ h/ h ∼ 6γ 0 /µ 0 = 1.0-2.6 is order one.However, the forcing term |P 0 |/µ 0 = (4/3)|P 0 |/ε 0 = 0-0.07 is much smaller than unity and is a weak perturbation to the sloping-bottom KdV equation with its sech 2 solitary-wave and shelf solution.
Perturbed KdV equations similar to Eq. ( 37) received some attention in past literature.In particular, the unforced, shoaling case can be recast in a number of asymptotically equivalent forms including Eqs. ( 25) or (37) with P = 0. Previous authors applied different mathematical techniques to these various asymptotic forms to derive closed-form approximate solutions.For instance, Knickerbocker and Newell [61] solved the unforced analog of Eq. ( 25) using conservation laws.Alternatively, Grimshaw [63] used a multiple-scale analysis to solve a windless, modified form of Eq. ( 25).Similarly, Ablowitz and Segur [64] solved the P = 0 version of Eq. ( 37) using a direct perturbation solution while Newell [46] solved it using the inverse scattering transform.Since Eq. (37) shows that shoaling and wind forcing can be treated on an equal footing, it should be possible to extend these analyses to wind-forced solitary waves.However, such an analysis is outside the scope of the current work, but is suitable for future work.
Reverting back to dimensional variables, we isolate the effect of wind by separating out the sech 2 solitary wave and the Miles rear shelf using the unforced P = 0 normalized profiles to represent shoaling and rear-shelf generation.We define a normalized tail ζ as the difference between the forced and unforced P = 0 normalized profiles of Fig. 7: For constant depth, the height H and width L of unforced, sech 2 solitary waves always satisfy HL 2 = const.Since our numerical results (e.g., Fig. 7) are dominated by the sech 2 solitary wave profile, scaling the wave profile by H requires that we scale the spatial coordinate by L ∝ 1/ √ H to respect this symmetry and enable comparison of waves with different heights.We replace h → h(x peak ) in the expression for the flat-bottom solitary wave width L Eq. (35) to yield the wave width for a slowly varying depth as We normalize the spatial coordinate as x/L to compare the normalized tails ζ in Fig. 9.We show the normalized tail ζ versus (x − x pb )/L for different pressures P and bottom slopes β in Fig. 9. First, increasing the pressure magnitude |P | increases the tail's amplitude and wavelength.For example, the wavelength with β = 0.01 is approximately 5L for P/(ρ w gL 0 ε 0 ) = −0.025and 7.5L for P/(ρ w gL 0 ε 0 ) = −0.05.This amplitude increase is expected, as higher pressures put more energy into the tail, causing growth.Additionally, increasing the bottom slope β decreases the shelf's width and the tail's amplitude without noticeably changing its wavelength.We can explain the narrower shelf and smaller amplitude by recognizing that larger β's cause the wave to reach prebreaking (when these profiles are compared) earlier, decreasing the time over which the wind (tail) and shoaling (shelf) act.The wavelength's independence of the beach slope β also implies that the width L of the solitary wave sets the tail's wavelength.Additionally, we note that onshore (P > 0) and offshore (P < 0) winds change the polarity of the tail, consistent with Zdyrski and Feddersen [25].Lastly, wind induces a small, bound wave in front of the prebreaking solitary wave with minimum near (x − x pb ) = 0 and extremum near (x − x pb )/L ≈ 2 of the same polarity as the rear shelf (Fig. 9), similar to the flat-bottom results of Zdyrski and Feddersen [25].
Hence, the numerically calculated wave profiles (Fig. 7) are a superposition of the sech 2 solitary wave, Miles' shelf [45], and a wind-induced bound, dispersive, and decaying tail [25].Furthermore, this decomposition of the full wave enables us to understand the effects of wind and shoaling from previous studies.The sech 2 solitary wave grows and narrows due to wave shoaling [45] and wind forcing [25].Miles' shelf is generated by the mass flux of the growing wave.The shelf's absence from the normalized tails in Fig. 9 implies its shape is largely unchanged by the wind, and its amplitude for a given bottom slope β is approximately proportional to the sech 2 solitary wave.Finally, the amplitude and wavelength of the bound, dispersive, and decaying tail grow with the sech 2 solitary wave [25].This decomposition relies on the assumption that the tail and shelf are both small compared to the solitary wave and do not influence each other or the solitary wave.This is only possible when the wind-forcing P 0 is weak and the wave width-to-beach width ratio γ 0 is small.Miles [45] analyzed a vKdV equation requiring the same weak-slope assumption γ 0 ∼ ε 0 , though his adiabatic results required an even smaller γ 0 = O 10 −2 [31].The realistic beach widths we utilized yield a γ 0 = 3 × 10 −2 -6 × 10 −2 somewhat larger than this adiabatic regime, and the γ 0 /ε 0 term in Eq. ( 37) is not as small as the pressure-forcing term, implying some nonlinear interactions between the shoalinginduced shelf and the sech 2 solitary are possible.For this reason, we subtracted off the unforced solitary wave and shelf rather than approximate them analytically.Nevertheless, the pressure forcing |P 0 |/µ 0 = 0-0.07we used was sufficiently small that the weak wind forcing can be considered to interact linearly, as seen in the clean separation between pressure-induced tail and sech 2 plus shelf in Fig. 9. Thus, we show that there are physically reasonable parameter regimes wherein wave shape and evolution are the superpositions of the previously understood shoalingand wind-induced effects.

C. Breakpoint location comparison to prior laboratory experiments and models
Numerical studies investigated the effect of wind on the breaking of shoaling solitary [40] and periodic [41] waves using a RANS k-ε model to simulate both the air and water.Xie [40] considered solitary waves with initial height H 0 /h 0 = 0.28 on a beach slope of 0.05 with onshore winds of up to U/ √ gh 0 = 3, while Xie [41] investigated periodic waves with initial height H 0 /h 0 = 0.3 and initial inverse wavelength h 0 /λ 0 = 0.1 on a beach slope of 0.03 forced by onshore winds up to U/ √ gh 0 = 2.These studies determined that the (absolute) maximum wave heights max t (η)/h 0 increased with increasing onshore wind at each location x < x b prior to breaking at x b , consistent with our findings in Fig. 5(b).Furthermore, Xie [40] found the effect of wind on breaker depth is significant while the effect on breaker height is minor, again consistent with our prebreaking findings.Finally, comparing wave profiles in Xie [40] shows that onshore winds increase the wave slope at a fixed location, which is consistent with our Fig.5(c).
For periodic waves, previous laboratory experiments also investigated wind's effect on the breaking characteristics of shoaling waves [38,39].Douglass [38] considered waves with initial height H 0 /h 0 = 0.3 and initial inverse wavelength h 0 /λ 0 = 0.1 under wind speeds of up to U/ √ gh 0 = ±2.3 on a beach with slope 0.04 while King and Baker [39] considered waves with initial height H 0 /h 0 = 0.2 and initial inverse wavelength h 0 /λ 0 = 0.3 with wind speeds of up to U/ √ gh 0 = ±1.1 on a beach with slope 0.05.Douglass [38] measured how wind speed affects wave parameters and changes the surf zone width for periodic waves.Directly comparing our Fig. 8 to Fig. 2 of Douglass [38], we see many qualitative similarities, including the prebreaking zone width's flatter response near U = 0 and a stronger response for offshore winds (U < 0) than the corresponding onshore winds (U > 0), with our change roughly four times smaller than theirs.Furthermore, the laboratory studies also found that the relative breaking height H b /h b , normalized by the breaking depth, decreased by as much as 40% for offshore wind speeds of U/ √ gh b = 4 and increased by up to 10% for onshore wind speeds of U/ √ gh b = −2 compared to the unforced case [38,39].By comparison, over those same wind speed ranges of U/ √ gx pb = 1 ± 3, our simulations found that the relative prebreaking height h pb /x pb varied by approximately 1% between onshore and offshore winds [Fig.5(b)], with the same polarity as the laboratory experiments.Finally, Douglass [38] observed only a slight dependence of the breaking wave height on wind speed, which is consistent with our finding that offshore-forced waves are only 1% larger than onshore-forced waves at prebreaking.On the contrary, King and Baker [39] measured no statistically significant change with wind speed in the ratio between the change in the fractional change (H b − H b U=0 )/ H b U=0 of (absolute) breaking height H b compared to the unforced case H b U=0 .This is surprising, as both our results and those of Douglass [38] had near-constant relative (pre)breaking height ∆(H b /h b ).Indeed, constant ∆(H b /h b ) implies ∆H b / H b U=0 must vary with wind (we find approximately ±5% at prebreaking) to counteract the varying of ∆h b / h b U=0 = ∆L b / L b U=0 with wind that Douglass [38], King and Baker [39], and we all find.
Our results qualitatively agree with prior numerical results on solitary waves [40] as well as experimental and numerical results on periodic waves [38,39,41], and the quantitative mismatch can be partly explained by the different nondimensional parameters.All three studies mentioned used larger initial waves (ε 0 ≈ 0.3), so nonlinear effects were likely more important.They also used steeper beach slopes, enhancing the shoaling effect.Additionally, while the surf zone width change is roughly four times larger for Douglass [38] than for our simulations over the same windspeed range, Douglass [38] investigated waves that were actually breaking.In contrast, we stopped our simulations at prebreaking max x (Fr) = Fr pb = 1/3, significantly before actual breaking max x (Fr) ≈ 1 [33], thus we expect smaller changes to the fractional surf prebreaking zone width (cf.Sec.III B).

D. Wave shape comparison to prior laboratory experiments
Feddersen and Veron [23] experimentally examined the effect of wind on the temporal shape of shoaling, periodic waves at a fixed location for no-wind (U/c = 0) and onshore-wind (U/c = 6) cases [cf.Fig. 10(a)].The waves began shoaling at a depth of h 0 = 0.37 m and were observed at a depth h * /h 0 = 0.62.To enable direct comparison with our results, we normalized the profile by the height of the no-wind profile H * u at depth h * .We also normalized the time coordinate with an unforced solitary wave's temporal width τ * = L * u / √ gh * at depth h * , where the unforced spatial width L * u is related to the wave height via an equation similar to Eq. ( 39).We observe that the onshore wind increases the wave height and causes earlier peak arrival, relative to the arrival time of η/H * u = 0.25 (Fig. 10).For comparison, our numerical results [Fig.10(b)] over different wind speeds, calculated using Eq.(36), are shown at x/h = 20.1 corresponding to h * /h 0 = 0.68.This corresponds to x pb for the strongest onshore wind U/ √ gh * = 4 (i.e., P = 0.05) case.We truncated the time series at Fr = 1/2, instead of Fr pb = 1/3, to enable qualitative comparison of the strong-wind numerical case (U/ √ gh * = 4) with the laboratory onshore-wind case (U/c = 6).However, we note that this result is strictly out of the asymptotic range.
Our numerical results [Fig.10(b)] are qualitatively similar to the experimental results of Feddersen and Veron [23] [Fig.10(a)], despite different experimental and model conditions.Specifically, Feddersen and Veron [23] presented shallower (h * /h 0 = 0.62) measurements of shoaling periodic waves (h 0 /λ 0 = 0.27) forced by stronger laboratory winds (U/c = 6) in contrast to our deeper (h * /h 0 = 0.68) simulations of solitary waves (h 0 /λ 0 → 0) for weaker winds (U/ √ gh * = 4).Therefore, we display the results side by side as opposed to overlaid to emphasize the qualitative comparison.Onshore wind causes wave growth and narrowing for both laboratory observations and numerical results.The peaks of the onshore-forced waves occur earlier (relative to the time of η/H * u = 0.25) due to wave narrowing, as seen in Fig. 5(d).Finally, the rear faces of the onshore-forced waves (∆t/τ * from 1 to 2.5) dip below the no-wind cases.Numerical results with offshore wind continue the pattern.The qualitative similarities between our results and that of Feddersen and Veron [23] suggests that the our vKdV-Burgers equation captures the essential aspects of wind-induced effects on shoaling wave shape.No other theoretical model has yet shown such qualitative similarity with wind-forced wave shape experiments, despite the differences (e.g., periodic versus solitary) between the laboratory and numerical studies.

V. CONCLUSION
While shoaling-induced changes to wave shape are well understood, the interaction of wind-induced and shoalinginduced shape changes has been less studied and lacked a theoretical framework.Utilizing a Jeffreys-type windinduced surface pressure, we defined four nondimensional parameters that controlled our system: the initial wave height ε 0 , the inverse wavelength squared µ 0 , the pressure strength P 0 , and the wave width-to-beach width ratio γ 0 .We leveraged these small parameters to reduce the forced, variable-bathymetry Boussinesq equations to a variable-coefficient Korteweg-de Vries-Burgers equation for the wave profile η.We also extended the convective breaking criterion of Brun and Kalisch [35] to include pressure and shoaling.A third-order Runge-Kutta solver determined the time evolution of a solitary wave initial condition up a planar beach under the influence of onshore and offshore winds.
Stopping the simulations at a prebreaking Froude number of 1/3 revealed that the prebreaking relative height and maximum slope are largely independent of wind speed, but onshore winds cause a narrowing of the waves.The width of the prebreaking zone is strongly modulated by wind speed, with offshore wind decreasing the prebreaking zone width by approximately 5% for the mildest beach slopes.Investigating the wave shape at prebreaking revealed that the front of the wave is relatively unchanged and matches an unforced solitary wave, while the rear shelf is strongly affected by wind speed and bottom slope.We isolated the effect of wind from the effect of shoaling and revealed a bound, dispersive, and decaying tail similar to wind-induced tails on flat bottoms.By leveraging the relationship between surface pressure P and wind speed U , we directly compared our results to existing experimental and numerical results.We found qualitative agreement in surf width changes and wave height changes and expect better quantitative agreement as the waves propagate closer to breaking.These results suggest that wind significantly impacts wave breaking, and our simplified model highlights the relevant physics and changes to wave shape.Future avenues of research could include calculating asymptotic, closed form solutions to Eq. ( 37) or deriving coupled equations for both the water and air motions to more accurately predict the surface pressure distribution.

FIG. 3 .
FIG.3.Shoaling solitary-wave η evolution under [(a),(b)] onshore P > 0, [(c),(d)] unforced P = 0, and [(e),(f)] offshore P < 0 wind-induced surface pressure versus nondimensional distance x/h0 as the wave propagates up the [(g),(h)] planar bathymetry.The profile times shown depend on the Froude number(34) and therefore vary between the panels.The first profile (purple) occurs when the peak is located at x = −L0 where the pressure begins turning on and the time is defined so t = 0 here.The last profile (green) occurs when the convective prebreaking condition maxx(Fr) = Fr pb = 1/3 is met (cf.Sec.II F), and the middle profile (blue) occurs at a time halfway between the first and last profiles.Both columns have ε0 = 0.2 and µ0 = 0.15, and [(a),(e)] the left-column forced cases have |P/(ρwgL0ε0)| = 0.05 and β = 0.015 while [(b),(f)] the right-column forced cases use |P/(ρwgL0ε0)| = 0.025 and β = 0.025.The ×'s denote the locations with the highest Froude number(34), and the ×'s on the last profiles (green) are the prebreaking locations x pb .We only display a subset of the full spatial domain.Note that the aspect ratio is chosen to highlight the wind's effect on the shoaling solitons.
FIG.3.Shoaling solitary-wave η evolution under [(a),(b)] onshore P > 0, [(c),(d)] unforced P = 0, and [(e),(f)] offshore P < 0 wind-induced surface pressure versus nondimensional distance x/h0 as the wave propagates up the [(g),(h)] planar bathymetry.The profile times shown depend on the Froude number(34) and therefore vary between the panels.The first profile (purple) occurs when the peak is located at x = −L0 where the pressure begins turning on and the time is defined so t = 0 here.The last profile (green) occurs when the convective prebreaking condition maxx(Fr) = Fr pb = 1/3 is met (cf.Sec.II F), and the middle profile (blue) occurs at a time halfway between the first and last profiles.Both columns have ε0 = 0.2 and µ0 = 0.15, and [(a),(e)] the left-column forced cases have |P/(ρwgL0ε0)| = 0.05 and β = 0.015 while [(b),(f)] the right-column forced cases use |P/(ρwgL0ε0)| = 0.025 and β = 0.025.The ×'s denote the locations with the highest Froude number(34), and the ×'s on the last profiles (green) are the prebreaking locations x pb .We only display a subset of the full spatial domain.Note that the aspect ratio is chosen to highlight the wind's effect on the shoaling solitons.

9 FIG. 4 .
FIG. 4. Shoaling solitary-wave [(a),(b)] nondimensional profile η/h0, [(c),(d)] slope ∂η/∂x, and [(e),(f)] nondimensional wave velocity profile u/ √ gh0 under [(a),(c),(e)] onshore and [(b),(d),(f)] offshore wind-induced surface pressure as the wave propagates up the [(g),(h)] planar bathymetry.All panels use the same [(g),(h)] bathymetry and differ only in the sign of the pressure forcing.Values are shown versus nondimensional distance x/h0 for ε0 = 0.2, µ0 = 0.15, |P/(ρwgL0ε0)| = 0.05, β = 0.015, and nondimensional times t √ gh/L0 indicated in the legends.The red lines in [(e),(f)] represent the phase speed c adi(28) at each location multiplied by the prebreaking Froude number Fr pb = 1/3.The ×'s denote the locations with the highest Froude number and the ×'s on the last (green) profiles are the prebreaking locations x pb .The squares are the locations of the maximum slope magnitude |∂η/∂x| and the upside-down triangles represent the locations of the maximum wave velocity profile.We only display a subset of the full spatial domain.Note that the aspect ratio is chosen to highlight the wind's effect on the shoaling solitons.

Figure 5 (
Figure 5(c) shows the evolution of the maximum wave slope magnitude max t |∂ x η|, corresponding to the front face's slope [Figs.4(c) and 4(d)].Like the relative height [Fig.5(b)], the steepness is enhanced by onshore wind (P > 0), suppressed for offshore wind (P < 0), and approaches nearly the same prebreaking value of 0.15 for all wind speeds, being only 1% larger for onshore winds than offshore winds.The maximum slope at prebreaking is nearly constant because solitary waves have a fixed relationship between the wave height and wave width (and hence slope) as discussed in Sec.II E. And the height at prebreaking is approximately constant since Fr ∝ η to leading order and Fr = Fr pb is constant.However, this relationship is only approximate on a slope, with deviations due to nonlinearity, dispersion, shoaling, and wind forcing(32).Finally, we examine the FWHM L W , normalized by the local water depth h(x) [Fig.5(d)].While L W /h(x) ultimately decreases from its initial value of 4.56 for all pressure magnitudes, there is significant variation in the prebreaking value.For our parameters, L W /h(x) changes 215% more for onshore wind (P = 0.05) than offshore wind (P = −0.05)from start to prebreaking.Figures4(a) and 4(b) show that the rear shelf does not rise to half the wave height, so the FWHM does not incorporate the shelf's width.Instead, the onshore-forced narrowing is occurring in the top region above the shelf.Hence, while the relative height and slope at prebreaking are largely similar for all the wind speeds, the FWHM at prebreaking is strongly affected by the wind speed indicating wind effects on shoaling shape.We expect the wind-induced changes to maximum wave slope and FWHM to be even more stronger approaching wave breaking (Fr ≈ 1).

FIG. 5 .
FIG. 5. Shoaling solitary-wave shape statistics under onshore and offshore pressure forcing versus nondimensional distance x/h0.The (a) Froude number maxt(Fr) (34), (b) maximum height normalized by the local water depth maxt(η)/h(x), (c) maximum slope maxt(|∂η/∂x|), and (d) full width at half maximum normalized by the local water depth LW /h(x) (cf.Sec.II H) are displayed at each location along the (g) planar bathymetry.Results are shown for ε0 = 0.2, µ0 = 0.15, β = 0.015, and pressure magnitude |P/(ρwgL0ε0)| up to 0.05, as indicated in the legend.The solid black line is the unforced case, P = 0.The light gray line in (a) represents the convective prebreaking Froude number Fr pb = 1/3 at which the simulations were stopped.

FIG. 8 .
FIG. 8.The fractional change in prebreaking zone width ∆Lpz compared to the unforced case Lpz U =0 (cf.Sec.II H) versus the nondimensional wind speed U/ gh(x pb ) normalized by the local, shallow-water phase speed gh(x pb ) and evaluated at a height of half the solitary wave width L. The results are shown for beach slopes β = 0.01-0.025.

FIG. 10 .
FIG. 10.Nondimensional wave profile η/H * u at a fixed location versus nondimensional time ∆t/τ * with η normalized by the no-wind wave height H * u .The time difference ∆t := t − t0.25 is relative to the time when the profile reaches η/H * u = 0.25 and is normalized by the unforced solitary wave temporal width τ * = L * u / √ gh * at a depth h * , with L * u the unforced spatial width [Eq.(39)].(a) Periodic wave laboratory results from Feddersen and Veron[23] with ε0 ≈ 0.28, µ0 = 0.035, and β = 0.125.The solid line represents the no-wind case U/c = 0 and the dashed line corresponds to an onshore wind U/c = 6, as indicated in the legend.The waves began shoaling at a depth of h0 = 0.37 m and were measured at a nondimensional depth h * /h0 = 0.62.The wave height η is plotted relative to the wave trough mint(η).(b) Numerical results for ε0 = 0.2, µ0 = 0.15, and β = 0.015 at a nondimensional water depth h * /h0 = 0.69 and for nondimensional wind speeds U/ √ gh * [cf.Eq. (36)] as indicated in the legend.

TABLE I :
Range of nondimensional parameters simulated.