Mixing Efficiency for Breaking Internal Solitary Waves

We propose a semi-analytic approach to estimate mixing induced by internal solitary waves (ISWs) breaking over a sloping boundary. The theoretical framework we develop describes the energetics of a stratified fluid flow during the interaction of ISWs with a slope. In particular, through both Ozmidov and Thorpe lengthscales we derive a heuristic expression for mixing efficiency, valid for plunging and plunging-collapsing breakers. On the other hand, we perform laboratory experiments that also provide mixing efficiency by estimating both changes of the background potential energy and the incident and reflected ISWs energy. We then compare those results with the ones derived by our theoretical model, obtaining a power law that relates the two quantities. This function allows to estimate mixing efficiency when the change in the background potential energy due to the ISWs breaking is unknown. We finally successfully apply our model to estimate the mixing efficiency under real field conditions.

parcels, entrained into the intermediate layer, are redistributed by chaotic advection and stretched in response to the large strain rates produced by turbulence (Peltier & Caulfield, 2003).
Mixing is usually quantified by the mixing efficiency, that is, the amount of kinetic energy irreversibly converted into background potential energy as well as the amount of incident wave energy contributing to irreversible mixing of the density field (Gregg et al., 2018;Osborn, 1980;Winters et al., 1995). Still today, semi-empirical parametrizations assume a constant mixing efficiency for the whole ocean, of about 0.2 (Nikurashin & Ferrari, 2013), although there is clear consensus from laboratory and direct numerical simulations that mixing efficiency is highly variable and changes with mechanism, evolutionary stage of a turbulent event, and location in the domain (Barry et al., 2001;Bouffard & Boegman, 2013;Ivey et al., 2008;Shih et al., 2005;Smyth et al., 2001). In particular, a variety of laboratory and DNS experiments (e.g., Barry et al. [2001] and Shih et al. [2005]; for a complete review, we refer the reader to Ivey et al. [2008] and Gregg et al. [2018]) investigated the variability of mixing efficiency as function of the Reynolds buoyancy number (Re b ), which represents the ratio between turbulent stirring and effects of buoyancy and viscosity. These studies showed a nonlinear dependency of mixing efficiency on Re b , reaching the maximum value of 0.2 for Re b of O(10 2 ). In the ocean, Re b varies from O(1) (e.g., Wuest & Lorke, 2005) to O(10 5 ), in deep regions (e.g., Ferron et al. [1998]) and, consequently, the use of a constant value for mixing efficiency is inaccurate. Furthermore, observations show that the oceans are not uniformly and steadily agitated, but there is a relatively quiescent interior and turbulent hot spots, often near boundaries (Artale et al., 2018;Wunsch & Ferrari, 2004). Thus, it is important to determine the geography of mixing hot spots and the dynamics responsible for their spatial and temporal distribution. In particular, mixing is enhanced through the water column in those regions that are characterized by rough bathymetry (Kokoszka et al., 2019;Ledwell et al., 2000;Mashayek et al., 2017;Polzin et al., 1997;St. Laurent et al., 2012). However, seafloor-induced turbulence is expected to occur within the bottom boundary layer. This suggests that internal waves are likely responsible for transport energy up from the bottom, due to their interaction with bathymetric features such as seamounts, sills, ridges and continental slopes (Cavaliere et al., 2021;Polzin et al., 1997).
Field observations estimated energetics of shoaling and breaking internal waves (e.g., Davis & Monismith, 2011;Klymak et al., 2008;Scotti et al., 2006;Shroyer et al., 2010;Walter et al., 2012). However, these results are limited by the spatial and temporal resolution of field data. For this, laboratory and numerical studies have been widely carried out by using small-scale idealized domains. Helfrich (1992) performed laboratory experiments on shoaling and breaking ISWs interacting with a sloping boundary, where the interface intersected the bottom and breaking occurred in all cases. These experiments allowed to measure the mixing efficiency, defined as the ratio between the change in potential energy near the breaking point (i.e., ΔE p ) and the net energy into the breaking region: where S w = A w /λ w , with A w wave amplitude and λ w = S/A w characteristic wavelength, with S wave surface. Boegman et al. (2005) modeled the reflectance coefficient R as a function of Ir, also attempting to relate the Iribarren number to the values for ϵ obtained by Michallet and Ivey (1999).
Broadening the parameter range that was explored in previous experiments, La Forgia, Adduce, and Falcini (2018) found that, for plunging breakers (Ir < 1), a dominant steepening of the trailing edge of the wave (Sutherland et al., 2013), followed by a quick clockwise overturning in the onshore direction, induces strong local mixing (Figure 1a). For collapsing breakers (Ir between 1 and 1.5), the trapped dense fluid leaves its original position with a fast downward motion in the adverse pressure gradient region, and a turbulent separated bolus forms and quickly dissipates ( Figure 1b). Then, part of the incident wave is reflected and a gravity current composed by the denser fluid flows up the slope, until hydrostatic conditions are reestablished (Aghsaee & Boegman, 2015;la Forgia, Tokyay, Adduce, & Constantinescu, 2020). An intermediate breaking mechanism, that is, the plunging-collapsing breaker, was also observed when the two main shoaling processes occur in the breaking location (i.e., for Ir ≃ 1). Finally, in the case of surging breakers (Ir > 1.5), ISWs are not subject to any observable large-scale instability during the shoaling, until the wave trough reaches the sloping bottom ( Figure 1c); the wave is almost reflected by the right wall of the tank and a gravity current composed of denser fluid moves up the slope causing mixing (la Forgia, Tokyay, Adduce, & Constantinescu, 2018). Each breaking mechanism is characterized by different effects in terms of mixing, fluid entrainment, and shear stress over the bottom (Aghsaee & Boegman, 2015;la Forgia, Tokyay, Adduce, & Constantinescu, 2020).
Numerical modeling has been used extensively to study shoaling and breaking internal waves on slopes both at the field Lamb, 2002;Vlasenko & Stashchuk, 2007;Walter et al., 2012) and laboratory scales (Aghsaee et al., 2010;Arthur & Fringer, 2014;Vlasenko & Hutter, 2002;Venayagamoorthy & Fringer, 2007). In particular,  revisited the laboratory experiments of Michallet and Ivey (1999) using two-dimensional numerical simulations and they proposed a parameterization for R as a function of the Iribarren number Ir for laboratory scale waves without side-wall effects: with ξ 0 = 0.78 ± 0.02, least squares fitting parameter. Aghsaee et al. (2010) investigated the breaking of fully nonlinear ISWs of depression, shoaling upon a uniformly sloping boundary in a smoothed two-layer density field by using high-resolution two-dimensional simulations. They proposed a parameterization for R similar to that of  by also considering the energy dissipations between the toe of the slope and the breaking location.
In general, to quantify mixing, all experimental and numerical studies we briefly summarized above analyze the evolution of density and velocity fields. Here, through the Ozmidov and Thorpe lengthscales, we develop a novel expression for the mixing efficiency, valid for plunging, and plunging-collapsing breakers, that is, the breakers type mostly expected in the continental shelf region (  Adduce, & Falcini, 2018). By means of experimental results, we relate the mixing efficiency derived by our heuristic model with the one obtained by a canonical, theoretical definition. This relation allows us to directly estimate the change in the background potential energy induced by each breaking event, as a function of the ISWs features and the inclination of the sloping boundary.

Theoretical Background
To mix a stably stratified fluid, energy is required to lift heavy fluid elements and lower light elements. For closed fluid systems, only irreversible, diabatic processes can change the probability density function (p.d.f.) of water column density (Winters et al., 1995). In this context, by defining as adiabatic a process in which there is no heat or molecular mass transfer, and as diabatic a process that is not adiabatic, Winters et al. (1995) partitioned the changes in potential energy due to diabatic mixing from changes due to adiabatic processes, in order to properly quantify the energetics of mixing. By assuming that the state of the flow is known within a fixed two-dimensional domain D, the instantaneous potential and kinetic energies of a fluid are where ρ is the local instantaneous density field, ρ 0 is a constant reference density, z is the vertical spatial coordinate, u, w are the velocities in the x, z directions, respectively. For adiabatic processes, changes in the fluid potential energy result from the switching of the kinetic energy into potential energy without any diffusive mixing (i.e., no heat or mass transfer occurs). Otherwise, diabatic processes produce the change of the total potential energy of the fluid, induced by irreversible molecular diffusion. The instantaneous volume-integrated background potential energy is defined as where ρ*(x, z, t) is the density field in the configuration of the minimum potential energy, obtained by sorting the fluid parcels by an adiabatic volume-conserving rearrangement of ρ. For this reason, it uniquely depends on the pdf of density and thus it is not affected by the instantaneous spatial distribution of density in the flow domain. Mixing of the density field ρ*(x, z, t) corresponds to a change in the p.d.f. induced by mass diffusion, leading to a reduction of the density variance. Blending of fluid parcels is enhanced by turbulence, which causes the steepening of the scalar gradients and the increasing of the iso-scalar surfaces (Winters & D'Asaro, 1996). Changes in the background potential energy are thus associated with the energy consumed in mixing the fluid and can be used to characterize this process.
The difference between the total potential energy and the background potential energy quantifies the amount of potential energy released in the adiabatic transition from ρ(x, z, t) to ρ*(x, z, t) without altering the p.d.f. of density. This amount of energy is called available potential energy (APE), since it represents the amount of potential energy stored in the fluid when it is not in gravitational equilibrium (Lorenz, 1955): The APE (E a ) represents the potential energy released if the fluid were to be adiabatically rearranged to the state of minimum potential energy (Thorpe, 1977). This approach estimates empirically the length scales of turbulent overturning in a stratified turbulent flow and is useful for the analysis of vertical density profiles when density inversions are the result of turbulent stirring. The method consists of rearranging or ordering an observed potential density profile, which may contain inversions, into a stable monotonic profile with no inversions. In practice, by considering n samples of density ρ n , each of which was observed at depth z n , the Thorpe displacement is d n = z m − z n , where the sample at depth z n had to be moved to depth z m to generate the stable profile (Dillon, 1982). The resultant Thorpe length scale is the root mean square value of the distances d n : L T can be estimated from fine-scale density profiles (Galbraith & Kelley, 1996;Gargett & Garner, 2008;Park et al., 2014), by using a reordering routine or a sorting algorithm that converts the observed profile into one in which density increases downwards everywhere (Thorpe, 2005).
Strictly related to L T is the Ozmidov length scale, defined as (Dillon, 1982;Ozmidov, 1965;Lumley, 1964;Thorpe, 2005) where ɛ is the turbulent dissipation rate and N is the buoyancy frequency. L O provides a measure of the vertical size of the largest eddies that may overturn in stably stratified water. The Ozmidov scale can be derived on dimensional grounds when viscosity is negligible (Dillon, 1982;Osborn, 1980;Ozmidov, 1965;Thorpe, 2005). Dillon (1982) suggested a linear relationship between L T and L O (Park et al., 2014): On the other hand, the Kolmogorov length scale, given by represents the scale at which viscosity dominates and the turbulent kinetic energy is dissipated into heat. Considering the standard assumptions for stratified turbulence (Davidson, 2015;Davidson et al., 2013;Gregg et al., 2018;Imberger & Boashash, 1986), we notice that Equation 11 relates L K to L O through the Reynolds number of the large turbulent motion Re = c w L O /ν, where c w is the ISW phase speed.
Thus, we can infer that the sizes of the largest and smallest eddies in high Reynolds number turbulence potentially differ by many orders of magnitude.
We finally note that, in order to eliminate the dissipation rate ɛ, by combining the Kolmogorov length scale in Equation 11 and the Ozmidov scale (Equation 9), we can write: Therefore, we write the Ozmidov scale as

Heuristic Model for Mixing Efficiency Estimation
The potential energy change ΔE p in Equation 1 represents the increase in the irreversible potential energy caused by wave breaking on the slope . Hereafter, we indicate it as ΔE b , that is, the change in the background potential energy.
We calculate the wave APE considering finite and infinite boundaries of the domain for the case of ISWs in two-layer fluid (see also Blokhina [2009]). We assume that the wave of amplitude A w is propagating in a two-layer fluid with depths h 1 and h 2 , and densities ρ 1 and ρ 2 , with ρ 2 > ρ 1 (Figure 2), with a symmetric structure with respect to the wave; the horizontal domain along which the ISW propagates is 2l 0 (corresponding to l t − x 0 in Figure 2). We place the coordinate system at the interface between the two layers, in order to set the center of mass of each layer at half-depth. It results that the background potential energy (Equation 6) is The potential energy E p of the ISW calculated by Equation 4 is where η(x) is the ISW profile. Thus, the APE of the wave is Michallet and Ivey (1999) where t 0 and t 1 are appropriate times, chosen to compute the wave energy. Equation 17 can be rewritten through a change of variable as where we used relation (Equation 16). We notice in Equation 18 that the integration domain extents to infinity, in order to reach a true APE value (see Hebert [1988] and Lamb [2008] for details).
The APE of an overturn can be defined as  2 2 a T E N L (Imberger & Boashash, 1986;Ivey & Imberger, 1991). Given the relationship (10) between L T and L O , we find that E a is highly correlated with the generation of the largest eddies that may overturn in stably stratified water. Therefore, by assuming that most of the incident wave energy is expended to generate an overturn, that is,    (a) Sketch of the initial experimental setting, and density profiles associated to both the lock region, on the left, and the fluid ambient, on the right. The Perspex wave tank is 3 m long, 0.3 m high, and 0.2 m wide. The lower layer is filled with a solution of sodium chloride (NaCl) with density ρ 2 , and the upper layer with fresh water of uniform density ρ 1 < ρ 2 . the fresh water is dyed by a controlled quantity (0.3 ml/l) of Methylthioninium chloride (C 16 H 18 ClN 3 S). Density of the saline mixture is measured by a density meter (Anton Paar DMA 4100M), with an accuracy of 10 −1 kg/ m 3 . A charge-coupled device (CCD) camera placed at a fixed distance from the front wall of the tank is used to record the flow evolution with a frequency of 25 Hz and a spatial resolution of 1,024 × 668 pixels. The resolution of each pixel is approximately 3 × 3 mm.
where we introduced the wave mass per unit of width, M w ≃ A w λ w ρ 1 , written in terms of simple geometrical quantities, in order to be dimensionally consistent. In Equation 19, R can be parameterized through Equation 3, as done by  and Aghsaee et al. (2010).
) still requires the knowledge of ΔE b , which we will be evaluated in the next section.We remark that ϵ O , although dimensionally consistent with ϵ, cannot be rigorously considered as a mixing efficiency. Indeed, it is not obtained in the canonical fashion by estimating the ISW total mechanical energy (e.g., as in Davies Wykes and Dalziel [2014]), but as a function of the main wave geometrical parameters.

Experimental and Numerical Assessment
To diagnose Equation 19 we performed laboratory experiments, generating ISWs by the standard lock-release method in a two-layer stratified system (Kao et al., 1985;La Forgia, Adduce, & Falcini, 2018). Initially, at a distance x 0 from the left wall of the tank, a 4 mm thick, vertical and removable Perspex gate separates the lock fluid, on the left hand-side of the tank, from the ambient fluid on the right (l t − x 0 in Figure 2). For each side of the tank, a stratified two-layer distribution is obtained by filling the lower layer with a solution of density ρ 2 and depth h 2 , and the upper layer with fresh water of uniform density ρ 1 < ρ 2 and depth h 1 < h 2 . This particular setting is aimed to reproduce ISWs of depression, the most common type of ISWs detected in the coastal areas (Alpers et al., 2008;Duda et al., 2004;Helfrich & Melville, 2006;Jackson et al., 2013;Osborne & Burch, 1980), where the stratification shows a relatively thin upper layer. Before each run, we measure the density of the saline mixture in order to make sure that the density difference between the two uniform layers is of 30 kg/m 3 . Indeed, in the coastal ocean, the typical value of the Boussinesq parameter, that is, the ratio of the density difference between the two layers and a reference density, is about 0.03 (Chen et al., 2007). We observe that, for this value, the effects of the disturbances at the free surface are negligible, and do not have significant impact on the induced mixing (Kodaira et al., 2016;la Forgia & Sciortino, 2019). The different depths of the pycnoclines of the two stratifications produces the interfacial displacement η P ; a sloping boundary, making an angle θ with the horizontal, is placed on the right hand-side of the tank to induce the ISWs breaking and partial reflection ( Figure 2). All experiments start with the gate removal and the consequent gravity collapse: an ISW of depression generates and propagates downstream.
We perform seven runs, characterized by different initial conditions in order to generate ISWs having different geometric and kinematic features (see Table 1 ca, 2016). The evaluation of the instantaneous density fields has an accuracy of about 0.1 kg/m 3 (i.e., the ratio between the maximum range of density and the gray scale levels).
The normalized density ρ*(x, z, t) is given by: where x and z represent the streamwise and vertical directions (Figure 2), t is the time since the gate removal, and ρ* ranges between 0 (fresh upper layer) and 1 (salty lower layer).
The analysis of the density fields allows us to infer the pycnocline thickness δ P (Figure 2) and the instantaneous pycnocline position that is associated to the iso-density level ρ*(x, z, t) = 0.5. For all cases, we derive an initial pycnocline thickness δ Pi of 1 cm. However, in general, the pycnocline thickness does not affect the ISW geometrical features and the breaking mechanism (la Forgia, Ottolenghi, et al., 2020;Sutherland et al., 2013). Then, we derive the ISW amplitude A w , surface S w , and celerity c w (i.e., the first derivative of the trough's position). The characteristic wavelength λ w is estimated as (Michallet & Ivey, 1999): where η(x) is the pycnocline displacement along the domain. Being L = 3 m the length of the tank, H = h 1 + h 2 the total water depth, g′ = g(ρ 2 − ρ 1 )/ρ the reduced gravity, and   b u g H the buoyancy velocity, we herein consider the normalized quantities x* = x/L, z* = z/H, and t* = t H/u b to analyze the spatiotemporal evolution of the performed experiments in a dimensionless form.
The gravity collapse, induced by the gate removal, causes the generation of an ISW of depression ( Figure 3) that propagates with approximately constant celerity, amplitude, and wavelength (as proved by the white strip at 1.5 < t* < 3.5 in Figure 4a for Case 5). The incident ISW approaches the sloping boundary, where it shoals and breaks (Figure 3). The wave breaking develops with the trailing edge overturning, which induces mixing and partial reflection: a train of four, rank-ordered ISWs are indeed observed to propagate upstream, toward the lock region (Figure 4c for Case 5). The Hovmöller diagram shows that the reflected ISWs celerity is proportional to their size (see Figure 4a at 5 < t* < 9). For all experiments, the train of reflected waves is composed by four ISWs.
To estimate the energy released at the sloping region, and thus to obtain the change of background potential energy (ΔE b ) in Equation 19, we evaluate the difference between the incident and the reflected ISWs energy. To avoid that experimental inaccuracies (e.g., camera resolution) may affect the evaluation of the instantaneous pycnocline position, and thus of wave geometrical features and ISW energies, we derive the ρ* = 0.5 iso-density line from the fully nonlinear and strongly dispersive DJL model (Dubreil-Jacotin, 1937;Dunphy et al., 2011;Long, 1953;Turkington et al., 1991;Xu & Stastna, 2019). This synergy between experimental runs and the DJL model allows us to obtain the instantaneous pycnocline position also for smaller waves, characterized by pycnocline displacements comparable to the pixel size.   Table 1).
Solutions of Equation 22a represents exact solutions of Euler's equations, and they are derived by a generalization of a variational technique and a numerical algorithm (Dunphy et al., 2011;Turkington et al., 1991). In this algorithm neither the wave amplitude, nor the wave propagation speed are specified. Whereas, the kinetic energy of the disturbance is minimized under the constraint that the available potential energy (scaled by ρ 0 gH) is held fixed (Dunphy et al., 2011;Xu & Stastna, 2019), with  density profile scaled by the reference density ρ 0 .
For each ISW, we set the DJL model by imposing the density profile measure in the tank experiment, and by tuning the wave APE in order to obtain, numerically, the experimental ISW amplitudes (e.g., see red line in Figure 4b for Case 5 incident ISW and red and yellow lines in Figure 4c for the first two reflected ISWs). From the DJL solutions for the pycnocline displacements, we evaluate ISWs energies through Equation 18. The total reflected ISWs energy E R , associated to the train of n reflected waves, is obtained as the sum of the single energetic contribution of each ISW, that is, E R = ∑ i E Ri , with i = 1, …, n.
Once derived the energy released over the breaking location, we estimate the effective amount of this energy lost for mixing processes. To quantify mixing induced by ISWs breaking, we adopt the energy budget method of Winters et al. (1995), which has been widely used in literature (Fragoso et al., 2013  the change in the background potential energy (ΔE b ) between the initial condition, associated to the ISW that is approaching the sloping boundary, and a final condition associated to the release of all the reflected waves. This allows us to exclude mixing induced by the generation mechanism, and to consider the effective change in the background potential energy produced by the breaking process.

Estimation of the Mixing Efficiency
For each experiment, we estimate both mixing efficiency that is obtained from laboratory experiments (ϵ), that is, directly from Equation 1, and the one we derive from the heuristic model (ϵ O ) in Equation 19. In particular, we estimate the ISWs features, that is, A w , M w , c w , and λ w , and their energies E 0 and E R through image analysis of the laboratory experiments and the DJL model. The change of the background potential energy (ΔE b ) is estimated from the resulting density fields by using image analysis (Winters et al., 1995).
To validate and extend our results to a larger number of cases, we consider also experimental results obtained by Michallet and Ivey (1999), who followed a laboratory technique similar to the one we adopt for our cases. These authors used two probes, placed in the wave tank to measure the density profiles nearby the sloped surface. Then they evaluated the change in the total potential energy (ΔE p ) occurring in the trapezoidal region upstream the sloping boundary. They compared, in particular, the total potential energy before and after the breaking event, such as the fluid confined within the considered domain is at rest, and the total available potential energy is zero (i.e., ΔE b = ΔE p ). Considered the non-uniform density distribution along the x-direction, we argue that the estimation of E b from experimental density fields provides larger values of potential energy with respect to those obtained by local probe measurements.
Despite these methodological differences, the reflectance parameter R = E R /E 0 obtained by our experimental and numerical analysis is related with the internal Iribarren number, and it follows the same trend discussed by Michallet and Ivey (1999) (Figure 5b). tion of the leading off-shore propagating ISW, assuming that the energy related to the following waves is negligible. The large coefficient of correlation ( However, in agreement with the results obtained by Michallet and Ivey (1999) and Arthur and Fringer (2014), there is not a feasible relation between Ir and the experimental mixing efficiency (Figure 5a). The internal Iribarren number is, indeed, a dimensionless parameter able to identify the different breaking mechanisms that may occur when an ISW shoals upon a sloping boundary, but it is not suitable to identify the amount of mixing induced by wave breaking. Parameters that define Ir only consider ISWs and sloping boundary main geometrical properties (i.e., A w , λ w , and s) and they do not account for either any energetic parameter or stratification effects. Our results confirm that mixing efficiency is not solely related to the Iribarren number. This means that mixing cannot be uniquely predicted from wave geometry but also ISW energetics must be considered.
In seeking to define a predictive tool for mixing efficiency, we compare results from laboratory experiments to those we derived from heuristic model (Figure 5a). Since our model takes into account only plunging breakers, besides all our cases (black solid dots in Figure 5), we consider a selected group of runs performed by Michallet and Ivey (1999), characterized by Ir < 0.85 (gray solid dots in Figure 5). The resulting relation between ϵ and ϵ O can be described by the power law which confirms that our parametrization can be considered a valid approximation to describe mixing efficiency induced by plunging and plunging-collapsing breakers.
Data associated to larger ϵ show a larger scattering that could be associated to the different approach adopted by Michallet and Ivey (1999) to estimate ΔE b and to evaluate the reflected energy.
Interestingly, although we obtain the power law from data that are characterized by Ir < 0.85, also plunging-collapsing breakers (i.e., the ones with Ir > 0.85) show a relatively good agreement with Equation 25. This suggests that the bulk mixing is largely affected by the overturning of the ISWs trailing edge, rather than by the boundary layer separation occurring close to the sloping bottom, associated to the collapsing breakers. For this reason, our results can be considered valid for a slightly wider range of Ir values (i.e., for Ir → 1.5).
The relation between ϵ and ϵ O allows us to directly derive the change in the background potential energy once the undisturbed stratification, ISWs geometry, and topographic features (Ir, M w , and E 0 ) are known: By estimating ΔE b , it is therefore possible to derive the expected mixing efficiency as ϵ = ΔE b /[E 0 (1 − R)] without directly evaluating the evolution of the density field during the breaking event. Although the variation of pycnocline thicknesses does not affect the breaking mechanism (i.e., ISWs geometrical features), we observe that for weaker stratifications, slower waves are expected. The wave celerity, indeed, represents a scaling factor for the wave energy (see Equation 17) and, consequently, for the mixing efficiency (see Equation 26). This means that for faster waves, lower values for the mixing efficiency will occur. Our heuristic model takes into account stratification effects in terms of wave celerity in the definition of the pseudo-mixing efficiency (Equation 19). We thus can consider the pseudo-mixing efficiency as representative of all possible density distributions characterized by layer depths h 1 < h 2 (i.e., for ISWs of depression).
Internal waves over topographic constraints enhance fine-scale shear and strain over rough bathymetry (Artale et al., 2018;Polzin et al., 1997). In particular, internal waves breaking induces an expansion of energy that results in mixing the stratified fluid and, consequently, causes upwelling of abyssal waters along sloping boundaries (Ferrari et al., 2016;Wunsch & Ferrari, 2004). We note that our estimation of mixing efficiency, based on the knowledge of both ISW and topographic geometries a priori, expands the ability to parametrize fine-scale shear and strain processes that are crucial for the diagnosis of diapycnal mixing and, in general, ocean heat content conundrum (Artale et al., 2018;McDougall & Ferrari, 2017).
Fields observations show strong evidence of large amplitude ISWs propagating in two-layer configurations, interacting with a sloped seafloor. To compare our results with those obtained in previous studies we finally test our model for some real-field cases. In particular, we consider ISWs breaking occurring in the following locations: Seneca Lake (Hunkins & Fliegel, 1973), State of New York, US; Lock Ness (Thorpe, 1974), Scottish Highlands; Sulu Sea (Apel et al., 1985), southwestern area of the Philippines. For these cases, Michallet and Ivey (1999) provide values of mixing efficiency associated to individual breaking events, using physical parameters described in previous literature. In particular, they obtained the following values for ϵ: 0.18 (Seneca Lake), 0.13 (Loch Ness), 0.10 (Sulu Sea). We test our model on these cases, evaluating ϵ = ΔE b / (E 0 (1 − R)), with ΔE b as in Equation 26, and E 0 calculated from Equation 18 with η(x) derived from the DJL model. We obtain a mixing efficiency of 0.17, 0.08, and 0.03, respectively. This suggests that our heuristic model reproduce a similar trend although underestimating the mixing efficiency with respect to the results obtained by Michallet and Ivey (1999). Our model could represent a promising tool to evaluate the mixing efficiency trend, although a rigorous comparison with real-field observations, aimed at deriving more precise values for ϵ, should be carried out. However, we expect our approach to conveniently capture the variability in mixing efficiency for ISWs with different geometrical features. In particular, this can be useful for studying the mixing effects of ISWs with different features breaking in a single specifying location.

Conclusions
For stratified flows, prediction of changes in kinetic and background potential energies represents the main challenge for estimating mixing, without directly evaluate instantaneous density and velocity fields by experimental studies, numerical modeling or real field observations.
In the present paper, we develop a synergy of laboratory experiments, numerical modeling and theoretical analysis, deriving a predictive tool for the mixing efficiency induced by breaking ISWs over a sloping boundary. The different ISWs breaking mechanisms are usually identified by the internal Iribarren number. However, it is not possible to derive a relation between waves and topography geometrical features and the associated mixing. Our theoretical approach is based on the available potential energy associated to the ISW and on specific turbulent scales, that is, Ozmidov and Thorpe lengthscales. Starting from the definition provided by Helfrich (1992) and Michallet and Ivey (1999), we obtain a heuristic expression for the mixing efficiency (pseudo-mixing efficiency, ϵ O ), assuming that all the incident ISW energy is involved in the overturning event. This is the expected breaking behavior for plunging and plunging-collapsing breakers, in response to the verticalization of the trailing wave edge as the ISW shoals over the sloping boundary. The relationship shows that the pseudo-mixing efficiency depends on both the change of the background potential energy (ΔE b ) and on the net wave energy (i.e., E 0 (1 − R)) involved in the mixing process. In particular, the incident energy E 0 appears to be well represented by the product of the wave mass and the squared wave celerity. We estimate mixing by both the canonical and the heuristic expressions for a set of laboratory experiments performed in order to widen the range of Iribarren number investigated by Michallet and Ivey (1999). The relation between ϵ and ϵ O is well approximated by a power law, with a good degree of correlation. This relationship allows us to directly estimate the change in the background potential energy, without a direct evaluation of the density field evolution. Indeed, it is possible to derive ΔE b once the undisturbed stratification, ISWs geometry, and topographic features are known (Ir, M w , and E 0 ). By testing our results for real field conditions, we found that our approach provides values of mixing efficiency close to those previously estimated for oceanographic cases, and, in particular, it captures the same trend. Therefore, the developed model represents a suitable predicting tool for estimation of ISWs-driven diapycnal mixing in the coastal oceans, expanding the ability to parametrize fine-scale shear and strain processes (Artale et al., 2018;McDougall & Ferrari, 2017).