1. Introduction
The study of wind and ocean wave interactions began with Jeffreys (Reference Jeffreys1925) and continues to be an active field of research (e.g. Janssen Reference Janssen1991; Donelan et al. Reference Donelan, Babanin, Young and Banner2006; Sullivan & McWilliams Reference Sullivan and McWilliams2010). Many theoretical studies (e.g. Jeffreys Reference Jeffreys1925; Miles Reference Miles1957; Phillips Reference Phillips1957) focus on calculating wind-induced growth rates and often employ phase-averaging techniques. However, experimental (e.g. Leykin et al. Reference Leykin, Donelan, Mellen and McLaughlin1995; Feddersen & Veron Reference Feddersen and Veron2005) and theoretical (e.g. Zdyrski & Feddersen Reference Zdyrski and Feddersen2020) studies have shown wind can also influence wave shape, quantified by third-order shape statistics such as skewness and asymmetry, corresponding to vertical and horizontal asymmetry, respectively. Furthermore, while many numerical studies on coupled wind and waves employ sinusoidal water waves and therefore neglect wind-induced shape changes (e.g. Hara & Sullivan Reference Hara and Sullivan2015; Husain et al. Reference Husain, Hara, Buckley, Yousefi, Veron and Sullivan2019), some recent numerical studies have incorporated wind-induced changes to the wave field using coupled air–water simulations (e.g. Liu et al. Reference Liu, Yang, Guo and Shen2010; Hao & Shen Reference Hao and Shen2019) or direct numerical simulations of two-fluid flows (e.g. Zonta, Soldati & Onorato Reference Zonta, Soldati and Onorato2015; Yang, Deng & Shen Reference Yang, Deng and Shen2018). Wave shape influences sediment transport, affecting beach morphodynamics (e.g. Drake & Calantoni Reference Drake and Calantoni2001; Hoefel & Elgar Reference Hoefel and Elgar2003), while wave skewness affects radar altimetry signals (e.g. Hayne Reference Hayne1980) and asymmetry influences ship responses to wave impacts (e.g. Soares, Fonseca & Pascoal Reference Soares, Fonseca and Pascoal2008).
Waves in shallow water, where 
$kh \ll 1$ (with 
$h$ the water depth, 
$k=2{\rm \pi} /\lambda$ the wavenumber and 
$\lambda$ the wavelength), differ qualitatively from those in intermediate (
$kh \sim 1$) to deep (
$kh \gg 1$) water. For waves with small amplitudes 
$a_0 \ll h$, leveraging the small parameters 
$a_0/h \sim (kh)^2 \ll 1$ yields the Boussinesq equations with weak dispersion and nonlinearity. When dispersion balances nonlinear focusing, a special class of waves, known as solitary waves, are formed and appear in environments ranging from nonlinear optical pulses (e.g. Kivshar Reference Kivshar1993) to astrophysical dusty plasmas (e.g. Sahu & Tribeche Reference Sahu and Tribeche2012). These well-understood waves are often used to study fluid dynamical (e.g. Munk Reference Munk1949; Hammack & Segur Reference Hammack and Segur1974; Miles Reference Miles1979; Lin & Liu Reference Lin and Liu1998) and engineering (e.g. Monaghan & Kos Reference Monaghan and Kos1999; Lin Reference Lin2004; Xu, Chen & Chen Reference Xu, Chen and Chen2018) contexts owing to their simplicity. One of the simplest equations displaying solitary waves is the Korteweg–de Vries (KdV) equation, which incorporates dispersion and nonlinearity. When augmented with a dissipative term, this becomes the KdV–Burgers equation, with applications to damped internal tides (e.g. Sandstrom & Oakey Reference Sandstrom and Oakey1995), electron waves in graphene (e.g. Zdyrski & McGreevy Reference Zdyrski and McGreevy2019) and viscous flow in blood vessels (e.g. Antar & Demiray Reference Antar and Demiray1999). While field observations (e.g. Cavaleri & Rizzoli Reference Cavaleri and Rizzoli1981) have investigated the wind-induced growth of shallow-water waves, the interaction of wind and shallow-water waves has not yet been formulated into a simple equation such as the KdV–Burgers equation.
The influence of wind on wave shape has been previously investigated in intermediate and deep water (Zdyrski & Feddersen Reference Zdyrski and Feddersen2020). However, the coupling between wind and wave shape has not yet been investigated in shallow water. To investigate wind and surface wave interactions in shallow water over a flat bottom, we introduce a wind-induced pressure term to the Boussinesq equations in § 2. The resulting KdV–Burgers equation governs a solitary wave's evolution, which we solve numerically to yield the wave's energy, skewness and asymmetry in § 3. We calculate the wind speed, discuss the asymmetry and compare our results to intermediate- and deep-water waves in § 4.
2. Derivation of the KdV–Burgers equation
2.1. Governing equations
We treat the flow as irrotational and inviscid and neglect surface tension. Furthermore, we restrict ourselves to planar wave propagation in the 
$+x$ direction. Finally, we choose a coordinate system with 
$z=0$ at the mean water level and a horizontal, flat bottom located at 
$z=-h$. Then, the incompressibility condition and standard boundary conditions are
\begin{gather} 0 = \phi_{xx} + \phi_{zz} \quad \text{on}\ -h < z < \eta , \end{gather}
\begin{gather}0 = \phi_{z} \quad \text{on}\ z={-}h , \end{gather}
\begin{gather}\phi_{z} = \eta_{t} + \phi_{x} \eta_{x} \quad \text{on}\ z = \eta , \end{gather}
\begin{gather}0 = \frac{p}{\rho_w} + g\eta + \phi_{t} + \frac{1}{2} \left[\phi_{x}^2 + \phi_{z}^2\right] \quad \text{on}\ z= \eta . \end{gather}
Here, 
$\eta (x,t)$ is the wave profile, 
$\phi (x,z,t)$ is the flow's velocity potential related to the velocity 
$\boldsymbol{u} = {\boldsymbol{\phi}}$, 
$p(x,t)$ is the surface pressure, 
$g$ is the gravitational acceleration and 
$\rho _w$ is the water density. We used the 
$\phi$ gauge freedom to absorb the Bernoulli ‘constant’ 
$C(t)$ in the dynamic boundary condition. We seek a solitary, progressive wave which decays at infinity, 
$\eta (\boldsymbol {x},t) \to 0$ as 
$\left |\boldsymbol {x}\right | \to \infty$, with similar conditions on 
$\boldsymbol {u}$. We choose a coordinate system where the average bottom horizontal velocity vanishes,
\begin{equation} \overline{\frac{\partial \phi}{\partial x}} = 0 \quad \text{on}\ z={-}h , \end{equation}
with the overline a spatial average 
$\bar {f} := \lim _{L\to \infty } \int _{-L}^{L} f \,\mathrm {d}x / (2L)$. Additionally, we assume the surface pressure 
$p(x,t)$ is a Jeffreys-type forcing (Jeffreys Reference Jeffreys1925),
\begin{equation} p(x,t) = P \frac{\partial \eta(x,t)}{\partial x} . \end{equation}
Here, 
$P$ is proportional to 
$(U-c)^2$, with 
$c$ the wave's nonlinear phase speed and 
$U$ the wind speed (cf. § 4.1). Note that 
$P>0$ corresponds to (‘onshore’) wind in the same direction as the wave while 
$P<0$ denotes (‘offshore’) wind opposite the wave. We use a Jeffreys forcing for its analytic simplicity and clear demonstration of wind–wave coupling. Jeffrey's separated sheltering mechanism is likely only relevant in special situations (e.g. near breaking, Banner & Melville Reference Banner and Melville1976, or for steep waves under strong winds, Touboul & Kharif Reference Touboul and Kharif2006; Tian & Choi Reference Tian and Choi2013). Additionally, numerical simulations of sinusoidal waves suggest the peak surface pressure is shifted approximately 
$135^{\circ }$ from the wave peak, while Jeffreys would give a 
$90^{\circ }$ shift (Husain et al. Reference Husain, Hara, Buckley, Yousefi, Veron and Sullivan2019). However, a fully dynamic coupling between wind and waves – necessary for an accurate surface pressure over a non-sinusoidal, dynamic water surface – is outside the scope of this paper. Furthermore, the applicability of Jeffreys forcing to extreme waves means our theory could apply to the wind forcing of rogue waves in shallow water (Kharif et al. Reference Kharif, Giovanangeli, Touboul, Grare and Pelinovsky2008).
2.2. Non-dimensionalization
We non-dimensionalize our system with the known characteristic scales: the horizontal length scale 
$L$ over which 
$\eta$ changes rapidly, expressed as an effective wavenumber 
$k_E := 2 {\rm \pi}/L$; the (initial) wave amplitude 
$a_0 = H_0/2$ (i.e. half the wave height 
$H_0$); the depth 
$h$; the gravitational acceleration 
$g$; and the wind speed 
$U$, expressed as a pressure magnitude 
$P \propto \rho _a (U-c)^2$, with 
$\rho _a \approx 1.225\times 10^{-3} \rho _w$ the density of air. Denoting non-dimensional variables with primes, we have
\begin{equation} \left. \begin{gathered} x = \frac{x'}{k_E} = h \frac{x'}{\sqrt{\mu_E}}, \quad t = \frac{t'}{k_E c_0} = \frac{t'}{\sqrt{\mu_E}} \sqrt{\frac{h}{g}} , \quad \eta = a_0 \eta' = h \varepsilon \eta' ,\\ z = h z' ,\quad P = \varepsilon P' \frac{\rho_w g}{k_E} = \frac{\varepsilon}{\sqrt{\mu_E}} P' \rho_w c_0^2 ,\quad \phi = \phi'\frac{a_0}{k_E}\sqrt{\frac{g}{h}} = \frac{\phi'\varepsilon}{\sqrt{\mu_E}} c_0 h , \end{gathered} \right\} \end{equation}
with linear, shallow-water phase speed 
$c_0 = \sqrt {g h}$. Our system's dynamics is controlled by three small, non-dimensional parameters: 
$\varepsilon := a_0/h$, 
$\mu _E := (k_E h)^2$ and 
$P k_E/(\rho _w g)$. We will later require 
${O}(\varepsilon ) = {O}(\mu _E) = {O}(Pk_E/(\rho _w g))$. Now, our non-dimensional equations take the form
\begin{gather} 0 = \mu_E \phi'_{x'x'} + \phi'_{z'z'} \quad \text{on}\ {-}1 < z' < \varepsilon \eta' , \end{gather}
\begin{gather}0 = \phi'_{z'} \quad \text{on}\ z'={-}1 , \end{gather}
\begin{gather}\phi'_{z'} = \mu_E \eta'_{t'} + \varepsilon \mu_E \phi'_{x'} \eta'_{x'} \quad \text{on}\ z' = \varepsilon \eta' , \end{gather}
\begin{gather}0 = \varepsilon P' \eta'_{x'} + \eta' + \phi'_{t'} + \frac{1}{2} \left(\varepsilon \phi_{x'}^{\prime \, 2} + \frac{\varepsilon}{\mu_E} \phi_{z'}^{\prime \, 2}\right) \quad \text{on}\ z'= \varepsilon \eta' . \end{gather}We will drop the primes throughout the remainder of this section for readability.
2.3. Boussinesq equations, multiple-scale expansion, KdV equation and initial condition
Here, we modify the Boussinesq equation's derivation provided by Mei, Stiassnie & Yue (Reference Mei, Stiassnie and Yue2005) or Ablowitz (Reference Ablowitz2011) by including the surface pressure forcing in (2.4). Taylor expanding the velocity potential 
$\phi$ about the bottom, 
$z=-1$, and applying Laplace's equation (2.8) and the bottom boundary condition (2.9) yields an expansion of 
$\phi$ in terms of 
$\mu _E \ll 1$ and the velocity potential at the bottom, 
$\varphi := \left .\phi \right |_{z=-1}$. This 
$\phi$ expansion can be substituted into the two remaining boundary equations, (2.10) and (2.11), to give the Boussinesq equations with a pressure forcing term,
\begin{gather} \partial_t \eta + \partial_x^2 \varphi + \varepsilon \partial_x \left(\eta \partial_x \varphi\right) -\tfrac{1}{6}\mu_E \partial^4_x \varphi = {O}(\mu_E^2) , \end{gather}
\begin{gather}\partial_t \varphi + \varepsilon P \partial_x \eta + \eta - \tfrac{1}{2}\mu_E \partial_t \partial_x^2 \varphi + \tfrac{1}{2}\varepsilon\left(\partial_x \varphi\right)^2 = {O}(\mu_E^2) . \end{gather}
Further, we will now assume 
${O}(\varepsilon ) = {O}(\mu _E) \ll 1$.
We now expand 
$t$ using multiple time scales 
$t_n = \varepsilon ^n t$ for 
$n= 0,1$, so all time derivatives become 
$\partial _t \to \partial _{t_0} + \varepsilon \partial _{t_1}$. Then, we write 
$\eta$ and 
$\varphi$ as asymptotic series of 
$\varepsilon$,
\begin{equation} \eta(x,t) = \sum_{k=0}^{\infty} \varepsilon^k \eta_{k}(x,t_0,t_1) \quad \text{and}\quad \varphi(x,t) = \sum_{k=0}^{\infty} \varepsilon^k \varphi_{k}(x,t_0,t_1,) . \end{equation}
Now, we will reduce the Boussinesq equations, (2.12) and (2.13), to the KdV equation following a similar method to Mei et al. (Reference Mei, Stiassnie and Yue2005) and Ablowitz (Reference Ablowitz2011). Collecting order-one terms 
${O}(\varepsilon ^0)$ from (2.12) and (2.13) gives a wave equation for 
$\eta _0$ and 
$\phi _0$. The right-moving solutions are
\begin{equation} \varphi_0 = f_0(x-t_0,t_1) \quad \text{and}\quad \eta_0 = f'_0(x-t_0,1) \quad \text{with}\ f_0' := \left.\frac{\partial f_0(\theta,t_1)}{\partial \theta}\right|_{\theta = x-t_0} . \end{equation}
Continuing to the next order of perturbation theory, we retain terms of 
${O}(\varepsilon )$,
\begin{gather} \frac{\partial \eta_1}{\partial t_0} + \frac{\partial^2 \varphi_{1}}{\partial x^2} ={-}\frac{\partial \eta_0}{\partial t_1} - \frac{\partial}{\partial x} \left(\eta_0\frac{\partial \varphi_0}{\partial x}\right) + \frac{1}{6} \frac{\mu_E}{\varepsilon} \frac{\partial^4 \varphi_0}{\partial x^4} , \end{gather}
\begin{gather}\eta_1 + \frac{\partial \varphi_1}{\partial t_0} ={-}P \frac{\partial \eta_0}{\partial x} -\frac{\partial \varphi_0}{\partial t_1} + \frac{1}{2} \frac{\mu_E}{\varepsilon} \frac{\partial^3 \varphi_0}{\partial t_0 \partial^2 x} - \frac{1}{2} \left(\frac{\partial \varphi_0}{\partial x}\right)^2. \end{gather}
Inserting our leading-order solutions for 
$\eta _0$ and 
$\varphi _0$, eliminating 
$\eta _1$ and preventing resonant forcing of 
$\varphi _1$ gives the KdV–Burgers equation,
\begin{equation} \frac{\partial \eta_0}{\partial t_1} + \frac{3}{2} \eta_0 \frac{\partial \eta_0}{\partial x} + \frac{1}{6} \frac{\mu_E}{\varepsilon} \frac{\partial^3 \eta_0}{\partial x^3} ={-}P \frac{1}{2} \frac{\partial^2 \eta_0}{\partial x^2} . \end{equation}
Note that (2.18) has a rescaling symmetry, with 
$\mu _E \to \lambda ^2 \mu _E$ equivalent to taking 
$(x,t_0,t_1,P) \to (x,t_0,t_1,P)/\lambda$. Therefore, we fix the length scale (equivalently, 
$k_E$) by choosing 
$\mu _E = 6 \varepsilon$. Note that incorporating slowly varying bottom bathymetry 
$\partial _x h = {O}(\varepsilon )$ can yield an equation of the form (2.18) with spatially varying coefficients (e.g. Johnson Reference Johnson1972; Ono Reference Ono1972), although such an analysis is outside the scope of this study.
For offshore wind, the pressure term 
$P \partial ^2_x \eta _0$ acts as a positive viscosity causing damping, and (2.18) is the (forward) KdV–Burgers equation with 
$P<0$. However, for onshore wind, the viscosity is negative and causes wave growth, yielding the backward KdV–Burgers equation with 
$P>0$. The backward KdV–Burgers equation is ill posed in the sense of Hadamard because the solution is highly sensitive to changes in the initial condition (Hadamard Reference Hadamard1902). While a finite-time singularity (i.e. wave breaking) is likely, the multiple-scale expansion used to derive (2.18) is only valid for time intervals of 
${O}(1/\varepsilon )$, and we limit our analysis to short times removing the need to regularize the solution.
The solitary-wave solutions of the unforced (
$P=0$) KdV equation exist due to a balance of dispersion 
$\partial _x^3 \eta _0$ with focusing nonlinearity 
$\eta _0 \partial _x \eta _0$ and have the form (e.g. Mei et al. Reference Mei, Stiassnie and Yue2005)
\begin{equation} \eta_0 = H_0 \,\mathrm{sech}^2 \left(\frac{x}{\varDelta}\right) \quad \text{with}\ \varDelta = \sqrt{\frac{8}{H_0}} , \end{equation}
in a co-moving frame with 
$H_0>0$ an order-one parameter. For reference, unforced solitary waves travel relative to the laboratory frame with non-dimensional, nonlinear phase speed (e.g. Mei et al. Reference Mei, Stiassnie and Yue2005)
\begin{equation} c = 1 + \varepsilon \frac{H_0}{2}. \end{equation}
We use (2.19) for our initial condition and choose 
$H_0 = 2$ so the initial, dimensional amplitude 
$a_0$ is half the wave height (cf. § 2.2). Note that the unforced KdV equation also has periodic solutions known as cnoidal waves. For a fixed height, these cnoidal waves have a smaller characteristic wavelength 
$1/k_E$ than solitary waves and can be studied by choosing larger 
$\mu _E > 6 \varepsilon$ (cf. § 4.3). However, wind-induced shape changes are more readily understood when considering solitary waves owing to their reduced number of free parameters (i.e. 
$\mu _E$). Furthermore, since solitary waves are well understood and highly relevant to fluid dynamical systems (e.g. Hammack & Segur Reference Hammack and Segur1974; Miles Reference Miles1979; Lin & Liu Reference Lin and Liu1998; Monaghan & Kos Reference Monaghan and Kos1999), we will restrict our analysis to solitary waves for brevity and clarity. The wind-forcing term 
$P \partial _x^2 \eta _0$ in (2.18) disrupts the solitary wave's balance of dispersion and nonlinearity, inducing growth/decay and shape changes. The KdV–Burgers equation has no known solitary-wave solutions, so we will solve it numerically.
2.4. Numerics and shape statistics
To solve (2.18) numerically, we will use the Dedalus spectral solver (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020) which implements a generalized tau method with a Chebyshev basis. Since the onshore wind, 
$P>0$ case is ill posed, we require an implicit solver, so time stepping is done with coupled four-stage, third-order diagonally implicit Runge–Kutta and explicit Runge–Kutta schemes. The spatial domain has a length of 
$L=80$, and we require 
$\eta _0 = 0$ at 
$x' = -40$ and 
$\eta _0 = \partial _x \eta _0 = 0$ at 
$x' = 40$. We employ 
$N_c = 1600$ Chebyshev coefficients and zero padding with a scaling factor of 
$3/2$ to prevent aliasing of nonlinear terms. This corresponds to 
$N_x = 2400$ spatial points with spacing 
$\Delta x = 7.7\times 10^{-5}\text { to }7.9\times 10^{-2}$ for an average spacing of 
$\Delta x = 0.05$. The simulation runs from 
$t_1 = 0$ to 
$t_1 = T = 10$, since the multiple-scale expansion of § 2.3 is only accurate for times of 
${O}(1/\varepsilon )$. Adaptive time stepping is employed such that the Courant–Friedrichs–Lewy number is 
$(\Delta t) \max (\eta _0)/(\Delta x) = 1$. For the unforced case, this corresponds to 
$\Delta t \approx 7.86\times 10^{-3}$, increases to 
$1.04\times 10^{-2}$ for 
$P=-0.25$ and decreases to 
$4.73\times 10^{-4}$ for 
$P = 0.25$. We found that linearly ramping up 
$P$ from 
$0$ at 
$t_1=0$ to its full value at 
$t_1 = \varepsilon$, or full, dimensional time 
$T_0 = 1/(\sqrt {gh} k_E)$ (i.e. the time required to cross the inverse, effective wavenumber 
$1/k_E$, or ‘wave-crossing time’) did not qualitatively modify the results, so we do not utilize such a ramp-up here. The spectral solver results in high numerical accuracy, with the normalized root-mean-square difference between the unforced (
$P=0$) profile 
$\eta _0$ at 
$t'_1=10$ and the initial condition 
$\eta ^{(0)}_0$ is 
$2\times 10^{-13}$, and the normalized wave height change is 
$1-\left [\max (\eta _0) - \min (\eta _0)\right ]/ \left [\max (\eta ^{(0)}_0) - \min (\eta ^{(0)}_0)\right ] = -1\times 10^{-13}$.
We quantify the wave shape with the wave's energy 
$E$, skewness 
$\textrm {Sk}$ and asymmetry 
$\textrm {As}$,
\begin{align} E := \langle \eta_0^2 \rangle , \quad \textrm{Sk} := \frac{\langle \eta_0^3 \rangle}{\langle \eta_0^2 \rangle^{3/2}} \quad \text{and}\quad \textrm{As} := \frac{\langle \mathcal{H} \left\{\eta_0^3\right\} \rangle}{\langle \eta_0^2 \rangle^{3/2}} , \quad \text{with}\ \langle \,f \rangle := \frac{1}{L} \int_{{-}L/2}^{L/2} f \,\mathrm{d}x . \end{align}
Here, 
$\mathcal {H}(f)$ is the Hilbert transform of 
$f$, defined as the imaginary part of 
$\mathcal {F}^{-1} \left ( \mathcal {F} \left (f\right ) 2U\right )$ with 
$U$ the unit step function and 
$\mathcal {F}$ the Fourier transform. Since these definitions depend on the domain size 
$L$, we normalize the energy 
$E$ and skewness 
$\textrm {Sk}$ by their initial values.
3. Results
We study the pressure magnitude's effect on solitary-wave evolution and shape by varying the KdV–Burgers equation's (2.18) one free parameter, 
$P k_E/(\rho _w g \varepsilon )$, with emphasis on the contrast between onshore (
$P > 0$) and offshore wind (
$P < 0$). We revert to denoting non-dimensional variables with primes and dimensional ones without.
The wave profile 
$\eta /h$ snapshots in figure 1 qualitatively show how the wave shape evolves over non-dimensional slow time 
$t'_1 = t \varepsilon \sqrt {g h} k_E$ in the unforced solitary wave's frame. The onshore wind generates wave growth, apparent at the wave crest (figure 1a), whereas the offshore wind causes decay (figure 1b). The wind also changes the phase speed, with the wave's acceleration (deceleration) under an onshore (offshore) wind visible by the advancing (receding) of the crest. This is expected due to the (unforced) solitary wave's nonlinear phase speed (2.20) dependence on the wave height 
$H$.
Figure 1. Solitary-wave evolution under (a) onshore and (b) offshore wind-induced surface pressure in the frame of the unforced solitary wave. Non-dimensional wave height 
$\eta /h$ versus non-dimensional distance 
$x/h$ for 
$\varepsilon =0.1$, 
$\mu _E = 0.6$, 
$\left |P k_E/(\rho _w g \varepsilon )\right | = 0.25$ and non-dimensional slow times 
$t'_1 = t \varepsilon \sqrt {gh} k_E = 0$, 
$5$ and 
$10$, as indicated in the legend. Only a subset of the full spatial domain is shown. The arrows denote the wave propagation (phase speed) and wind direction.
In shallow water, wave growth/decay and phase speed changes are well-known wind effects (e.g. Miles Reference Miles1957; Cavaleri & Rizzoli Reference Cavaleri and Rizzoli1981), but wind-induced wave shape changes (Zdyrski & Feddersen Reference Zdyrski and Feddersen2020) have not been previously studied for shallow-water systems. Such changes are visible in figure 1 where, despite the wave starting from a symmetric, solitary-wave initial condition, the wind induces a horizontal asymmetry in the wave shape, particularly on the rear face (
$x<0$) of the wave. The offshore wind (figure 1b) raises the rear base of the wave (near 
$x/h = -5$) relative to its initial profile (purple line), but the onshore wind (figure 1a) depresses the rear face and forms a small depression below the still water level at 
$t\varepsilon \sqrt {gh} k_E=5$ (blue line) which widens and deepens at 
$t\varepsilon \sqrt {gh} k_E=10$ (green line). Finally, the onshore wind (figure 1a) increases the maximum wave-slope magnitude with time while the offshore wind (figure 1b) decreases it, although the windward side of the wave becomes steeper than the leeward side for both winds (up to 
$8\,\%$ steeper for the time period shown). Although the equation is ill posed in the sense of Hadamard, the smooth solutions show that our solution is acceptable up to the current time and thus we are justified in neglecting a regularization scheme.
To further examine the wind-induced wave asymmetry, we fit 
$\eta$ to a reference solitary-wave profile 
$\eta _{{ref}}$ (2.19) by minimizing the 
$L_1$ difference, yielding the reference height 
$H_{{ref}}(t_1)$ and peak location 
$x_{{ref}}(t_1)$. The profile change is defined as 
$\Delta \eta (x) := \eta - \eta _{{ref}}$ and is shown as a function of the reference wave-centred distance 
$\tilde {x} := x - x_{{ref}}$ in figure 2. Notice that the profile change begins near the front face of the wave and has extrema for negative 
$\tilde {x}'$ but with opposite signs for onshore and offshore winds. Additionally, the magnitude of the extrema decay with distance in the 
$-\tilde {x}$ direction. Finally, note that the onshore (offshore) wind generates a small peak (trough) at 
$\tilde {x} = 0$ and two small troughs (peaks) near 
$\tilde {x}/h = {\pm }3$, with the 
$\tilde {x}<0$ extrema larger than the 
$\tilde {x}>0$ one. This is analogous to a dispersive tail, well known in KdV-type systems (e.g. Hammack & Segur Reference Hammack and Segur1974), and its appearance here helps explain the pressure-induced shape change (cf. § 4.2).
Figure 2. The non-dimensional profile change 
$\Delta \eta /h$ between the surface profile and reference solitary wave (2.19) under (a) onshore and (b) offshore Jeffreys forcing versus non-dimensional reference wave-centred distance 
$\tilde {x}/h$. Results are shown for 
$\varepsilon =0.1$, 
$\mu _E = 0.6$, 
$\left |P k_E/(\rho _w g \varepsilon )\right | = 0.25$ and non-dimensional slow times 
$t'_1 = t \varepsilon \sqrt {gh} k_E = 0$, 
$5$ and 
$10$, as indicated in the legend. Only a subset of the full spatial domain is shown. The arrows denote the direction of wave propagation (phase speed) or wind direction.
The effect of wind on wave shape is quantified by the time evolution of wave shape statistics – energy, skewness and asymmetry – for onshore and offshore wind (figure 3). We plot all cases for initial steepness 
$\varepsilon = 0.1$ up to slow time 
$t \varepsilon \sqrt {g h} k_E = 10$, corresponding to 
$10/\varepsilon = 100$ wave-crossing times, 
$T_0 = 1/(\sqrt {gh} k_E)$. The unforced case (
$P=0$) displays constant shape statistics and zero asymmetry, as expected. The normalized energy 
$E/E_0$ shows different growth/decay rates: the onshore wind (
$P>0$) causes accelerating wave growth while the offshore wind (
$P<0$) causes slowing wave decay (figure 3a). The energy of the unforced wave is virtually unchanged, with a normalized energy change of 
$1-E/E_0 = -1\times 10^{-13}$ at 
$t'_1 = 10$. The onshore (offshore) wind causes the wave to become more (less) skewed over time, with the normalized skewness nearly symmetric about unity with respect to 
$\pm P$. Finally, the onshore wind causes a backwards tilt and negative asymmetry while the offshore wind increases the asymmetry and causes a forward tilt, which was also seen in figure 1. Notice that 
$\left |\textrm {As}\right |$ is larger for onshore winds than offshore winds. Since the definitions of the skewness and asymmetry are insensitive to waveform scaling 
$\eta \to \lambda \eta$, this effect is not simply caused by the wave's growth/decay. Instead, the onshore wind generates a larger dispersive tail (figure 2), which is the asymmetric wave component.
Figure 3. Solitary-wave shape statistics under onshore and offshore Jeffreys forcing versus non-dimensional slow time 
$t'_1 = t \varepsilon \sqrt {gh} k_E = 0\text {--}10$. The (a) energy (normalized by the initial energy), (b) skewness (normalized by the initial skewness) and (c) asymmetry are defined in (2.21a–c). Results are shown for 
$\varepsilon =0.1$, 
$\mu _E = 0.6$ and pressure magnitude 
$|P k_E/(\rho _w g \varepsilon )|$ up to 
$0.25$, as indicated in the legend. The solid black line is the unforced case, 
$P = 0$, and shows no growth or asymmetry and a constant skewness.
4. Discussion
4.1. Wind speed estimation
We now relate the non-dimensional pressure magnitude 
$P k_E/(\rho _w g) = {O}(\varepsilon )$ to the wind speed. First, we need a relationship between the surface pressure and wave energy 
$E$ (2.21a–c), which we can approximate using the standard procedure (e.g. Mei et al. Reference Mei, Stiassnie and Yue2005) of multiplying the (non-dimensional, denoted by primes) KdV–Burgers equation (2.18) by 
$\eta '_0$ and integrating from 
$x'=-\infty$ to 
$\infty$ to obtain
\begin{equation} \frac{\partial}{\partial t'_1} \int_{-\infty}^{\infty} {\eta'}_0^2 \,\mathrm{d}x' = \int_{-\infty}^{\infty} P' \left(\frac{\partial \eta'_0}{\partial x'}\right)^2 \,\mathrm{d}x' . \end{equation}
The left integral is the non-dimensional energy (2.21a–c), so re-dimensionalizing and converting back to the full time 
$t$ gives the energy growth rate 
$\gamma$,
\begin{equation} \frac{\gamma}{c_0 k_E} := \frac{1}{c_0 k_E E} \frac{\partial E}{\partial t} = \frac{P k_E}{\rho_w g} \frac{\langle (\partial_x \eta)^2 \rangle} {\langle (k_E \eta)^2 \rangle} = \frac{1}{5} \frac{P k_E}{\rho_w g} , \end{equation}
with 
$\langle (\partial _x \eta )^2 \rangle / \langle (k_E \eta )^2 \rangle = 1/5$ evaluated with the initial, solitary-wave profile (2.19) and the linear, shallow-water phase speed 
$c_0 = \sqrt {gh}$ coming from the re-dimensionalization of 
$t' = t c_0 k_E$ (2.7). Alternatively, a secondary multiple-scale approximation of the forward KdV–Burgers equation has been used previously to derive the energy growth rate for solitary waves as (Zdyrski & McGreevy Reference Zdyrski and McGreevy2019)
\begin{equation} E \propto \frac{1}{\left(1 - \gamma t\right)^2} \quad \text{with}\ \gamma := b \left[\frac{P k_E}{\rho_w g}\right] c_0 k_E , \end{equation}
with analytically derived 
$b = 2/15$. Numerically fitting (4.3) to our calculated energy instead yields 
$b = 0.10081 \pm 0.00003$, similar to the analytic approximation. Note that the exponential energy growth (4.2) correctly approximates (4.3) for small times 
$\gamma t \ll 1$, and both expressions are consistent with the observed accelerating (decelerating) energy change for 
$P>0$ (
$P<0$) in figure 3.
Next, Jeffreys's (Reference Jeffreys1925) theory relates the growth rate of periodic waves to the wind speed 
$U_{\lambda /2}$, measured at a height equal to half the wavelength 
$z=\lambda /2$, as
\begin{equation} \frac{\gamma}{c k} = S_{\lambda/2} \frac{\rho_a}{\rho_w} \left(\frac{U_{\lambda/2}}{c}-1\right) \left|\frac{U_{\lambda/2}}{c}-1\right| , \end{equation}
with 
$S_{\lambda /2}$ a small, non-dimensional sheltering parameter potentially dependent on 
$\varepsilon$, 
$\mu _E$ and 
$U_{\lambda /2}/c$. For simplicity, we approximate the nonlinear phase speed 
$c$ (given non-dimensionally in (2.20)) by its leading-order term 
$c_0 = \sqrt {gh}$, yielding an error of only 
$10\,\%$ in the subsequent calculations. Combining this approximation of (4.4) with (4.2) gives
\begin{equation} U_{\lambda/2} = c_0 \left(1 \pm \sqrt{\frac{1}{5} \left|\frac{P k_E}{\rho_w g}\right| \frac{\rho_w}{\rho_a} \frac{1}{S_{\lambda/2}}}\right) . \end{equation}
Here, the 
$\pm$ corresponds to onshore (
$+$) or offshore (
$-$) winds. Note that changing the wind direction (i.e. 
$\pm$ sign) while holding the surface pressure magnitude 
$\left |Pk_E/(\rho _w g)\right |$ constant means onshore wind speeds 
$\left |U_{\lambda /2}\right |$ will be larger than offshore wind speeds.
We can evaluate (4.5) for the parameters of § 3: 
$\varepsilon =0.1$, 
$\mu _E = 0.6$ and 
$Pk/(\rho _w g \varepsilon ) = 0.25$. Donelan et al. (Reference Donelan, Babanin, Young and Banner2006) parameterized 
$S_{\lambda /2}$ for periodic shallow-water waves with a dependence on airflow separation: 
$S_{\lambda /2} = 4.91 \varepsilon \sqrt {\mu }$ for our non-separated flow (according to their criterion), with 
$\mu := (kh)^2$. Assuming this holds approximately for solitary waves, we choose 
$\lambda = 2 {\rm \pi}/k_E = 20\ \textrm {m}$ to calculate the wind speed at 
$z = \lambda /2 = 10\ \textrm {m}$. This choice corresponds to a depth of 
$h = 2.5\ \textrm {m}$ and initial wave height 
$H_0 = 0.5\ \textrm {m}$ and yields a wind speed of 
$U_{10} = 22\ \textrm {m}\,\textrm {s}^{-1}$, a physically realistic wind speed for strongly forced shallow-water waves. Weaker wind speeds will induce smaller surface pressures and thus take longer to change the wave shape.
4.2. Physical mechanism of asymmetry generation
Our initial, symmetric solitary waves (2.19) are permanent-form solutions of the unforced KdV equation. More generally, any initial solitary wave which does not exactly solve the KdV equation will evolve into a solitary wave and a trailing, dispersive tail according to the inverse scattering transform (e.g. Mei et al. Reference Mei, Stiassnie and Yue2005). In our system, the pressure continually perturbs the system away from the unforced KdV soliton solution resulting in a trailing, bound, dispersive tail (figure 2), which is responsible for the wave asymmetry. To see this, consider an initial, symmetric profile 
$\eta$. The pressure forcing term 
$P \partial _x^2 \eta$ preserves the initial symmetry and induces a symmetric bound wave after a short time 
$\Delta t'_1 \ll 1$. This is apparent when considering the non-dimensional KdV–Burgers equation (2.18) in the unforced solitary wave's frame (figure 1) at the initial time,
\begin{gather} \left.\frac{\partial \eta'_0}{\partial t_1}\right|_{t'_1=0} ={-}P' \frac{\partial^2}{\partial x^2} \left[\mathrm{sech}^2\left(\frac{x'}{2}\right)\right] \end{gather}
\begin{gather}\implies \eta'_0(x', \Delta t'_1) = (2-P'\Delta t'_1) \,\mathrm{sech}^2\left(\frac{x'}{2}\right) + P' \Delta t'_1 \frac{3}{2}\,\mathrm{sech}^4\left(\frac{x'}{2}\right) . \end{gather}
The 
$P'\Delta t'_1$ terms generate a small bound wave with a peak (trough) at 
$x'=0$ and troughs (peaks) symmetrically in front and behind the wave peak for onshore (offshore) wind. As time increases, the continual pressure forcing causes the bound wave to grow and lengthen behind the wave, as is apparent in figure 2 (e.g. 
$\tilde {x}/h = -20\text { to }3$ for 
$P'=0.25$ and 
$t'_1=10$).
The small numerical value 
$\left |P'\right | = 0.25 \ll 1$ used in § 3 allows us to consider the wave's evolution as two steps with time scale separation. First, the pressure generates a bound wave (4.7) on the slow time scale, and then the wave evolves a dispersive tail on the fast time scale according to the inverse scattering transform of the unforced KdV equation. The dispersive tail in figure 2 (e.g. located left of 
$\tilde {x}/h=-20$ for 
$P'=0.25$ and 
$t'_1=10$) is analogous to the ubiquitous dispersive tails in prior studies on shallow-water solitary waves, such as figures 8(b) and 8(c) of Hammack & Segur (Reference Hammack and Segur1974). However, unlike dispersive tails generated from initial conditions which fail to satisfy the KdV equation, our tail is continually forced and lengthened by the wind forcing. Finally, interactions with the trailing, dispersive tail are responsible for lengthening the bound wave (4.7) behind, rather than ahead, of the solitary wave. Hence, the disturbance induced by the pressure forcing (4.7) has two effects on the wave. First, the wind slowly generates a bound wave which changes the height and width of the initial solitary wave, which is reflected in the growth (decay) and narrowing (widening) under onshore (offshore) winds in figure 1. Second, it quickly generates an asymmetric, dispersive tail behind the wave (figure 2), producing a greater shape change on the wave's rear face (figure 1). Finally, the different wind directions (i.e. pressure forcing signs) change the sign of the bound wave and dispersive tail and, hence, the sign of the asymmetry in figure 3.
4.3. Comparison to intermediate and deep water
Zdyrski & Feddersen (Reference Zdyrski and Feddersen2020) investigated the effect of wind on Stokes-like waves in intermediate to deep water. This study, with wind coupled to waves in shallow water, finds qualitative agreement with those intermediate- and deep-water results. The shallow-water asymmetry magnitude increases as the pressure magnitude 
$P$ increases (figure 3), and figure 4(a) of Zdyrski & Feddersen (Reference Zdyrski and Feddersen2020) displayed a similar trend for the corresponding Jeffreys pressure profile, with positive (negative) pressure increasing (decreasing) the asymmetry. Although Zdyrski & Feddersen (Reference Zdyrski and Feddersen2020) compared their theoretical predictions to limited experimental results with 
$kh > 1$, there are no appropriate experiments on wind-induced changes to wave shape in shallow water for comparison with our results. In addition to the Jeffreys pressure profile employed here, Zdyrski & Feddersen (Reference Zdyrski and Feddersen2020) also utilized a generalized Miles profile, only applicable to periodic waves, wherein the pressure was proportional to 
$\eta$ shifted by a distance parameter 
$\psi _P/k$. Future investigations could couple a higher-order Zakharov equation (e.g. Dommermuth & Yue Reference Dommermuth and Yue1987) to a Jeffreys-type pressure forcing or to an atmospheric large eddy simulation, as was done for deep water by Hao & Shen (Reference Hao and Shen2019). Although this analysis focuses on solitary waves, we also investigated the effect of wind on periodic waves using the cnoidal-wave KdV solutions as initial conditions. Wind-forced cnoidal waves displayed qualitatively similar shape changes with stronger onshore (offshore) wind causing the energy and skewness to increase (decrease) while the asymmetry decreased (increased) with time. Furthermore, results were qualitatively similar across multiple classes of cnoidal waves with different values of 
$\mu _E$, implying that these results apply rather generally.
5. Conclusion
Prior results (Zdyrski & Feddersen Reference Zdyrski and Feddersen2020) in intermediate and deep water demonstrated that wind, acting through a wave-dependent surface pressure, can generate shape changes that become more pronounced in shallower water. Here, we produced a novel analysis of wind-induced wave shape changes in shallow water using a multiple-scale analysis to couple weak wind with small, shallow-water waves, i.e. 
$a_0/h \sim (k_E h)^2 \sim P k/(\rho _w g) \ll 1$. This analysis produced a KdV–Burgers equation governing the wave profile 
$\eta$, which we then solved numerically with a symmetric, solitary-wave initial condition. The deviations between the numerical results and a reference solitary wave had the form of a bound, dispersive tail, with differing signs for onshore and offshore wind. The tail's presence and shape are the result of a symmetric, pressure-induced shape change evolving under the inverse scattering transform. We also estimated the energy, skewness and asymmetry as functions of time and pressure magnitude. For onshore wind (positive 
$P$), the wave's energy and skewness increased with time while asymmetry decreased, while offshore wind produced the opposite effects. Furthermore, these effects were enhanced for strong pressures, and they reduced to the unforced case for 
$P=0$. The shape statistics found here show qualitative agreement with the results in intermediate and deep water. Finally, the wind speeds corresponding to these pressure differences were calculated and found to be physically realistic.
Acknowledgements
We are grateful to D.G. Grimes and M.S. Spydell for discussions on this work. Additionally, we thank the anonymous reviewers for their suggestions and comments.
Funding
We thank the National Science Foundation (OCE-1558695) and the Mark Walk Wolfinger Surfzone Processes Research Fund for their support of this work.
Declaration of interests
The authors report no conflict of interest.
Open access
