Articles | Volume 14, issue 6
Research article
05 Dec 2023
Research article |  | 05 Dec 2023

The effect of strong nonlinearity on wave-induced vertical mixing

Maciej Paprota and Wojciech Sulisz

A semi-analytical solution to an advection–diffusion equation is coupled with a nonlinear wavemaker model to investigate the effect of strong nonlinearity on wave-induced mixing. The comparisons with weakly nonlinear model predictions reveal that in the case of waves of higher steepness, enhanced mixing affects the subsurface layer of the water column. A fully nonlinear model captures the neglected higher-order terms from a weakly nonlinear solution and provides a reliable estimation of the time-mean velocity field. The corrected wave-induced mass-transport velocity leads to improved estimates of subsurface mixing intensity and ocean surface temperature.

1 Introduction

Mass-transport processes associated with the propagation of non-breaking ocean surface waves strongly affect the mixing of oceanic waters and the global exchange of heat at the air–water interface. Surface waves transfer energy into turbulence, modifying the mixing intensity of the upper ocean (Qiao et al.2004, 2010, 2013). Therefore, the correct identification and quantification of mass-transport processes associated with water waves leading to mixing of subsurface ocean waters is of practical importance for short-term and long-term weather forecasts.

Small-scale and large-scale climate modelling equally benefits from including the wave-induced mixing predictions into the general ocean circulation simulations (e.g. Qiao et al.2004; Xia et al.2006). Including the parameterized wave-induced mixing in ocean circulation models confirms that its contribution qualitatively improves the reliability of numerical results (e.g. Song et al.2007; Shu et al.2011). Despite the fact that some field measurements lead to the conclusion that surface waves contribute to heat exchange of the upper ocean (Matsuno et al.2006), it is difficult to collect reliable in situ data confirming the role of non-breaking surface waves in vertical mixing processes. An alternative approach to the problem is to investigate the wave-induced mixing based on the physical tests in a wave flume or a wave basin.

Laboratory experiments are the source of valuable information on wave-induced mixing processes as they provide data enabling the evaluation of wave-induced mixing coefficient (e.g. Babanin and Haus2009; Dai et al.2010; Sulisz and Paprota2015). In a wave flume, the mixing originated from breaking and non-breaking waves may be investigated to evaluate the real contribution from the non-breaking surface waves to general ocean circulation, which is of fundamental importance for calibration of wave-induced mixing models (Sulisz and Paprota2019). It should be noted that the experimental data are affected by undesirable spurious laboratory effects requiring special attention when analysing the results of mass-transport processes driven by waves (e.g. Paprota and Sulisz2018; Paprota2020).

The study extends the analysis of wave-induced vertical mixing performed by Sulisz and Paprota (2019), which used a weakly nonlinear theory applied to mechanically generated water waves. The improvements cover the more exact calculations of the wave velocity field using a fully nonlinear wavemaker model, which was successfully verified with respect to kinematics of regular waves against laboratory measurements collected in the flume (Paprota and Sulisz2018; Paprota2020). This method allows modelling of non-breaking waves with strong nonlinearities and admits amplitude dispersion, nonlinear wave–wave interactions in deep and intermediate waters, and solitary wave propagation, which go beyond the applicability of weakly nonlinear approaches. The model belongs to a family of wave solutions based on a pseudo-spectral approach, which consider highly nonlinear non-breaking waves and originate from methods derived by Dommermuth and Yue (1987) and West et al. (1987) (see also Paprota and Sulisz2019, for a review).

The application of the fully nonlinear model should lead to the more accurate estimation of the phase-averaged wave velocity field and, hopefully, the more reliable evaluation of the evolving water temperature field under regular waves. Since the numerical model provides a full description of the evolution of the velocity field, it is possible to separate the Stokes drift from Lagrangian and Eulerian mean velocities (see Paprota et al.2016; Paprota and Sulisz2018; Paprota2020). Hence, the results presented in this study may improve simple models based on Stokes drift applicable to ocean waves of random (spectral) nature (Myrhaug et al.2018), while the derived modelling framework may be modified to cover open-ocean hydrodynamics with other forms of introducing waves (Paprota2019).

New contributions are also presented with regard to methods of calculations of time-mean flows, which are essential for the present analyses. Namely, a new Eulerian procedure is reported, which is different to the Lagrangian-based method applied in earlier studies (Paprota et al.2016; Paprota and Sulisz2018; Paprota2020). The study also highlights the differences between two methods of calculation of velocity distribution of time-mean wave-induced flows based on either Lagrangian particle tracking or approximated Eulerian averaging. In the case of the former, an improved and more accurate procedure of estimating the phase-averaged velocity is developed.

This study presents the methods of an analysis of wave-induced vertical mixing from the wavemaker perspective. The derived mathematical modelling approach may be directly applied to experimental activities aiming at recognizing, quantifying and parameterizing wave-induced vertical mixing effects (see e.g. Dai et al.2010; Sulisz and Paprota2015) for a better design of physical model set-up and a wave parameter selection in the laboratory. The awareness of the influence of laboratory side effects on the experimental outcome is essential for the correct interpretation of the results. This issue was raised in earlier studies (Sulisz and Paprota2015, 2019).

The paper is composed as follows. First, an outline of the coupled theoretical model describing the mixing processes under mechanically induced water waves is presented. Then, a comprehensive comparison between weakly and fully nonlinear approaches is given using numerical results to evaluate the effects of strong nonlinearity on vertical mixing processes. Finally, a discussion on the major results is provided together with remarks on accuracy and reliability of mathematical and numerical methods together with further discussions on putting the results of the study in a broader context of earth system modelling with respect to general ocean circulation. The paper is then completed by a summary and conclusions.

2 Materials and methods

The considered numerical approach to the modelling of vertical mixing induced by mechanically generated waves is realized through a procedure involving two fundamental steps referring to the solution of the wavemaker problem and advection–diffusion balance, respectively.

2.1 Particle kinematics of mechanically generated waves

First, the problem of the generation of waves in a numerical wave flume is formulated and solved. In the present study, a potential flow wave theory is used to obtain the solution within the Eulerian frame of reference. Weakly nonlinear analytical and higher-order numerical methods are employed to determine the wave fields in the rectangular domain, in which the water elements are defined by the horizontal x and vertical z coordinates of the Cartesian system. The origin of the system is located at the intersection of the wavemaker zero position and free-surface level corresponding to the hydrostatic conditions. The mechanically driven oscillation of the free surface represented by elevation function η(x,t) is induced by the piston-like motion of the wavemaker paddle according to the displacement function χ(t). The flume bottom is assumed horizontal, and the water depth is h=const. The general presentation of the computational domain and the location of the coordinate system are depicted in Fig. 1.

Figure 1Schematic view of mechanically generated waves and the coordinate system.


According to the potential flow assumptions, the motion of an inviscid and incompressible fluid is irrotational. Moreover, the solid boundaries are impervious. The scalar velocity potential function ϕ(x,z,t) may be introduced to determine the velocity vector field v(x,z,t) such that v=∇ϕ. The wavemaker boundary-value problem is then formulated as


where g is the acceleration due to gravity.

The first part of a pseudo-spectral solution involves expanding the kinematic free-surface boundary condition (2), the dynamic free-surface boundary condition (3) and the kinematic wavemaker boundary condition (5) in a Taylor series about a mean position corresponding to the still water level (z=0) for Eqs. (2)–(3) and wavemaker paddle zero position (x=0) for Eq. (5):


In this way, a simple rectangular form of the computational domain is preserved, and the solution procedure is advanced further using either the perturbation or spectral approach. In the present study both weakly nonlinear (perturbation expansions) and fully nonlinear (spectral expansions) solutions are briefly presented and applied to calculate the velocity field for wave-induced mixing calculations. The order of nonlinearity depends on the order of nonlinear terms of Taylor series expansions of boundary conditions. Second-order terms are retained for weakly nonlinear solutions, while the fully nonlinear solution is provided with no upper limit (arbitrary order).

2.1.1 Weakly nonlinear solution

The perturbation expansion is first applied to solve the wavemaker problem defined by Eqs. (1)–(5) according to the solution derived by Hudspeth and Sulisz (1991) and Sulisz and Hudspeth (1993). The monochromatic wavemaker paddle displacement of amplitude s is assumed to be

(9) χ ( t ) = s sin ( σ t + φ ) ,

which generates periodic waves of the first-order amplitude of a, the angular frequency of σ and the phase φ in the semi-infinite flume domain. Due to the fact that x goes to infinity, the radiation condition (outward propagating waves) is imposed at the far-end lateral boundary of the domain (Hudspeth and Sulisz1991).

Using the expansions of the boundary conditions Eqs. (6)–(8) and retaining the terms up to the second order, the weakly nonlinear boundary conditions become


Additionally, the following small steepness parameter ϵ=ak (where k is the first-harmonic wave number) perturbation expansions of the angular frequency ω (σ=ω0), the free-surface elevation and the velocity potential functions are used (Hudspeth and Sulisz1991):


Substituting the perturbation forms Eqs. (13)–(15) into the boundary conditions correct up to the second order, Eqs. (10)–(12), leads to the final weakly nonlinear solution to the velocity potential and free-surface elevation functions (Hudspeth and Sulisz1991). The formulas for the horizontal U(x,y) and vertical W(x,y) components of the time-independent mass-transport velocity V(x,y) for the case of the piston-type wavemaker of full-depth draught are calculated from nondimensional forms reported as Eqs. (50a) and (50b) by Hudspeth and Sulisz (1991). In order to get dimensional values, the results calculated by Eqs. (50a) and (50b) in the work of Hudspeth and Sulisz (1991) are multiplied by agk.

Far away from the wavemaker paddle (x>3h), the vertical component of the time-independent velocity vanishes, and the time-independent horizontal velocity profile along the water depth UL converges to the sum of the Stokes drift US and return current UE velocities (Longuet-Higgins1953; Dean and Dalrymple1984), i.e.

(16) U L ( z ) = U S ( z ) + U E ,

where the Stokes drift profile is calculated as

(17) U S ( z ) = k σ a 2 cosh ( 2 k ( z + h ) ) 2 sinh 2 ( k h ) ,

and the return current value takes the form

(18) U E = - k σ a 2 sinh ( 2 k h ) 4 k h sinh 2 ( k h ) = - a 2 k g 2 σ h .

2.1.2 Fully nonlinear solution

An alternative approach, which admits higher-order wave components, is based on the spectral method applied to the wavemaker problem (Paprota and Sulisz2018). Contrary to the presented perturbation approach, the waves are generated by an arbitrary function χ(t) in a finite domain of length b, which is sufficiently large to exclude the effect of wave reflection from the far-end wall of the flume. In this regard, only progressive waves are considered to facilitate the comparisons with the weakly nonlinear solution.

The free-surface elevation function is expanded in a Fourier cosine series as

(19) η ( x , t ) = i = 0 a i cos ( λ i x ) ,

while the corresponding expansion of the velocity potential function is coupled with an additional term satisfying the wavemaker boundary condition (5) to give

(20) ϕ ( x , z , t ) = i = 0 A i cosh ( λ i ( z + h ) ) cosh ( λ i h ) cos ( λ i x ) + B 0 ( ( x - b ) 2 - ( z + h ) 2 ) + j = 1 B j cosh ( μ j ( x - b ) ) cosh ( μ j b ) cos ( μ j ( z + h ) ) ,

where λi=iπ/b and μj=jπ/h are the eigenvalues of the expansions, and the solution coefficients are the functions of time (ai(t), Ai(t) and Bj(t)). The Bj coefficients are (Paprota and Sulisz2018)

(21) B 0 = χ t h + m = 1 M χ m m ! sin m π 2 i = 0 λ i m A i tanh ( λ i h ) 2 h ( χ - b ) ,

(22) B j = - 2 cos ( μ j h ) m = 1 M χ m m ! sin ( m π 2 ) i = 0 λ i m + 2 λ i 2 + μ j 2 A i tanh ( λ i h ) h m = 0 M μ j m + 1 χ m m ! cos ( m π ) tanh ( μ j b ) | cos m π 2 | , j > 0 .

The coefficients Ai, Bj and ai are determined in an iterative solution procedure from the kinematic free-surface boundary condition (6), the dynamic free-surface boundary condition (7) and the kinematic wavemaker boundary condition (8). For a given time t, the unknown coefficients ai and Ai are calculated using a Fourier transform of η and ϕ as

(23) a i = 2 b 0 b η ( x , t ) cos λ i x d x , A i = 2 b 0 b ϕ ( x , 0 , t ) cos λ i x d x ,

while Bj values are determined using Eqs. (21) and (22). The coefficients are then used to calculate free-surface and velocity potential values using inverse Fourier transforms to advance the solution in time. Hence, time stepping is applied in a physical domain to obtain values of ϕ and η at a new time level. A fourth-order Adams–Bashforth–Moulton predictor–corrector approach is preferred as a time-marching scheme (see e.g. Press et al.1988), with initial values of ϕ(x,0,0)=0 and η(x,0)=0. The wavemaker model used in the present study is reported in the work of Paprota and Sulisz (2018).

The instantaneous wave velocity field is derived by expressing the horizontal u and vertical w velocity components in terms of the spatial derivatives of the velocity potential ϕx and ϕz, respectively, as


After a fully developed wave motion is achieved, the time-independent wave velocity field may be approximated by time averaging of the instantaneous wave velocity (see Hudspeth and Sulisz1991) within the range limited by two in-phase states of regular wave motion (over one wave period, which is either given or calculated from the dispersion relation) as


and is referred to here and after as the Eulerian-mean transport velocity (EMTV). The pair of triangle brackets 〈 〉 denotes the operator of time averaging over one wave period. The velocities u and w are calculated using Eqs. (24) and (25), while their derivatives are evaluated analytically. The integrals of u and w are determined directly from Eqs. (24) and (25) upon replacing the time-dependent coefficients Ai and Bj with their integrals Aidt and Bjdt. The integration is achieved by expanding Ai and Bj into a Fourier series with respect to time and integrating the resulting Fourier expansions analytically.

An alternative approach, which leads to the time-independent wave velocity field, is the procedure involving the Lagrangian particle tracking: Lagrangian-mean transport velocity (LMTV). The time-averaged velocity is calculated based on the displacement of a water particle moving between its two successive in-phase positions along the particle trajectory (Paprota et al.2016). The trajectory of a water particle is determined by numerical integration of the system of differential equations,


for the set of initial particle locations.

In the present study, the improvements to the method of evaluation of mass-transport velocity based on the Lagrangian particle tracking (Paprota and Sulisz2018) are introduced. In the previous works (Paprota et al.2016; Paprota and Sulisz2018; Paprota2020), the procedure relied on distributing artificial tracers in a water column under fully developed regular waves. Initially, the tracers are uniformly distributed along the water depth, and they coincide with the zero down-crossing phase of the wave for a given distance from the wavemaker. After one Lagrangian wave period (see e.g. Longuet-Higgins1986; Chen et al.2009), the tracers move from their original position due to mass transport. The particle tracking procedure is employed to determine the mass-transport velocity profile. The procedure is repeated for subsequent longitudinal positions to cover the accepted region of interest.

Here, in order to get the better estimation of the time-independent velocity field, the two hydrodynamic states corresponding to both zero up- and down-crossings of the regular wave are used to start up the tracking procedure – contrary to the previous method based only on the zero down-crossings (e.g. Paprota et al.2016; Paprota and Sulisz2018). The procedure is analogous but involves more tracers and covers two phase positions. The resulting LMTV is calculated as the mean of the time-independent velocity fields corresponding to both zero-crossing initial states.

2.2 Wave-induced vertical mixing intensity

The numerical solution to a two-dimensional advection–diffusion equation with relevant boundary conditions is used to predict the evolution of the temperature field under mechanically generated regular waves in the flume. Assuming that the changes in temperature are solely due to wave-induced and diffusion processes, the following boundary-value problem is formulated as in the work of Sulisz and Paprota (2019):


where T is the water temperature, and κ is the diffusion coefficient represented as a sum of the molecular (κm) and wave-induced (κv) diffusivities.

(35) κ = κ m + κ v

In the advection–diffusion equation (30) the velocities U an W are, in fact, Lagrangian velocities which describe the mean transport velocity field (the Stokes drift and a return current in this case). The inclusion of the Stokes drift (without return current) to the advection equation was discussed explicitly for the case of ocean waves by McWilliams and Sullivan (2000).

The molecular diffusivity κm=1.4×10-7 m s−2, while the wave-induced diffusivity is calculated by the formula reported by Sulisz and Paprota (2015) as

(36) κ v ( z ) = α a 3 k σ sinh 2 ( k ( z + h ) ) cosh ( k ( z + h ) ) sinh 3 ( k h ) ,

where α=0.002 is the dimensionless coefficient, which is evaluated based on measurements (Sulisz and Paprota2015). The derivation of κv is based on the weakly nonlinear theory. More information is provided in the work of Sulisz and Paprota (2015), where the comparison between the weakly nonlinear and more general forms is provided. The parameter α was estimated based on the experiments only for the presented form of κv. In order to use a more general formula, new values of α must be determined using experimental data from a wider range of wave conditions.

The advection–diffusion equation (Eq. 30) holds in the entire fluid domain, while the heat radiation is assumed zero at the water surface (Eq. 31), the bottom (Eq. 33) and the lateral boundaries (Eqs. 32 and 34). The length b is sufficiently long in order to reduce the effect of the finite domain on the results in the area of interest, which is limited to the region of several water depths from the wavemaker paddle. The omitted procedure of incorporating the non-zero heat radiation at the water surface is discussed in the paper by Sulisz and Paprota (2019).

The solution of the advection–diffusion equation (Eq. 30) is achieved by employing a similar methodology to that used to solve the wavemaker problem admitting higher-order nonlinearities presented in the previous section. Accordingly, the scalar temperature field function is expanded into a double Fourier series of the form (Sulisz and Paprota2019)

(37) T ( x , z , t ) = i = 1 j = 1 d i j cos ( λ i ( x - b ) ) cos ( μ j ( z + h ) ) .

Again, the solution of Eq. (30) is achieved by a time-stepping procedure, which, in this case, consists of an application of the spectral expansion method to describe T and the Runge–Kutta formulas to proceed in time (see e.g. Press et al.1988). Accordingly, the wave-induced diffusivity is determined using Eq. (36), while the velocity field is calculated by means of either the weakly nonlinear solution (Hudspeth and Sulisz1991) or a higher-order approach (Paprota and Sulisz2018) for selected wave parameters. The initial condition of given vertical temperature distribution T0(z) is used to start up the time-stepping procedure. The time derivative of temperature T is calculated from Eq. (30) with the aid of the expansion Eq. (37) after the coefficients dij are determined by applying a two-dimensional cosine fast Fourier transform. The implementation of this procedure is freely available in Paprota (2023).

3 Results and discussion

The evaluation of temperature evolution of an oscillating water body is analysed in a numerical wave flume environment. This approach provides a basis for the straightforward verification of the major outcome of the modelling procedures, which are presented in the previous section, against measurements in the real laboratory. On the other hand, the numerical model may be modified to cover pseudo-random ocean waves in a large periodic domain (see e.g. Paprota2019) to analyse the wave-induced vertical mixing processes in offshore conditions.

3.1 Numerical test cases

The waves are generated by a monochromatic wavemaker motion in a numerical flume. The transitional and shallow water wave cases are considered for three values of the depth-relative dimensionless parameter kh of 0.5, 1 and 2. The effect of wave height H is also studied and corresponds to the dimensionless steepness parameter H/h of 0.05, 0.1 and 0.2. In the case of the higher-order solution, the waves are generated starting from rest with the ramp function applied to the stroke of the wavemaker motion for the first five wave periods. As previously stated, the longitudinal size of the flume b is sufficiently large to exclude wave reflections from the analysis, and it corresponds to 10 wavelengths (L) for kh=2 and 20L in the remaining longer wave cases. The selected progressive wave parameters are the same as in previous studies on the modelling of wave-induced mixing in wave flumes (Sulisz and Paprota2019). After the fully developed state is achieved, the time-independent wave velocity field is determined. The resulting horizontal and vertical velocity components U and W are then substituted in Eq. (30), and the temperature field evolution is predicted. A summary of wave parameters is provided in Table 1 together with wave-induced diffusivity values at the surface in relation to its molecular counterpart.

Table 1Basic dimensionless parameters of mechanically generated regular waves.

Download Print Version | Download XLSX

It should be noted that the model may be applied to deep-water as well as shallow-water conditions, which are presented in the work of Sulisz and Paprota (2015). Here, the cases correspond to a limited range of kh from 0.5 to 2.0. In this way, the results from the previous work (Sulisz and Paprota2019) may be directly compared. Moreover, for deep-water conditions corresponding to kh>π, the mixing would completely sweep away the warmer water from the analysed part of the fluid domain for the selected time frame, and this case was omitted in the present study.

3.2 Phase-averaged velocity distribution

The accurate assessment of the time-independent velocity field is of significant importance in the modelling of temperature changes under an undulating water surface. The increasing steepness of generated waves intensifies wave-induced mixing processes but also changes the structure of heat flux distribution in the region occupied by the fluid. The analysis of phase-averaged velocity distribution in the direct vicinity of the wavemaker paddle helps to identify the effect of advection terms of Eq. (30) on the temperature evolution driven by waves.

Figure 2Time-independent velocity field under mechanically generated regular waves characterized by kh=0.5, (a) H/h=0.05, (b) H/h=0.1 and (c) H/h=0.2. Weakly nonlinear theory – black arrows and line, higher-order-theory EMTV – green arrows, LMTV – blue arrows. Higher-order vectors are shifted upwards for convenience of comparison.


Figure 3Time-independent velocity field under mechanically generated regular waves characterized by kh=1, (a) H/h=0.05, (b) H/h=0.1 and (c) H/h=0.2. Weakly nonlinear theory – black arrows and line, higher-order-theory EMTV – green arrows, LMTV – blue arrows. Higher-order vectors are shifted upwards for convenience of comparison.


Figure 4Time-independent velocity field under mechanically generated regular waves characterized by kh=2, (a) H/h=0.05, (b) H/h=0.1 and (c) H/h=0.2. Weakly nonlinear theory – black arrows and line, higher-order-theory EMTV – green arrows, LMTV – blue arrows. Higher-order vectors are shifted upwards for convenience of comparison.


In Figs. 24, time-independent velocity fields calculated by weakly nonlinear and higher-order methods are presented. The arrows represent the vectors of a phase-averaged velocity corresponding to mass transport induced by waves. Black arrows correspond to a weakly nonlinear solution, while green and blue arrows correspond to a higher-order solution and two methods of averaging – EMTV and LMTV, respectively (Figs. 24). The way the results are presented is in line with the expected velocity vector field pattern in the part of the fluid domain adjacent to the wavemaker (see Paprota and Sulisz2018, for a more detailed description). The fluid flow forms half of the circulation cell limited by the air–water interface (z/h=0), the wavemaker paddle (x/h=0) and the bottom (z/h=-1). The water mass flows with the direction of wave propagation near the surface, while the adverse flow pushes the water to the wavemaker paddle in the lower part of the water column. At the wavemaker paddle, the water flows vertically, forming the lateral boundary of the general circulation in the flume. Far away from the paddle, the vertical velocity components vanish, and the typical mass-transport velocity profile over depth emerges as a sum of the Stokes drift Eq. (17) and the return current velocity Eq. (18). Thus, the figures are complemented by the weakly nonlinear horizontal mass-transport velocity profile UL(z) calculated outside the direct vicinity of the wavemaker paddle Eq. (16), where the effect of evanescent standing waves generated by the wavemaker (see e.g. Dean and Dalrymple1984) may be neglected (x>3h). UL(z) is plotted at the right outer edge of the graph at x=3h for convenience of comparison. Vector plots are complemented by bubble charts of relative differences between corresponding vector magnitudes as presented in Figs. 57, where black circles denote the differences between LMTV and weakly nonlinear results, while filled green circles refer to the differences between LMTV and EMTV.

Figure 5Relative differences for time-independent velocity field under mechanically generated regular waves characterized by kh=0.5, (a) H/h=0.05, (b) H/h=0.1 and (c) H/h=0.2. LMTV vs. weakly nonlinear theory – black circles, LMTV vs. EMTV – filled green circles.


Figure 6Relative differences for time-independent velocity field under mechanically generated regular waves characterized by kh=1, (a) H/h=0.05, (b) H/h=0.1 and (c) H/h=0.2. LMTV vs. weakly nonlinear theory – black circles, LMTV vs. EMTV – filled green circles.


Figure 7Relative differences for time-independent velocity field under mechanically generated regular waves characterized by kh=2, (a) H/h=0.05, (b) H/h=0.1 and (c) H/h=0.2. LMTV vs. weakly nonlinear theory – black circles, LMTV vs. EMTV – filled green circles.


The results presented in Fig. 2 provide information on the time-independent velocity field predicted in the direct vicinity of the wavemaker paddle and the differences between the weakly nonlinear approach and the solution admitting higher-order terms in the most extreme nonlinear regime of shallow water (kh=0.5). Due to the fact that the considered wave cases are characterized by the highest values of the Ursell number, the differences between the two wavemaker models immediately appear even for the lowest amplitude of free-surface oscillations of H/h=0.05 (Fig. 2a). The highest differences in the range between 10 % and 20 % correspond to the velocities near the wavemaker paddle and the surface, respectively (Fig. 5a). By increasing the magnitude of the surface oscillations to H/h=0.1 and H/h=0.2 (Fig. 2b and c), the differences become higher in the larger area of the wavemaker paddle vicinity. The higher-order model predicts more intensive mass transport near the bottom and the surface. The increase in the subsurface velocity and the magnitude of the current near the bottom may even reach from 15 % to 30 % for H/h=0.1 and from 20 % to 40 % for H/h=0.2 (Fig. 5b and c, respectively).

Some important information may also be acquired from the comparison between two methods of averaging corresponding to either Lagrangian particle tracking (LMTV) or Eulerian averaging (EMTV) based on Eqs. (26) and (27), respectively. Although both methods of averaging provide consistent results for lower wave heights of the generated waves, EMTV results are less reliable in the case of the highest waves H/h=0.2 (Fig. 2c). It can be seen that the velocity vectors near the corner point determined by the intersection of the wavemaker paddle mean position and the surface (x/h=0, z/h=0) are unnaturally large as the velocity in this area (green arrows around the point (0,0) in Figs. 2c and 3c) is expected to vanish (see Grilli and Svendsen1990, for a comprehensive discussion on numerical methods' inaccuracy at corner points of a fluid domain with moving boundaries). The problem still persists when considering the subsurface velocity at some distance away from the wavemaker. The differences between the LMTV and EMTV may even reach 20 % in the velocity magnitude near the intersection between the wavemaker and the surface, while an even greater discrepancy is visible at corner points in most extreme longer-wave cases (Fig. 5c).

With increasing relative water depth (Figs. 3 and 4), the differences between weakly nonlinear and higher-order time-independent velocity are less pronounced but still significant for steeper waves. In the case characterized by kh=1, the discrepancies are generally less then 10 % for H/h=0.05 and 0.1, ranging from 10 % to 20 % for the highest waves (Fig. 6). Again the more intensive mass circulation is apparent when higher-order terms are taken into account for the case of kh=1 and H/h=0.2 (Fig. 3c). Similar conclusions may be drawn for the deeper water kh=2 (Fig. 4). However, in the case of the mass-transport velocity profile relatively far away from the wavemaker paddle (x=3h) the differences are twice as high as in the wave cases corresponding to kh=1. It should be noted that, for higher kh (Fig. 7), the results corresponding to two methods of averaging are more consistent; however the highest difference may still reach 20 % as in the case of longer waves (Fig. 5c).

Putting aside the problems with averaging of the velocity field, the general conclusion from the analysis of the results depicted in Figs. 27 is that admitting higher-order terms of the Taylor series expansions of the boundary conditions imposed at the moving boundaries leads to more accurate prediction of the time-independent velocity near the wavemaker as expected. This is of significant importance especially for the modelling of wave-induced mixing and the evolution of the water temperature field as the higher-order approach predicts enhanced streaming in layers adjacent to the bottom and the surface. This leads to more intense mixing processes and heat exchange.

3.3 Evolution of a temperature field

The modelling of wave-induced vertical mixing in terms of the temperature redistribution in the water column relies on the solution to the two-dimensional boundary-value problem defined by Eqs. (30)–(34). The input parameters being introduced to the governing advection–diffusion equation (Eq. 30) are the diffusion coefficient κv calculated using Eqs. (35) and (36) and the time-independent wave velocity field (U, W). The modelling procedure remains in accordance with previously published results (Sulisz and Paprota2019) with different velocity field input. The temperature T0(z) determines the initial thermal state of the fluid (Fig. 8).

Figure 8Initial temperature in the fluid domain. (a) Spatial distribution, (b) vertical profile.


In Figs. 911, the temperature spatial distributions representing the thermal states of the undulating water body after 100 s is provided for the test cases listed in Table 1. The results in the plots correspond to a weakly nonlinear solution and higher-order model predictions of the time-independent wave velocity field. Both methods are compared in order to highlight the effects of nonlinearity on vertical mixing due to waves.

Figure 9Changes in the temperature field due to waves after 100 s for kh=0.5, (a, b) H/h=0.05, (c, d) H/h=0.1 and (e, f) H/h=0.2. Weakly nonlinear theory (a, c, e; left), higher-order theory (b, d, f; right).


Figure 10Changes in the temperature field due to waves after 100 s for kh=1, (a, b) H/h=0.05, (c, d) H/h=0.1 and (e, f) H/h=0.2. Weakly nonlinear theory (a, c, e; left), higher-order theory (b, d, f; right).


Figure 11Changes in the temperature field due to waves after 100 s for kh=2, (a, b) H/h=0.05, (c, d) H/h=0.1 and (e, f) H/h=0.2. Weakly nonlinear theory (a, c, e; left), higher-order theory (b, d, f; right).


In the direct vicinity of the wavemaker paddle, the initial state of the water temperature is uniformly stratified according to T0. The moving wavemaker paddle generates regular waves, and the layers of water of equal temperature are deformed by the oscillatory motion of the water body according to the time-independent velocity field presented in Figs. 24. Advection plays an important role in wave-induced mixing processes as even small changes in the velocity strongly affect the evolution of the temperature field. This implies major differences between weakly nonlinear and higher-order predictions of the resultant temperature when the steepness of waves increases. This effect is clearly seen in Figs. 911. As previously stated, the discrepancy between weakly nonlinear velocities and the velocity predictions admitting higher-order terms in the free-surface boundary conditions grows with increasing steepness and wave nonlinearity (Ursell number) and may even reach 40 % in the subsurface layer.

Although the highest relative differences are relevant for the longest waves, the effect of including higher-order terms is most apparent for waves of the largest kh characterized by the highest time-independent velocities. It may be seen in Fig. 9 that after 100 s the colder water is moved to the surface close to the wavemaker only for the steepest waves (Fig. 9f) when predicted by the higher-order model, while for the corresponding weakly nonlinear temperature field this process is less intensive (Fig. 9e). Additionally, the higher-order results exhibit some variation in mixing along the flume (Fig. 9f).

In the case of waves characterized by kh=1, the wave-induced mixing is more intensive (see Table 1), and for the steepest waves (Fig. 10e, f) the warmer water is moved away from the wavemaker. It should be noted that the higher-order solution predicts enhanced streaming near the surface (Fig. 10f). In this way, the decreased water surface temperature affects regions located further away from the wavemaker when compared to weakly nonlinear results (Fig. 10). The similar temperature field modification affects the waves of kh=2, but for the lower analysed height (cf. Figs. 10e, f and 11c, d). In the case of the highest waves of kh=2, the warmer water is almost completely swept away from the region of direct wavemaker action (Fig. 11e, f).

3.4 Further discussions

Energy input from wind to ocean surface waves is tremendous and exceeds 60 TW (Wang and Huang2004). A large amount of energy dissipates, which implies stirring and mixing processes in the oceanic mixed layer. A parameterization scheme for surface wave-induced mixing was proposed, and numerical experiments indicate that this parameterization affects the performance and outcome of ocean circulation models (Qiao et al.2004, 2010; Xia et al.2006) as well as climate predictions (Song et al.2007; Huang et al.2008). Although field observations suggest that surface waves can generate vertical mixing (Matsuno et al.2006), it is still difficult to distinguish the mixing originating from the variety of processes that accompanied ocean waves. Although the problem of wave-induced vertical mixing is of significant importance for physical oceanographers and climatologists, the research on this subject is still in its infancy. Wave-induced mixing is a very important process from a practical point of view and is a challenging problem for theoretical investigations. The understanding of wave-induced mixing is of fundamental importance for the modelling and accurate prediction of ocean transport processes and climate changes. The problem is that the parameterization scheme applied in the modelling and prediction of ocean surface wave-induced mixing is based on drastic simplifications. Wave models derived using simplifying assumptions cannot be applied to predict ocean surface wave-induced mixing with sufficient accuracy. The present studies clearly show that it is necessary to apply at least a weakly nonlinear wave theory to obtain reasonable results (no wave-induced mass transport in linear theory means no heat exchange along the direction of wave propagation), while an increased accuracy is possible using advanced nonlinear wave models admitting higher-order effects. It should also be noted that in the present studies the effects of turbulence and the development of viscous boundary layers are neglected, while these processes may have some impact on larger timescales. For example, the mass-transport velocity profile changes due to the effects of viscosity near the boundaries (the undulating free surface and the bottom) for a sufficiently long time interval, which was studied theoretically by Longuet-Higgins (1953) and also investigated experimentally in a wave flume (see e.g. Swan1990; Grue and Kolaas2017).

There is one more critical outcome from the point of view of engineers and scientists working on the modelling of ocean transport processes, wave-induced mixing and climate changes. Namely, since it is difficult to distinguish the mixing that accompanied ocean waves, the only chance to provide reliable insight into the wave-induced mixing processes is to conduct laboratory experiments in a wave flume. The repeatable experiments in the well-controlled environment of a wave flume enable us to perform an accurate investigation that is essential in the analysis of the physics of the wave-induced mixing phenomenon. Moreover, laboratory investigations provide useful data for the analysis of the correlations between spectral and statistical characteristics of wave regimes and wave-induced mixing processes. Finally, laboratory experiments enable us to avoid various side effects and separate the mixing from other processes that accompanied ocean waves such as wave breaking. This is of fundamental importance for understanding of mixing and accurate calibration and verification of numerical models. The problem is that, in addition to progressive laboratory waves, the moving wavemaker represented by the kinematic wavemaker boundary condition enforces the return current that affects transport processes and wave-induced mixing.

It is the idea behind the present study, which employs a numerical wave flume model, to thoroughly analyse wave-induced mixing effects using the derived numerical approach and assist further works on experimental fluid mechanics aiming at better understanding of transport processes in the open ocean. To the best of the authors' knowledge, the derived wave-induced mixing approach is the only available numerical solution that, in addition to nonlinear free-surface boundary conditions, also satisfies the kinematic wavemaker boundary condition. Accordingly, the wavemaker model admits a return current and may be applied to quantify and separate the effects of the return flow on wave-induced mixing processes. The presented results should hopefully improve simple models based on the Stokes drift applicable to random ocean waves (Myrhaug et al.2018). As previously mentioned, the derived modelling framework may be modified to cover open-ocean conditions for the periodic domains and quasi-random sea states using other forms of wave excitation (Paprota2019).

4 Conclusions

The applicability of the wave-induced mixing model for waves generated in a wave flume is validated based on the solution admitting higher-order nonlinearities. In the range of wave conditions covering transitional and shallow waters, the weakly nonlinear results are in reasonable agreement with the more accurate pseudo-spectral solution in the case of waves of low to moderate steepness. The general discrepancy grows by increasing the wave height and the wavelength. Unlike the weakly nonlinear approach, the higher-order model is able to predict enhanced subsurface streaming affecting the evolution of the surface temperature for more severe sea states. It is due to the fact that the time-independent velocity field predicted by both methods differs especially in the subsurface and near-bottom layers of the oscillating water body.

General ocean circulation models admitting wave-induced vertical mixing but relying on simplified assumptions cannot predict input from mixing with sufficient accuracy. It is necessary to apply at least a weakly nonlinear correction to obtain reasonable approximation. For improved predictions, advanced highly nonlinear models are preferred, which is confirmed by the present study. Moreover, the derived method allows a return current to be correctly quantified in experimental investigations on wave-induced vertical mixing for better interpretation of laboratory results, giving more information for further improvements to parametrization schemes.

Data availability

Datasets for this research are available at (Paprota2023).

Author contributions

Conceptualization: MP and WS. Data curation: MP. Formal analysis: MP and WS. Investigation: MP and WS. Methodology: MP and WS. Software: MP and WS. Supervision: WS. Validation: MP and WS. Visualization: MP. Writing – original draft preparation: MP and WS. Writing – review and editing: MP and WS.

Competing interests

The contact author has declared that neither of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Review statement

This paper was edited by Ira Didenkulova and reviewed by two anonymous referees.


Babanin, A. V. and Haus, B. K.: On the existence of water turbulence induced by nonbreaking surface waves, J. Phys. Oceanogr., 39, 2675–2679, 2009. a

Chen, Y. Y., Hsu, H. C., and Hwung, H. H.: A Lagrangian asymptotic solution for finite-amplitude standing waves, Appl. Math. Comput., 215, 2891–2900, 2009. a

Dai, D., Qiao, F., Sulisz, W., Han, L., and Babanin, A. V.: An experiment on the non-breaking surface-wave-induced mixing, J. Phys. Oceanogr., 40, 2180–2188, 2010. a, b

Dean, R. G. and Dalrymple, R. A.: Water wave mechanics for engineers and scientists, vol. 2 of Advanced Series on Ocean Engineering, Prentice-Hall, Inc., Englewood Cliffs, 1984. a, b

Dommermuth, D. G. and Yue, D. K. P.: High-order spectral method for the study of nonlinear gravity waves, J. Fluid Mech., 184, 267–288, 1987. a

Grilli, S. T. and Svendsen, I. A.: Corner problems and global accuracy in the boundary element solution of nonlinear wave flows, Eng. Anal. Bound. Elem., 7, 178–195, 1990. a

Grue, J. and Kolaas, J.: Experimental particle paths and drift velocity in steep waves at finite water depth, J. Fluid Mech., 810,, 2017. a

Huang, C., Qiao, F., and Song, Z.: The effect of the wave-induced mixing on the upper ocean temperature in a climate model, Acta Oceanol. Sin., 27, 104–111, 2008. a

Hudspeth, R. T. and Sulisz, W.: Stokes drift in two-dimensional wave flumes, J. Fluid Mech., 230, 209–229, 1991. a, b, c, d, e, f, g, h

Longuet-Higgins, M. S.: Mass transport in water waves, Philos. T. R. Soc. A, 245, 535–581, 1953. a, b

Longuet-Higgins, M. S.: Eulerian and Lagrangian aspects of surface waves, J. Fluid Mech., 173, 683–707, 1986. a

Matsuno, T., Lee, J.-S., Shimizu, M., Kim, S.-H., and Pang, I.-C.: Measurements of the turbulent energy dissipation rate ε and an evaluation of the dispersion process of the Changjiang Diluted Water in the East China Sea, J. Geophys. Res., 111, C11S09,, 2006. a, b

McWilliams, J. C. and Sullivan, P. P.: Vertical mixing by Langmuir circulations, Spill Sci. Technol. B., 6, 225–237, 2000. a

Myrhaug, D., Wang, H., and Holmedal, L. E.: Stokes transport in layers in the water column based on long-term wind statistics, Oceanologia, 61, 522–526, 2018. a, b

Paprota, M.: On pressure forcing techniques for numerical modelling of unidirectional nonlinear water waves, Appl. Ocean Res., 90, 101857,, 2019. a, b, c

Paprota, M.: Experimental study on spatial variation of mass transport induced by surface waves generated in a finite-depth laboratory flume, J. Phys. Oceanogr., 50, 3501–3511, 2020. a, b, c, d, e

Paprota, M.: Wave-induced mixing in a numerical wave flume, Dryad [code, data set],, 2023. a, b

Paprota, M. and Sulisz, W.: Particle trajectories and mass transport under mechanically generated nonlinear water waves, Phys. Fluids, 30, 102101,, 2018. a, b, c, d, e, f, g, h, i, j, k, l

Paprota, M. and Sulisz, W.: Improving performance of a semi-analytical model for nonlinear water waves, J. Hydro-Environ. Res., 22, 38–49, 2019. a

Paprota, M., Sulisz, W., and Reda, A.: Experimental study of wave-induced mass transport, J. Hydraul. Res., 54, 423–434, 2016. a, b, c, d, e

Press, W. H., Flannery, B., Teukolsky, S. A., and Vetterling, W. T.: Numerical recipes in C, Cambridge University Press, Cambridge, UK, ISBN 0-521-43108-5, 1988. a, b

Qiao, F., Yuan, Y., Yang, Y., Zheng, Q., Xia, C., and Ma, J.: Wave-induced mixing in the upper ocean: distribution and application to a global ocean circulation model, Geophys. Res. Lett., 31, L11303,, 2004.  a, b, c

Qiao, F., Yuan, Y., Ezer, T., Xia, C., Yang, Y., Lu, X., and Song, Z.: A three-dimensional surface wave-ocean circulation coupled model and its initial testing, Ocean Dynam., 60, 1339–1355, 2010. a, b

Qiao, F., Song, Z., Bao, Y., Song, Y., Shu, Q., Huang, C., and Zhao, W.: Development and valuation of an Earth System Model with surface gravity waves, J. Geophys. Res., 118, 4514–4524, 2013. a

Shu, Q., Qiao, F., Song, Z., Xia, C., and Yang, Y.: Improvement of MOM4 by including surface wave-induced vertical mixing, Ocean Model., 40, 42–51, 2011. a

Song, Z., Qiao, F., Yang, Y., and Yuan, Y.: An improvement of the too cold tongue in the tropical Pacific with the development of an ocean-wave-atmosphere coupled numerical model, Prog. Nat. Sci., 17, 576–583, 2007. a, b

Sulisz, W. and Hudspeth, R. T.: Complete second-order solution for water waves generated in wave flumes, J. Fluids Struct., 7, 253–268, 1993. a

Sulisz, W. and Paprota, M.: Theoretical and experimental investigations of wave-induced vertical mixing, Math. Probl. Eng., 2015,, 2015. a, b, c, d, e, f, g

Sulisz, W. and Paprota, M.: On modeling of wave-induced vertical mixing, Ocean Eng., 194, 106622,, 2019. a, b, c, d, e, f, g, h, i

Swan, C.: Convection within an experimental wave flume, J. Hydraul. Res., 28, 273–282, 1990. a

Wang, W. and Huang, R. X.: Wind energy input to the surface waves, J. Phys. Oceanogr., 34, 1276–1280, 2004. a

West, B. J., Brueckner, K. A., and Janda, R.: A new numerical method for surface hydrodynamics, J. Geophys. Res., 92, 11803–11824, 1987. a

Xia, C., Qiao, F., Yang, Y., Ma, J., and Yuan, Y.: Three-dimensional structure of the summertime circulation in the Yellow Sea from a wave-tide-circulation coupled model, J. Geophys. Res., 111, C11S03,, 2006. a, b

Short summary
Ocean surface temperature is strongly affected by vertical mixing processes induced by waves. The correct estimation of a thermal state of the upper ocean layer is crucial for climate modelling and weather forecasting. The accuracy of predictions of exchange of heat between oceanic waters and the atmosphere depends on a reliable input from wave models admitting mixing processes. The study discusses the importance of higher-order approximation of vertical mixing due to high and steep waves.
Final-revised paper