1. Introduction
Heat and mass transport dynamics in the upper ocean, due to the presence of a boundary layer, are of practical importance in a range of applications. Examples include the initial stages of an oil spill, the motion of ocean biofilms, the transport of microplastic particles, and the movement of grease ice. These phenomena influence the properties of the free-surface boundary layer and its hydrodynamic behaviour, with direct consequences for wave-induced mass and heat transfer. Therefore, a better understanding of the complex thermal structure of the surface boundary layer may provide valuable insights into subsurface temperature patterns and their impact on the ocean environment.
It is well established that free-surface waves have important consequences for air–sea heat exchange processes. Witting (Reference Witting1971) and O’Brien (Reference O’Brien1967) analysed the effects of irrotational progressive oscillations and Gerstner waves, and demonstrated that air–ocean heat flux is sensitive to wave steepness. Later, Szeri (Reference Szeri2017) focused on mass and heat exchanges across a turbulent liquid–gas interface, whereas Hetsroni et al. (Reference Hetsroni, Mosyak and Yarin1997) showed that surface waves can have a significant effect on natural convection. Observations by Veron et al. (Reference Veron, Melville and Lenain2008) have shown that ocean waves could be very important in facilitating heat transport. Recently, Michele et al. (Reference Michele, Stuhlmeier and Borthwick2021, Reference Michele, Borthwick and van den Bremer2023) investigated the temperature field generated by heat sources beneath propagating and standing waves, and demonstrated that free-surface gravity waves play a crucial role in enhancing seabed heat transfer.
Thermal conduction in boundary layers is an established research field, with classical theoretical results and practical applications described by Landau & Lifshitz (Reference Landau and Lifshitz1989) and White (Reference White1991). Additional properties of heat transfer in a laminar boundary layer are thoroughly discussed by Lighthill (Reference Lighthill1950), whereas applications of heat and mass transfer driven by oscillatory flows are described by Pedley (Reference Pedley1972). However, to the authors’ knowledge, studies to date have been restricted to simplified irrotational or Gerstner wave velocity profiles that do not replicate the wave dynamics in the free-surface viscous boundary layer. Furthermore, these studies do not take into account the effect of the free-surface boundary layer on air–sea heat transfer or do not consider additional effects resulting from spatially dependent eddy viscosity, thermometric conductivity or heat sources profiles.
In this paper, we develop a mathematical theory to analyse heat transfer through the oceanic wave-driven free-surface boundary layer. Our weakly nonlinear model is based on solving the convection–diffusion equation for the temperature field in an incompressible fluid, where the velocity field is decoupled from temperature. We consider long-crested, two-dimensional surface gravity waves of amplitude greater than or comparable to the thickness of the boundary layer. Therefore, the Taylor series expansion of boundary conditions at the free surface is consistent only when using curvilinear coordinates (Longuet-Higgins Reference Longuet-Higgins1953) or alternatively, Lagrangian coordinates (Pierson Reference Pierson1962; Ünlüata & Mei Reference Ünlüata and Mei1970; Monin & Yaglom Reference Monin and Yaglom1971; Piedra-Cueva Reference Piedra-Cueva1995; Salmon Reference Salmon1998, Reference Salmon2020; Ng Reference Ng2004; Weber & Christensen Reference Weber and Christensen2019; Pizzo et al. Reference Pizzo, Lenain, Rømcke, Ellingsen and Smeltzer2023). Adopting a Lagrangian approach, we first derive a new form of Lagrangian velocity that is valid across the entire water depth, incorporating eddy viscosity profiles, and then use this to solve the convection–diffusion equation for temperature within the free-surface thermal boundary layer.
2. Governing equations
Consider two-dimensional water waves of amplitude
$A$
and frequency
$\omega$
propagating at the top of a viscous fluid domain, where
$x$
is horizontal distance from a fixed origin,
$z$
is distance vertically upwards from the undisturbed free surface, and
$t$
is time. The fluid has constant depth
$h$
, and density
$\rho$
, whereas the eddy kinematic viscosity
$\nu$
and thermometric conductivity
$\chi$
are assumed to be dependent on space and time. We further assume
$h$
to be comparable with the wavelength
$\lambda$
and much larger than the turbulent boundary layer thickness
$\delta$
, so we expect the effects due to viscosity and turbulence to be localised in regions close to the free surface and seabed, whereas in the core region, the fluid motion is dominated by irrotational potential flow (Mei et al. Reference Mei, Stiassnie and Yue2005). Constant fluid properties are reasonable solely for small variations in temperature
$\mathcal {T}\sim O(10)$
$^{\text {o}}$
C, and when the initial fluid temperature
$\mathcal {T}_i$
is far from freezing or boiling points. The resulting continuity equation, momentum equations and equation for convection and diffusion of relative temperature
$T=\mathcal {T}-\mathcal {T}_i$
in Eulerian form are (Landau & Lifshitz Reference Landau and Lifshitz1989; White Reference White1991)
\begin{align} \frac {\partial u}{\partial x}+\frac {\partial w}{\partial z}=0, \qquad\qquad\qquad\qquad\qquad\qquad\end{align}
\begin{align} \frac {\partial u}{\partial t}+ u\frac {\partial u}{\partial x}+ w\frac {\partial u}{\partial z}=-\frac {1}{\rho }\frac {\partial P}{\partial x}+\frac {\partial }{\partial x}\left (2\nu \frac {\partial u}{\partial x}\right )+\frac {\partial }{\partial z}\left [\nu \left (\frac {\partial w}{\partial x}+\frac {\partial u}{\partial z}\right )\right ], \quad\end{align}
\begin{align} \frac {\partial w}{\partial t}+ u\frac {\partial w}{\partial x}+ w\frac {\partial w}{\partial z}=-g-\frac {1}{\rho }\frac {\partial P}{\partial z}+\frac {\partial }{\partial z}\left (2\nu \frac {\partial w}{\partial z}\right )+\frac {\partial }{\partial x}\left [\nu \left (\frac {\partial w}{\partial x}+\frac {\partial u}{\partial z}\right )\right ], \end{align}
\begin{align} \frac {\partial T}{\partial t}+u\frac {\partial T}{\partial x}+ w\frac {\partial T}{\partial z}=\frac {\partial }{\partial x}\left (\chi \frac {\partial T}{\partial x}\right )+\frac {\partial }{\partial z}\left (\chi \frac {\partial T}{\partial z}\right ), \qquad\qquad\end{align}
where
$(u,w)$
are horizontal and vertical components of fluid velocity,
$g$
is the acceleration due to gravity, and
$P$
is pressure. Given that free-surface oscillations are much larger than
$\delta$
, standard Taylor series expansions of the free-surface boundary conditions about the horizontal plane
$z=0$
cannot be applied here. A convenient alternative approach is to adopt curvilinear coordinates (Longuet-Higgins Reference Longuet-Higgins1953) or Lagrangian coordinates, as in the case of the mass transport phenomena analysed by Ünlüata & Mei (Reference Ünlüata and Mei1970). Here, we use a Lagrangian approach. Let
$a$
and
$b$
be the initial horizontal and vertical coordinates of a fluid particle, where
$b=0$
denotes the free surface and
$b=-h$
the horizontal seabed. The governing equations (2.1)–(2.4) expressed in Lagrangian form are given by (Monin & Yaglom Reference Monin and Yaglom1971)
\begin{align} &\qquad\qquad\qquad\qquad\qquad\qquad\left [x,z\right ]\equiv \begin{vmatrix} \dfrac {\partial x}{\partial a} \dfrac {\partial x}{\partial b}\\[8pt] \dfrac {\partial z}{\partial a} \dfrac {\partial z}{\partial b} \end{vmatrix}=1, \end{align}
\begin{align} &\quad\frac {\partial ^2 x}{\partial t^2}=-\frac {1}{\rho }\left [P,z\right ]+\nu\, \nabla ^2\left (\frac {\partial x}{\partial t}\right )+2\left [\nu ,z\right ]\left [\frac {\partial x}{\partial t},z\right ]+\left [x,\nu \right ]\left (\left [\frac {\partial z}{\partial t},z\right ]+\left [x,\frac {\partial x}{\partial t}\right ]\right ), \end{align}
\begin{align} &\frac {\partial ^2 z}{\partial t^2}+g=-\frac {1}{\rho }\left [x,P\right ]+\nu\, \nabla ^2\left (\frac {\partial z}{\partial t}\right )+2\left [x,\nu \right ]\left [x,\frac {\partial z}{\partial t}\right ]+\left [\nu ,z\right ]\left (\left [\frac {\partial z}{\partial t},z\right ]+\left [x,\frac {\partial x}{\partial t}\right ]\right ), \end{align}
\begin{align} &\qquad\qquad\qquad\qquad\frac {\partial T}{\partial t}=\chi\, \nabla ^2T+\left [\chi ,z\right ]\left [T,z\right ]+\left [x,\chi \right ]\left [x,T\right ], \end{align}
where
$(a,b)$
are fixed as the fluid particle moves from place to place (Salmon Reference Salmon1998, Reference Salmon2020), and the Lagrangian operator
$\nabla ^2$
reads
\begin{equation} \nabla ^2f=\left [\left [f,z\right ],z\right ]+\left [x,\left [x,f\right ]\right ]. \end{equation}
We assign zero stress at the free surface
$b=0$
, no motion at the horizontal seabed
$b=-h$
, and a prescribed free-surface temperature
$T_s$
representing a heat source. We assume that the water temperature at depths greater than the boundary layer thickness
$z/\delta \rightarrow -\infty$
is fixed at
$\mathcal {T}_i$
for all time. The tangential and normal components of the stress tensor with respect to a material curve
$(x(t),z(t))$
are (Piedra-Cueva Reference Piedra-Cueva1995)
\begin{align} \boldsymbol {T_t}=\frac {\displaystyle\frac {\partial x}{\partial a}\frac {\partial z}{\partial a}\left (\tau _{zz}-\tau _{xx}\right )+\left [\left (\frac {\partial x}{\partial a}\right )^2-\left (\frac {\partial z}{\partial a}\right )^2\right ]\tau _{xz}}{|\boldsymbol {n}|^2}, \qquad\qquad\end{align}
\begin{align} \boldsymbol {T_n}=\frac {\displaystyle\left (\frac {\partial z}{\partial a}\right )^2\tau _{xx}+\left (\frac {\partial x}{\partial a}\right )^2\tau _{zz}-2\frac {\partial x}{\partial a}\frac {\partial z}{\partial a}\tau _{xz}-\left [\left (\frac {\partial x}{\partial a}\right )^2+\left (\frac {\partial z}{\partial a}\right )^2\right ]P}{|\boldsymbol {n}|^2}, \end{align}
where the normal and tangential vectors to the curve are given by
$ \boldsymbol {n}= (-\partial z/\partial a,\partial x/\partial a )$
,
$\boldsymbol {t}= (\partial x/\partial a,\partial z/\partial a )$
, and the stress components in Lagrangian form read (Kinsman Reference Kinsman1984)
\begin{equation} \tau _{xx}=2\nu \rho \left [\frac {\partial x}{\partial t},z\right ],\quad \tau _{zz}=2\nu \rho \left [x,\frac {\partial z}{\partial t}\right ],\quad \tau _{xz}=\nu \rho \left (\left [\frac {\partial z}{\partial t},z\right ]-\left [\frac {\partial x}{\partial t},x\right ]\right ). \end{equation}
Incompressible fluid motion is decoupled from temperature, so we first focus on the derivation of the oscillatory and wave-averaged mass transport velocity in Lagrangian form, and then use these solutions to solve the temperature problem.
3. Lagrangian velocity field
We adopt the following perturbation expansion up to second order
$O(A^2k^2)$
:
\begin{equation} \left \{x,z,P\right \}=\left \{a,b,P_0\right \}+\left \{x_1,z_1,P_1\right \}+\left \{x_2,z_2,P_2\right \}, \end{equation}
where
$\{x_1,z_1,P_1\}\sim \textit {O}(\epsilon )$
,
$\{x_2,z_2,P_2\}\sim \textit {O}(\epsilon ^2)$
are the first- and second-order components,
$\epsilon =Ak\ll 1$
represents small wave steepness, and
$k$
is the wavenumber (see Pierson Reference Pierson1962). Substituting the expansion (3.1) into the governing equations and boundary conditions derived in the previous section, we obtain a sequence of linearised boundary value problems up to order
$O(A^2k^2)$
. We now focus on the analytical solution of the velocity components and fluid pressure at each order.
3.1. Zeroth-order solution
The continuity equation (2.5) at zeroth order is identically satisfied, whereas the momentum equations (2.6)–(2.7) yield
\begin{align} \left [P_0,b\right ]=\frac {\partial P_0}{\partial a}=0, \qquad\end{align}
\begin{align} g+\frac {1}{\rho }\left [a,P_0\right ]=g+\frac {1}{\rho }\frac {\partial P_0}{\partial b}=0. \end{align}
Since we assumed zero stress at the free surface
$b=0$
, the foregoing equations give the straightforward hydrostatic solution in Lagrangian coordinates
$P_0=-\rho g b$
.
3.2. First-order solution
Solutions at first order are known for constant kinematic viscosity (Ünlüata & Mei Reference Ünlüata and Mei1970), whereas velocity components due to variable
$\nu$
are not readily available. In what follows, we use a procedure based on boundary layer theory to derive
$x_1,z_1$
.
The continuity and momentum equations at first order
$O(Ak)$
read
\begin{align} \frac {\partial x_1}{\partial a}+\frac {\partial z_1}{\partial b}=0, \qquad\qquad\qquad\qquad\qquad\qquad\end{align}
\begin{align} \frac {\partial ^2 x_1}{\partial t^2}=-g\frac {\partial z_1}{\partial a}-\frac {1}{\rho }\frac {\partial P_1}{\partial a}+\nu\, \nabla _L^2\left (\frac {\partial x_1}{\partial t}\right )+2\frac {\partial \nu }{\partial a}\frac {\partial ^2 x_1}{\partial t\, \partial a}+\frac {\partial \nu }{\partial b}\left (\frac {\partial ^2 z_1}{\partial t\, \partial a}+\frac {\partial ^2 x_1}{\partial t\, \partial b}\right ), \end{align}
\begin{align} \frac {\partial ^2 z_1}{\partial t^2}=-g\frac {\partial z_1}{\partial b}-\frac {1}{\rho }\frac {\partial P_1}{\partial b}+\nu\, \nabla _L^2\left (\frac {\partial z_1}{\partial t}\right )+2\frac {\partial \nu }{\partial b}\frac {\partial ^2 z_1}{\partial t\, \partial b}+\frac {\partial \nu }{\partial a}\left (\frac {\partial ^2 z_1}{\partial t\, \partial a}+\frac {\partial ^2 x_1}{\partial t\, \partial b}\right ), \end{align}
where the operator is
$\nabla ^2_L=\partial ^2/\partial a^2+\partial ^2/\partial b^2$
. The free-surface displacement due to monochromatic regular waves of frequency
$\omega$
and amplitude
$A$
is
\begin{equation} \eta =\textrm {Re}\left \{A\mathrm {e}^{\mathrm {i}(ka-\omega t)}\right \}, \end{equation}
where
$\mathrm {i}$
denotes the imaginary unit. Let us assume the following decomposition of the Lagrangian coordinates:
\begin{equation} \left \{x_1,z_1\right \}=\left \{x_1^p+x_1^r, z_1^p+z_1^r\right \}, \end{equation}
where the superscript
$p$
denotes the irrotational (potential) component, and
$r$
denotes the rotational part. Substitution of (3.8) into (3.4)–(3.6) gives (Mei et al. Reference Mei, Stiassnie and Yue2005)
\begin{align} \left (\frac {\partial x_1^p}{\partial t},\frac {\partial z_1^p}{\partial t}\right )=\nabla _L\Phi ,\quad \nabla ^2_L\Phi =0,\quad \Phi _t=-\frac {P_1}{\rho }-gz_1, \quad\end{align}
\begin{align} \frac {\partial x_1^r}{\partial a}+\frac {\partial z_1^r}{\partial b}=0, \qquad\qquad\qquad\qquad\quad\end{align}
\begin{align} \frac {\partial ^2 x_1^r}{\partial t^2}=\nu\, \nabla _L^2\left (\frac {\partial x_1^r}{\partial t}\right )+2\frac {\partial \nu }{\partial a}\frac {\partial ^2 x_1^r}{\partial t\, \partial a}+\frac {\partial \nu }{\partial b}\left (\frac {\partial ^2 z_1^r}{\partial t\, \partial a}+\frac {\partial ^2 x_1^r}{\partial t\, \partial b}\right ), \end{align}
\begin{align} \frac {\partial ^2 z_1^r}{\partial t^2}=\nu\, \nabla _L^2\left (\frac {\partial z_1^r}{\partial t}\right )+2\frac {\partial \nu }{\partial b}\frac {\partial ^2 z_1^r}{\partial t\, \partial b}+\frac {\partial \nu }{\partial a}\left (\frac {\partial ^2 z_1^r}{\partial t\, \partial a}+\frac {\partial ^2 x_1^r}{\partial t\, \partial b}\right ), \end{align}
where
$\Phi$
is the velocity potential for irrotational flow satisfying Laplace’s equation, whereas the rotational components
$(x_1^r,z_1^r)$
are affected by fluid viscosity.
The thickness of the viscous boundary layer is much smaller than the water depth, so the irrotational components
$x_1^p$
and
$z_1^p$
satisfy
$\partial \Phi /\partial b=0$
at
$b=-h$
, and
$\partial \Phi /\partial b=\partial \eta /\partial t$
at
$b=0$
. Hence
\begin{align} x_1\sim x_1^p=\textrm {Re}\!\left \{\!\frac {\mathrm {i}A\cosh [k(h+b)]\mathrm {e}^{\mathrm {i}(ka-\omega t)}}{\sinh (kh)}\!\right \}\!,\quad z_1\sim z_1^p=\textrm {Re}\!\left \{\!\frac {A\sinh [k(h+b)]\mathrm {e}^{\mathrm {i}(ka-\omega t)}}{\sinh (kh)}\!\right \}\!,\notag\\[2pt] \end{align}
which correspond to the inviscid solution of monochromatic waves propagating over horizontal water depth
$h$
(Mei et al. Reference Mei, Stiassnie and Yue2005).
The velocity components in the seabed boundary layer satisfy
$x_1=0$
at
$b=-h$
, and
$x_1\rightarrow x_1^p$
as
$(b+h)/\delta \rightarrow \infty$
, which is now known. In this region the potential part and the kinematic viscosity
$\nu$
do not change significantly, and we obtain the well-known solution (Longuet-Higgins Reference Longuet-Higgins1953)
\begin{align} x_1=\textrm {Re}\!\left \{\!\frac {\mathrm {i}A\left (1-\mathrm {e}^{(\mathrm {i}-1)({b+h})/{\delta }}\right )\mathrm {e}^{\mathrm {i}(ka-\omega t)}}{\sinh (kh)}\!\right \}\!,\!\!\quad z_1=\textrm {Re}\!\left \{\!\frac {\delta k\, A\left (\mathrm {e}^{(\mathrm {i}-1)({b+h})/{\delta }}-1\right )\mathrm {e}^{\mathrm {i}(ka-\omega t)}}{(1-\mathrm {i})\sinh (kh)}\!\right \}\!.\notag\\[2pt] \end{align}
Similarly, in the free-surface boundary layer, the potential part remains almost constant, and the hydrodynamics is dominated by viscous and turbulent effects. By assuming the boundary layer approximation, the
$x$
-momentum equation valid in the free-surface region can be approximated by
\begin{equation} \frac {\partial ^2 x_1^r}{\partial t^2}=\frac {\partial }{\partial b}\left (\nu\, \frac {\partial ^2 x_1^r}{\partial t\,\partial b}\right ). \end{equation}
The corresponding boundary and matching conditions needed for the solution of the above diffusion-like equation can be derived from (2.10). Specifically, by applying zero shear stress at
$b=0$
, and matching with the outer flow (3.13) as
$b/\delta \rightarrow -\infty$
, we obtain
\begin{align} \frac {\partial ^2 x_{1}}{\partial t\,\partial b}=-\frac {\partial ^2 z_{1}}{\partial t\,\partial a},\quad b=0, \end{align}
\begin{align} \frac {\partial x_{1}}{\partial t}=\frac {\partial x_{1}^p}{\partial t},\quad \frac {b}{\delta }\rightarrow -\infty . \end{align}
By using (3.8) and recognising that
$\partial x_{1}^r/\partial b\gg \partial z_{1}^r/\partial a$
from the boundary layer approximation, expressions (3.16)–(3.17) become
\begin{align} \frac {\partial ^2 x_{1u}^r}{\partial t\,\partial b}=\textrm {Re}\left \{-2\omega kA\,\mathrm {e}^{\mathrm {i}(ka-\omega t)}\right \},\quad b=0, \end{align}
\begin{align} \frac {\partial x_{1w}^r}{\partial t}=0,\quad \frac {b}{\delta }\rightarrow -\infty .\qquad\qquad \end{align}
By applying the boundary conditions (3.18)–(3.19), the horizontal rotational velocity component can, in principle, be determined by solving (3.15) once
$\nu$
has been prescribed.
Many researchers have investigated the effects of depth-dependent eddy viscosity with the aim of deriving closed-form expressions for practical applications. For instance, a linear depth-increasing eddy viscosity profile was recently adopted by Constantin et al. (Reference Constantin, Dritschel and Paldor2020) and Lewis & Belcher (Reference Lewis and Belcher2003), whereas Huang & Mei (Reference Huang and Mei2003) considered a quadratic variation of
$\nu$
throughout the entire water depth. Regarding numerical simulations, McWilliams et al. (Reference McWilliams, Sullivan and Moeng1997) applied a large-eddy simulations technique, and found that
$\nu$
had a convex profile, as did Sentchev et al. (Reference Sentchev, Yaremchuk, Bourras, Pairaud and Fraunié2023) from the acoustic Doppler current profiler. By contrast, Shen et al. (Reference Shen, Triantafyllou and Yue2000) studied turbulent diffusion characteristics beneath a free surface, reporting a rapid increase in eddy viscosity with depth, followed by a gradual decrease in the bulk region. For finite water depth, Svendsen (Reference Svendsen1984) assumed an exponential decrease in prescribed eddy viscosity, whereas other authors assumed constant values of eddy viscosity, largely due to uncertainties in available measurements and numerical computations (Dutykh Reference Dutykh2009; Tian et al. Reference Tian, Perlin and Choi2010; Higgins et al. Reference Higgins, Vanneste and van den Bremer2020; Klein et al. Reference Klein, Heifetz and Toledo2022). With respect to Lagrangian coordinates, Weber & Melsom (Reference Weber and Melsom1993) and Christensen & Weber (Reference Christensen and Weber2005) adopted an eddy viscosity profile that increased linearly with
$b$
in the region close to the free surface. A similar approach was previously employed by Miles (Reference Miles1993), who assumed that
$\nu$
is constant along each streamline.
Although the foregoing models are relatively simplistic, they can approximate the averaged motion without the need for complex numerical simulations. Consequently, since we expect strong turbulence near the free surface, and based on these studies, we adopt the following piecewise eddy viscosity profile:
\begin{equation} \nu =\left [\nu _s-\frac {b}{l}(\nu _m-\nu _s)\right ]H[b+l]+\nu _m\,{\rm e}^{\beta (b+l)}\,H[-b-l], \end{equation}
where
$H$
denotes the Heaviside step function,
$\nu _s$
is the eddy viscosity value at the free surface, and
$\nu _m \geqslant \nu _s$
is the maximum eddy viscosity within the free-surface boundary layer. The location where
$\nu = \nu _m$
is represented by
$b = -l$
, while
$\beta$
is a positive constant. In other words, we assume that the eddy viscosity initially increases linearly with depth as in Weber & Melsom (Reference Weber and Melsom1993), and for
$b \lt -l$
, it begins to decrease exponentially at a rate determined by the value of
$\beta$
because of decrease of turbulent effects. We point out that profile (3.20) is similar to the one used in Christensen & Weber (Reference Christensen and Weber2005).
Solution of (3.15) subject to (3.18)–(3.19) and (3.20) is finally given by
\begin{align} x_1^r&=\textrm {Re}\left \{\left \{c_1 I_0\left [\frac {(1-\mathrm {i})\sqrt {l\omega \left [l\nu _s+b\left (\nu _s-\nu _m\right )\right ]}}{\sqrt {2}(\nu _m-\nu _s)}\right ]\right .\right .\nonumber \\ &\left .\left .\quad{}+c_2 K_0\left [\frac {(1-\mathrm {i})\sqrt {l\omega \left [l\nu _s+b\left (\nu _s-\nu _m\right )\right ]}}{\sqrt {2}(\nu _m-\nu _s)}\right ]\right \}\mathrm {e}^{\mathrm {i}(ka-\omega t)}\right \},\quad -l\lt b\lt 0, \end{align}
\begin{align} x_1^r&=\textrm {Re}\left \{c_3\,{\rm e}^{{-a(b+l)}/{2}}\, K_1\left [-\frac {2(-1)^{{3}/{4}}\,{\rm e}^{-a({b+l})/{2}}}{a}\sqrt {\frac {\omega }{\nu _m}}\right ]\mathrm {e}^{\mathrm {i}(ka-\omega t)}\right \},\quad b\lt -l, \end{align}
where
$c_1, c_2, c_3$
are long complex constants not reported here for brevity, and
$I_n, K_n$
are the modified Bessel functions of the first and second kind and of order
$n$
(Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2007). The solution above can be simplified significantly if eddy viscosity is constant, i.e.
$\nu =\nu _s=\nu _m$
and
$\beta \rightarrow 0$
. In this case, we recover the well-known result (Ünlüata & Mei Reference Ünlüata and Mei1970)
\begin{equation} x_{1}^r=\textrm {Re}\big \{Ak(1-\mathrm {i})\delta\,\mathrm {e}^{{(1-\mathrm {i})b}/{\delta }}\,\mathrm {e}^{\mathrm {i}(ka-\omega t)}\big \}, \end{equation}
where
$\delta =\sqrt {2\nu /\omega }$
is the viscous boundary layer thickness. The vertical velocity component
$z_{1}^r$
can be found by applying the continuity equation (3.10). Using boundary layer scales, we obtain
\begin{equation} z_1^r\sim \textit {O}\big (Ak^2\delta ^2\big )\,\text {m}. \end{equation}
In a typical ocean environment,
$A\sim \textit {O}(1)$
m,
$k\sim \textit {O}(10^{-1})$
m
$^{-1}$
and boundary layer thickness is up to
$\delta \sim \textit {O}(10^{-1})$
m. Therefore,
$z_1^r\sim O (10^{-4} )$
m, which is taken to be negligible because of its very minor effect. The total first-order solution in the free-surface boundary layer is finally given by
\begin{equation} x_1=\textrm {Re}\big \{\mathrm {i}A\coth (kh)\,\mathrm {e}^{\mathrm {i}(ka-\omega t)}\big \}+x_{1}^r,\quad z_1=\textrm {Re}\big \{A\,\mathrm {e}^{\mathrm {i}(ka-\omega t)}\big \}. \end{equation}
We now evaluate the wavenumber
$k$
. The effect of a viscous layer on the free surface is to produce spatial damping of water waves, so the wavenumber
$k$
should be complex in principle (Weber Reference Weber1987; Jenkins & Jacobs Reference Jenkins and Jacobs1997). Using continuity of normal stress (2.11) at
$b=0$
at
$O(Ak)$
, we find
\begin{equation} -P_{1}+2\rho \nu\, \frac {\partial ^2 z_{1}}{\partial t\,\partial b}=0. \end{equation}
Given that
$z_1$
remains almost constant in the boundary layer, use of
$P_1=-\rho (\Phi _t+gz_1)$
gives
$\omega ^2=gk\tanh (kh)$
. In other words, we have recovered the standard linear dispersion relation for free-surface water waves propagating over horizontal seabed
$h$
(Ünlüata & Mei Reference Ünlüata and Mei1970). Higher-order wave-attenuation effects can therefore be neglected, as by Ng (Reference Ng2004).
3.3. Second-order solution
The horizontal and vertical momentum equations at second order for generalised eddy viscosity read
\begin{align} \frac {\partial ^2 x_2}{\partial t^2}=-g\frac {\partial z_2}{\partial a}-\frac {1}{\rho }\frac {\partial P_2}{\partial a}+\nu\, \nabla _L^2\left (\frac {\partial x_2}{\partial t}\right )+\frac {\partial \nu }{\partial a}\frac {\partial ^2 x_2}{\partial t\,\partial a}+\frac {\partial \nu }{\partial b}\left (\frac {\partial ^2 x_2}{\partial t\,\partial b}+\frac {\partial ^2 z_2}{\partial t\,\partial a}\right )+\mathcal {G}, \end{align}
\begin{align} \frac {\partial ^2 z_2}{\partial t^2}=-g\frac {\partial z_2}{\partial b}-\frac {1}{\rho }\frac {\partial P_2}{\partial b}+\nu\, \nabla _L^2\left (\frac {\partial z_2}{\partial t}\right )+\frac {\partial \nu }{\partial b}\frac {\partial ^2 z_2}{\partial t\,\partial b}+\frac {\partial \nu }{\partial a}\left (\frac {\partial ^2 z_2}{\partial t\,\partial a}+\frac {\partial ^2 x_2}{\partial t\,\partial b}\right )+\mathcal {H}, \end{align}
where the forcing terms
$\mathcal {G}$
and
$\mathcal {H}$
are given by quadratic products of first-order solutions, i.e.
\begin{align} \mathcal {G}&=-\frac {\partial ^2x_1}{\partial t^2}\frac {\partial x_1}{\partial a}-\frac {\partial ^2z_1}{\partial t^2}\frac {\partial z_1}{\partial a}+\nu \left [\frac {\partial x_1}{\partial a}\frac {\partial }{\partial t}\left (\frac {\partial ^2 x_1}{\partial a^2}+3\frac {\partial ^2 x_1}{\partial b^2}\right )+2\frac {\partial ^3 x_1}{\partial t\,\partial a^2}\frac {\partial z_1}{\partial b} \right. \nonumber \\ &\quad+\frac {\partial ^2 x_1}{\partial t\,\partial a}\left (\frac {\partial ^2 z_1}{\partial a\,\partial b}-\frac {\partial ^2 x_1}{\partial b^2}\right )-2\frac {\partial ^3 x_1}{\partial t\,\partial a\,\partial b}\left (\frac {\partial x_1}{\partial b}+\frac {\partial z_1}{\partial a}\right )\nonumber \\ &\quad\left .+\frac {\partial ^2 x_1}{\partial t\,\partial b}\left (\frac {\partial ^2 x_1}{\partial a\,\partial b}-\frac {\partial ^2 z_1}{\partial a^2}\right )+\frac {\partial z_1}{\partial a}\frac {\partial }{\partial t}\left (\frac {\partial ^2 z_1}{\partial a^2}+\frac {\partial ^2 z_1}{\partial b^2}\right )\right ]\nonumber \\ &\quad+\frac {\partial \nu }{\partial a}\left (2\frac {\partial z_1}{\partial b}\frac {\partial ^2 x_1}{\partial t\,\partial a}-\frac {\partial z_1}{\partial a}\frac {\partial ^2 x_1}{\partial t\,\partial b}-\frac {\partial x_1}{\partial b}\frac {\partial ^2 z_1}{\partial t\,\partial a}-\frac {\partial x_1}{\partial b}\frac {\partial ^2 x_1}{\partial t\,\partial b}+\frac {\partial z_1}{\partial a}\frac {\partial ^2 z_1}{\partial t\,\partial a}\right )\nonumber \\ &\quad +\frac {\partial \nu }{\partial b}\left (-3\frac {\partial z_1}{\partial a}\frac {\partial ^2 x_1}{\partial t\,\partial a}+3\frac {\partial x_1}{\partial a}\frac {\partial ^2 x_1}{\partial t\,\partial b}-\frac {\partial x_1}{\partial b}\frac {\partial ^2 x_1}{\partial t\,\partial a}+\frac {\partial x_1}{\partial a}\frac {\partial ^2 z_1}{\partial t\,\partial a}\right ), \end{align}
and
\begin{align} \mathcal {H}&=-\frac {\partial ^2x_1}{\partial t^2}\frac {\partial x_1}{\partial b}-\frac {\partial ^2z_1}{\partial t^2}\frac {\partial z_1}{\partial b}+\nu \left [\frac {\partial x_1}{\partial b}\frac {\partial }{\partial t}\left (\frac {\partial ^2 x_1}{\partial a^2}+\frac {\partial ^2 x_1}{\partial b^2}\right )+\frac {\partial z_1}{\partial b}\left (\frac {\partial ^3 z_1}{\partial t\,\partial a^2}-\frac {\partial ^2 z_1}{\partial b^2}\right )\right .\nonumber \\ &\quad \left .-\frac {\partial ^2 z_1}{\partial t\,\partial b}\left (\frac {\partial x_1}{\partial b}+\frac {\partial z_1}{\partial a}\right )-\frac {\partial ^2 z_1}{\partial t\,\partial a}\left (\frac {\partial ^2 x_1}{\partial a^2}+\frac {\partial ^2 x_1}{\partial b^2}\right )-\frac {\partial ^2 z_1}{\partial t\,\partial b}\left (\frac {\partial ^2 z_1}{\partial a^2}+\frac {\partial ^2 z_1}{\partial b^2}\right )\right ]\nonumber \\ &\quad+\frac {\partial \nu }{\partial b}\left (2\frac {\partial x_1}{\partial a}\frac {\partial ^2 z_1}{\partial t\,\partial b}-\frac {\partial x_1}{\partial b}\frac {\partial ^2 z_1}{\partial t\,\partial a}-\frac {\partial z_1}{\partial a}\frac {\partial ^2 z_1}{\partial t\,\partial a}-\frac {\partial z_1}{\partial a}\frac {\partial ^2 x_1}{\partial t\,\partial b}+\frac {\partial x_1}{\partial b}\frac {\partial ^2 x_1}{\partial t\,\partial b}\right )\nonumber \\ &\quad+\frac {\partial \nu }{\partial a}\left (-3\frac {\partial x_1}{\partial b}\frac {\partial ^2 z_1}{\partial t\,\partial b}+3\frac {\partial z_1}{\partial b}\frac {\partial ^2 z_1}{\partial t\,\partial a}-\frac {\partial z_1}{\partial a}\frac {\partial ^2 z_1}{\partial t\,\partial b}+\frac {\partial z_1}{\partial b}\frac {\partial ^2 z_1}{\partial t\,\partial b}\right ). \end{align}
Pierson (Reference Pierson1962) derived similar equations for constant
$\nu$
, whereas the above include terms multiplied by
$\partial \nu /\partial a$
and
$\partial \nu /\partial b$
.
The steady component
$U_L=\overline {\partial x_2/\partial t}$
at second order is the mass transport velocity and is determined by time averaging over a wave period. Following a procedure similar to that of Ng (Reference Ng2004), we first integrate the time-averaged vertical momentum (3.28) and use the time-averaged inviscid condition at the free surface
$\overline {P_2}=0$
to obtain (Ünlüata & Mei Reference Ünlüata and Mei1970)
\begin{equation} \overline {z_2}+\frac {\overline {P_2}}{\rho g}-\overline {\eta _2}=-\int _c^0\overline {\mathcal {H}}=\int _c^0\frac {\omega ^2Ak^2\sinh [2k(b+h)]}{2\sinh ^2(kh)}\,\mathrm {d}b, \end{equation}
where the overline denotes an averaged value, and
$\mathcal {H}$
is dominated by the inviscid terms. In the latter equation,
$\eta _2$
denotes the wave set-up and the right-hand side depends only on
$b$
because quadratic products are independent of
$a$
. Assuming
$\nu =\nu (b)$
(as in (3.20)) and zero vertical flux, and substituting (3.31) into (3.28), we obtain the following equation for horizontal mass transport velocity:
\begin{align} \frac {\partial }{\partial b}\left (\frac {\partial U_L}{\partial b}\nu \right )&=-\overline {\mathcal {G}}\nonumber \\ &=-\frac {k\omega }{2}\textrm {Re}\left \{\nu \left [k^2\left (\left |x_1\right |^2+\left |z_1\right |^2\right )-4x_1^*\frac {\partial ^2 x_1}{\partial b^2}-\mathrm {i}k\frac {\partial (x_1^*z_1)}{\partial b}\vphantom{-3\left |\frac {\partial x_1}{\partial b}\right |^2}\right .\right .\nonumber \\ &\left .\left .\quad-3\left |\frac {\partial x_1}{\partial b}\right |^2 -z_1^*\frac {\partial ^2 z_1}{\partial b^2}\right ]-4\frac {\partial \nu }{\partial b}\left [\mathrm {i}x_1^*z_1+x_1\frac {\partial x_1^*}{\partial b}\right ]\right \}, \end{align}
where a
$*$
denotes the complex conjugate. The above is valid across the entire depth and can be solved by deriving the asymptotic behaviour of the horizontal velocity in the boundary layers. By assuming
$O(\delta )\leqslant O(A)$
, in the free-surface boundary layer region, (3.32) can be approximated by
\begin{equation} \frac {\partial }{\partial b}\left (\frac {\partial U_L}{\partial b}\nu \right )=2A\omega k\coth {(kh)}\left [\textrm {Im}\left \{\frac {\partial }{\partial b}\left (\nu \frac {\partial x_1^r}{\partial b}\right )\right \}+2Ak\frac {\partial \nu }{\partial b}\right ], \end{equation}
and integrated by applying a condition of zero stress at the free surface,
\begin{equation} \frac {\partial U_{L}}{\partial b}=0,\quad b=0. \end{equation}
After algebraic manipulation, we get
\begin{equation} \frac {\partial U_{L}}{\partial b}\rightarrow 4A^2k^2\omega \coth {(kh)},\quad \frac {b}{\delta }\rightarrow -\infty , \end{equation}
at the outer edge of the free-surface boundary layer, as also obtained by Ünlüata & Mei (Reference Ünlüata and Mei1970). Equation (3.35) is the first upper boundary condition for the mass transport velocity in the inner core governed by (3.32). A second boundary condition can be inferred by applying the well known behaviour of mass transport velocity at the seabed (Longuet-Higgins Reference Longuet-Higgins1953)
\begin{equation} U_{L}\rightarrow \frac {5A^2k\omega }{4\sinh ^2(kh)},\quad \frac {b+h}{\delta }\rightarrow \infty . \end{equation}
Having determined the boundary conditions for the inner core velocity, the governing equation (3.32) admits the following inviscid solution (Ünlüata & Mei Reference Ünlüata and Mei1970; Ng Reference Ng2004):
\begin{equation} U_{Li}=\frac {A^2k \omega }{4}\left \{\frac {3+2\cosh [2k(b+h)]}{\sinh ^2{(kh)}}+\frac {8k}{\tanh (kh)}\left (h+b\right )\right \}. \end{equation}
Note that the Lagrangian mass transport outside the boundary layer is analogous to the case when the kinematic viscosity is constant. Conversely, the velocity in the free-surface boundary layer region is affected by
$\nu (b)$
, and can be found by solving the approximated governing equation (3.33) jointly with (3.34) and the matching condition for the outer inviscid flow velocity
$U_{Li}$
. This yields the following velocity profile that is valid throughout the water depth affected by the free-surface boundary layer:
\begin{equation} U_{L}=U_{Li}+2Ak\omega \coth (kh)\,\textrm {Im}\left \{x_1^r\right \}. \end{equation}
The expression above is new and includes the effects of spatial variations of
$\nu (b)$
into
$x_1^r$
. By comparing (3.38) against (3.37), we further note that
$U_L$
is always smaller than its inviscid approximation
$U_{Li}$
. This is of crucial significance, especially in view of practical applications devoted to free-surface ocean mass transport and heat transfer.
In the simple case of constant eddy viscosity, we can substitute (3.23) into (3.38) and obtain the explicit solution
\begin{align} U_{L}=\frac {A^2k \omega }{4}\!\left \{\frac {3+2\cosh [2k(b+h)]}{\sinh ^2{(kh)}}+\frac {8k}{\tanh (kh)}\!\left \{h+b-\delta\, \mathrm {e}^{b/\delta }\!\left [\cos \!\left (\!\frac {b}{\delta }\!\right )+\sin \!\left (\!\frac {b}{\delta }\!\right )\!\right ]\right \}\!\right \}\!,\notag\\[2pt] \end{align}
which is also new and shows the exponential decaying component due to the presence of the free-surface boundary layer. We point out that taking the deep-water limit, (3.38) becomes identical to expression (5.2) derived by Weber (Reference Weber1983). Having obtained
$U_L$
, we now turn to the convection–diffusion problem for fluid temperature advected by (3.38).
4. Free-surface temperature
Lagrangian coordinates
$(a,b)$
are now used to solve the Lagrangian temperature equation (2.8). One of the advantages of using (2.8) with respect to its Eulerian counterpart is the absence of convective terms. However, the diffusion term in Lagrangian form (2.9) becomes more complicated than the standard Laplacian operator, and approximations of (2.8) are necessary to find explicit analytical solutions of practical interest. Here, we focus on the free-surface thermal boundary layer because convection and diffusion effects are mainly concentrated in the upper ocean surface layer. By substituting
$x,z$
into the governing equation (2.8), and assuming boundary layer scales, i.e. by non-dimensionalising
$a,b$
as
\begin{equation} b'=b/\delta _T,\quad a'=ak,\quad a'\sim b'\sim O(1), \end{equation}
we find that at leading order,
\begin{equation} \frac {\partial T}{\partial t}=\frac {\partial }{\partial b}\left (\chi \frac {\partial T}{\partial b}\right ), \end{equation}
where
$\delta _T\ll \lambda$
denotes the thermal boundary layer thickness. Significant benefits accrue from using the parabolic equation (4.2), which is considerably simpler than (2.8) and enables us to obtain the analytical solutions reported later in this section.
We now define the boundary conditions for the temperature problem, and consider the effects of prescribed free-surface temperature
$T=T_s(a,t)$
. For simplicity, we assume the Dirichlet boundary condition
\begin{equation} T_s=T_0\left (H[a+U_0t]-H[a+U_0t-L]\right ), \quad b=0, \end{equation}
where
$T_0$
is the temperature of the heat source on the free surface, whereas
$U_0$
is the value of the Lagrangian mass transport (3.38) evaluated on
$b=0$
. From (3.1), we observe that (4.3) represents a heat source of finite Lagrangian width
$L$
, oscillating with the wave field at frequency
$\omega$
. This is further evidenced by performing the transformation from Lagrangian to Eulerian coordinates, as detailed later in this section.
The Lagrangian mass transport velocity is an
$O(\epsilon ^2)$
effect for the velocity field; therefore, it is natural to question whether this effect is important for the temperature problem. To this end, we first define the moving coordinate
$s = a + U_0 t$
. Expression (4.2) then converts into the following governing equation of convection–diffusion type:
\begin{equation} \frac {\partial T}{\partial t} +U_0 \frac {\partial T}{\partial s} = \frac {\partial }{\partial b}\left (\chi \frac {\partial T}{\partial b}\right ). \end{equation}
Meanwhile, the boundary condition (4.3) becomes
\begin{equation} T_s = T_0 \left ( H[s] - H[s - L] \right ), \quad b = 0. \end{equation}
Since
$U_0\sim O(A^2k\omega )$
, and the thickness of the thermal boundary layer is
$\delta _T\sim 3.64\sqrt {\chi L/U_0}$
(Michele et al. Reference Michele, Borthwick and van den Bremer2023), we find that the convective and diffusion terms in the expression above are
$O(\mathcal {T}_i U_0k)$
and
$O(\mathcal {T}_iU_0/3.64^2L)$
, respectively, i.e. diffusion effects are smaller than convection contributions because
$L\sim O(k^{-1})$
. Hence we find that the Lagrangian mass transport velocity results in convective effects that dominate over diffusion. This leads to the following steady-state boundary value problem:
\begin{align} \frac {\partial }{\partial b}\left (\chi \frac {\partial T}{\partial b}\right )-U_0\frac {\partial T}{\partial s}=0, \quad b\in (-\infty ,0], \end{align}
\begin{align} T=T_0\left (H[s]-H[s-L]\right ), \quad b=0, \end{align}
\begin{align} T\rightarrow 0, \quad b\rightarrow -\infty . \qquad\end{align}
Solution of (4.6)–(4.8) can be found by first applying the exponential Fourier transform (Mei Reference Mei1997):
\begin{equation} \tilde {T}(\alpha ,b)=\int _{-\infty }^\infty \mathrm {e}^{-\mathrm {i}\alpha s}\,T\,\mathrm {d}s,\quad T(s,b)=\frac {1}{2\unicode{x03C0} }\int _{-\infty }^{\infty }\mathrm {e}^{\mathrm {i}\alpha s}\,\tilde {T}\,\mathrm {d}\alpha . \end{equation}
Eddy viscosity and thermometric conductivity
$\chi$
are assumed equal owing to the Reynolds analogy; i.e. we assume Prandtl number
$Pr =1$
and that temperature is transported by turbulent eddies (Landau & Lifshitz Reference Landau and Lifshitz1989). Substitution of (4.9) into (4.6)–(4.8) yields the following solution of the transformed problem:
\begin{align} \tilde {T}&=\tilde {T}_s \!\left \{c_4 I_0\left [2\sqrt {\frac {\mathrm {i}lU_0\alpha }{\chi _m-\chi _s}\left (\frac {l\chi _s}{\chi _m-\chi _s}-b\right )}\right ]+c_5 K_0\left [2\sqrt {\frac {\mathrm {i}lU_0\alpha }{\chi _m-\chi _s}\left (\frac {l\chi _s}{\chi _m-\chi _s}-b\right )}\right ]\right \}\!,\nonumber \\ &\quad-l\lt b\lt 0, \end{align}
\begin{align} \tilde {T}&=\tilde {T}_s c_6\,{\rm e}^{-\beta({b+l})/{2}}\sqrt {\alpha }\, K_1\left [\frac {2(-1)^{{1}/{4}}\,{\rm e}^{-\beta({b+l})/{2}}}{\beta }\sqrt {\frac {U_0\alpha }{\chi _m}}\right ],\quad b\lt -l, \end{align}
where
$c_4$
,
$c_5$
and
$c_6$
are long complex constants,
$\chi _s=\nu _s$
,
$\chi _m=\nu _m$
, and
$\tilde {T}_s$
is the Fourier transform of the boundary condition at the free surface, i.e.
\begin{equation} \tilde {T}_s=T_0\int _{0}^L \mathrm {e}^{-\mathrm {i}\alpha s}\,\mathrm {d}s=\frac {\mathrm {i}T_0}{\alpha }\left (\mathrm {e}^{-\mathrm {i}L\alpha }-1\right ). \end{equation}
Substitution of (4.10) and (4.11) into (4.9) allows us to find the solution of the temperature field in integral form, then integration can be performed numerically as undertaken previously for transient dispersive waves (Michele et al. Reference Michele, Renzi, Borthwick, Whittaker and Raby2022).
In evaluating the free-surface heat flux, we derive the normal derivative with respect to a material curve expressed in Lagrangian form by first calculating partial derivatives with respect to the Eulerian coordinates (Monin & Yaglom Reference Monin and Yaglom1971; Piedra-Cueva Reference Piedra-Cueva1995):
\begin{align} \frac {\partial T}{\partial x}=\frac {\partial T}{\partial a}\frac {\partial z}{\partial b}-\frac {\partial T}{\partial b}\frac {\partial z}{\partial a},\quad \frac {\partial T}{\partial z}=\frac {\partial T}{\partial b}\frac {\partial x}{\partial a}-\frac {\partial T}{\partial a}\frac {\partial x}{\partial b}, \end{align}
where we apply conservation of mass (2.5). The normal derivative of
$T$
with respect to the normal vector
$\boldsymbol {n}$
is consequently given by
\begin{equation} \frac {\partial T}{\partial \boldsymbol {n}}=\left (\frac {\partial T}{\partial x},\frac {\partial T}{\partial z}\right ) {\boldsymbol \cdot} \frac {\boldsymbol {n}}{\left |\boldsymbol {n}\right |}=\frac {\dfrac {\partial T}{\partial z}\dfrac {\partial x}{\partial a}-\dfrac {\partial T}{\partial x}\dfrac {\partial z}{\partial a}}{\left |\boldsymbol {n}\right |}=\frac {\partial T}{\partial b}\left |\boldsymbol {n}\right |-\frac {1}{\left |\boldsymbol {n}\right |}\frac {\partial T}{\partial a}\left (\frac {\partial x}{\partial a}\frac {\partial x}{\partial b}+\frac {\partial z}{\partial a}\frac {\partial z}{\partial b}\right ). \end{equation}
By assuming the usual boundary layer scales, the corresponding normal heat flux through the free surface is
\begin{equation} q=-\kappa \frac {\partial T}{\partial \boldsymbol {n}}=-\kappa \frac {\partial T}{\partial b}\sqrt {\left (\frac {\partial x}{\partial a}\right )^2+\left (\frac {\partial z}{\partial a}\right )^2},\quad b=0, \end{equation}
where
$\kappa =\chi \rho c_{p}$
is the thermal conductivity, and
$c_{p}$
is the specific heat at constant pressure. Here, the square root represents the effects of the free-surface displacement on the thermal boundary layer thickness. The total heat flux
$Q$
exchanged across the heat source is given by
\begin{equation} Q=\int _{\left .x\right |_{s=0,b=0}}^{\left .x\right |_{s=L,b=0}}q[x(s,b=0),t]\,\sqrt {1+\left (\frac {\partial \eta [x(s),t]}{\partial x}\right )^2}\,\mathrm {d}x, \end{equation}
where quantities are expressed in Eulerian coordinates, and the square root is related to the definition of arc length. The integral above is evaluated by finding the relation
$s=s(x)$
on the free surface
$b=0$
. In the previous section, we derived
$x=x(a,b)$
,
$z=z(a,b)$
, i.e. the particle trajectory dependent on the initial position; hence using an asymptotic approach as suggested by Pizzo et al. (Reference Pizzo, Lenain, Rømcke, Ellingsen and Smeltzer2023), we get
\begin{equation} ks=kx+\epsilon \frac {\sin (kx-\omega t)}{\tanh (kh)}+\epsilon ^2\frac {\sin [2(kx-\omega t)]}{2\tanh (kh)^2}+O(\epsilon ^3), \end{equation}
and the following expression for the free-surface elevation in Eulerian coordinates:
\begin{equation} k\eta =\epsilon \cos (kx-\omega t)-\epsilon ^2\frac {\sin (kx-\omega t)^2}{\tanh (kh)}+O(\epsilon ^3). \end{equation}
We observe that the definition of
$\eta$
in (3.7) introduces a second harmonic component proportional to
$\epsilon ^2$
in the Eulerian framework. Substituting (4.15), (4.17) and (4.18) into (4.16), and performing a Taylor expansion about
$\epsilon \rightarrow 0$
, we obtain
\begin{equation} Q=-\kappa \int _0^L\left .\frac {\partial T}{\partial b}\right |_{b=0}\left \{\left [1-\epsilon \frac {\cos (ks-\omega t)}{\tanh (kh)}\right ]^2+\epsilon ^2\sin (ks-\omega t)\right \}\,\mathrm {d}s. \end{equation}
To better evaluate the effects of free-surface waves on heat transfer, we calculate the averaged total heat flux per unit source length
$L$
as
\begin{equation} \overline {Q}=\frac {\omega }{2\unicode{x03C0} L}\int _0^{{2\unicode{x03C0} }/{\omega }}Q\,\mathrm {d}t=-\frac {\kappa }{L}\left [1+\frac {\epsilon ^2}{2}\left (1+\frac {1}{\tanh (kh)}\right )\right ]\int _0^L\left .\frac {\partial T}{\partial b}\right |_{b=0}\,\mathrm {d}s, \end{equation}
where we assume
$T$
at steady state and independent of time. The term in brackets represents the positive
$O(\epsilon ^2)$
effect arising from wave steepness, indicating that water waves enhance heat transfer. A similar result was derived previously by Witting (Reference Witting1971) (see their expression (3.4)) for the case of deep water, where
$kh \rightarrow \infty$
. This outcome holds even if
$T$
is independent of Lagrangian coordinate
$s$
, which corresponds to scenarios involving a large heat source with a constant surface temperature, such as daylight heating over extensive areas. This also suggests that capillary waves or breaking waves could further increase the value of
$\overline {Q}$
due to the greater arc length and consequent increase in surface area (Fedorov & Melville Reference Fedorov and Melville1998; Melville & Fedorov Reference Melville and Fedorov2015).
The previous expressions simplify significantly if thermometric conductivity is constant across water depth. In this case, the temperature problem (4.6)–(4.8) can be solved easily by using a Fourier transform, as in Michele et al. (Reference Michele, Borthwick and van den Bremer2023). We obtain the expressions for temperature
\begin{equation} T=T_0\left \{H[s]\,\text {Erfc}\left [b\,\sqrt {\frac {U_0}{4\chi s}}\right ]-H[s-L]\,\text {Erfc}\left [b\,\sqrt {\frac {U_0}{4\chi (s-L)}}\right ]\right \} \end{equation}
and heat flux
\begin{equation} q=-\kappa T_0\sqrt {\frac {U_0}{\unicode{x03C0} \chi }}\left \{\frac {H[s]}{\sqrt {s}}-\frac {H[s-L]}{\sqrt {s-L}}\right \}\sqrt {\left [1+\epsilon \frac {\cos (ka-\omega t)}{\tanh (kh)}\right ]^2+\epsilon ^2 \sin (ka-\omega t)^2}. \end{equation}
We find that heat flux
$q$
is always negative and directed towards the fluid domain for
$s\in [0,L]$
, whereas
$q$
is positive for
$s\gt L$
. Note that
$q\rightarrow \pm \infty$
for
$s\rightarrow 0^-$
and
$s\rightarrow L+0^+$
, respectively, because of the discontinuous gradients in
$T$
at the edges of the heat source. The total averaged heat flux per source length
$L$
becomes
\begin{equation} \overline {Q}=-2T_0\kappa \left [1+\frac {\epsilon ^2}{2}\left (1+\frac {1}{\tanh (kh)}\right )\right ]\sqrt {\frac {U_0}{L\chi \unicode{x03C0} }}. \end{equation}
Since the above depends on the Lagrangian transport velocity (3.38), the effect of the free-surface viscous boundary layer is to reduce the heat flux relative to the inviscid velocity approximation (3.37), thereby enhancing the underwater temperature. Furthermore, we observe that as
$\omega$
increases,
$Q$
also increases. This occurs due to the corresponding increases in
$\epsilon$
and
$U_0$
with wave frequency. We point out that the opposite effect is observed for heat flux at the seabed (Michele et al. Reference Michele, Borthwick and van den Bremer2023), where the corresponding mass transport velocity decreases with
$\omega$
(Longuet-Higgins Reference Longuet-Higgins1953).
5. Results and discussion
In this section, we first represent the Lagrangian mass transport velocity for
$A=0.5$
m,
$h=5$
m, and different values of viscosity. Then we focus on local temperature fields and heat flux exchange at the free surface.
5.1. Lagrangian mass transport velocity
For waves of large amplitude, the laminar boundary layers at the seabed and at the free surface become unstable, and turbulence may be generated (Brocchini & Peregrine Reference Brocchini and Peregrine2001). The effect of ocean turbulence can be included by considering simplified models based on constant eddy viscosity with values much larger than the kinematic viscosity of water (
$10^{-6}$
m
$^2$
s
$^{-1}$
) – a well-established approach in physical oceanography (see e.g. Weber & Christensen Reference Weber and Christensen2019). We point out that
$\nu \sim O(10^{-2}){-}O(10^{-3})$
m
$^2$
s
$^{-1}$
corresponds to eddy viscosity values previously adopted in physical oceanography to model the turbulent free-surface ocean boundary layer (Huang Reference Huang1979; Dutykh Reference Dutykh2009; Higgins et al. Reference Higgins, Vanneste and van den Bremer2020). Similar values have been used to represent breaking-waves effects (Tian et al. Reference Tian, Perlin and Choi2010), leading to significant mass transport (Deike et al. Reference Deike, Pizzo and Melville2017), which is not considered in this work.
Figure 1. (a) Variation in near free-surface Lagrangian mass transport velocity
$U_L$
with vertical coordinate
$b$
for
$A=0.5$
m,
$h=5$
m,
$\omega =1.5$
rad s
$^{-1}$
and (constant) viscosity
$\nu _s=\nu _m=0,10^{-4},10^{-3},10^{-2}$
m
$^2$
s
$^{-1}$
. The black line represents the inviscid approximation
$U_{Li}$
. (b) Eddy viscosity profiles defined by (3.20) and listed in table 1. (c) Plot of
$U_L$
versus
$b$
for
$A=0.5$
m,
$h=5$
m,
$\omega =1.5$
rad s
$^{-1}$
for the eddy viscosity profiles in (b). (d) Behaviour of free-surface mass transport velocity
$U_{0}$
versus
$h$
,
$\omega$
and fixed
$A=0.5$
m,
$\nu _s=\nu _m= 10^{-2}$
m
$^2$
s
$^{-1}$
.
Figure 1(a) shows the behaviour of the Lagrangian velocity
$U_L$
(3.38) with respect to the Lagrangian coordinate
$b$
for
$\omega = 1.5 \, \mathrm {rad} \, \mathrm {s}^{-1}$
and constant eddy viscosity values
$\nu _s = \nu _m = 10^{-4}, 10^{-3}, 10^{-2} \ \mathrm {m}^2 \, \mathrm {s}^{-1}$
. In the same figure, we also depict the inviscid approximation
$ U_{Li}$
(3.37), represented by a black curve. We observe that the inviscid approximation is characterised by the highest velocity at the free surface, whereas for
$b \lt \delta$
, each profile quickly approaches
$U_{Li}$
. We now examine the impact of various eddy viscosity profiles (3.20) on the behaviour of
$ U_L$
for
$\omega = 1.5 \, \mathrm {rad} \, \mathrm {s}^{-1}$
. Table 1 provides a summary of nine profiles along with their respective parameters, where
$\delta _m = \sqrt {2\nu _m/\omega }$
. For clarity, these nine eddy viscosity profiles are also depicted in figure 1(b). Profile 1 corresponds to the case of constant eddy viscosity. Profiles 2–5 are characterised by an exponential decay starting from
$ b = -\delta _m/2$
, whereas the decay of eddy viscosity in profiles 6–9 occurs at a deeper value,
$ b = -3\delta _m$
. Profiles 3, 5, 7 and 9 are characterised by a linear increase in eddy viscosity from the free surface (
$ b = 0$
), whereas the remainder exhibit a constant value before the decay begins, i.e.
$\nu _s = \nu _m$
in
$-l \lt b \lt 0$
. Finally, the different values
$\beta = 4.88, 48.8 \ \mathrm {m}^{-1}$
determine the slow or rapid exponential decay in
$ b \lt -l$
, respectively. Figure 1(c) illustrates the Lagrangian mass transport velocity
$U_L$
for each of the profiles, along with the inviscid approximation
$U_{Li}$
. We immediately note that the constant profile 1 exhibits the smallest value of
$U_0$
, whereas
$U_{Li}$
has the greatest mass transport velocity value at the free surface. Profile 2 features constant eddy viscosity up to
$b = -\delta _m/2$
, beyond which it decays exponentially. Here, the velocity profile resembles that of profile 1, though with a slightly higher value in the range
$-l \lt b \lt 0$
. For
$b \ll -l$
, the profile mirrors the behaviour of
$U_{Li}$
due to the exponential decay of
$\nu$
. Profiles 3, 5, 7 and 9, all characterised by the smallest
$\nu _s$
, yield velocity profiles that closely align with the inviscid approximation. This indicates that when the free-surface eddy viscosity
$\nu _s$
is small, the Lagrangian mass transport velocity can be approximated as
$U_L \sim U_{Li}$
with good accuracy. The main difference occurs in a region close to the free surface, where viscosity effects are dominant. Specifically, profiles 7 and 9 allow greater mass transport velocity because of slower increase in
$\nu$
with depth. Also, we note that the exponential decay does not have a significant effect. Profile 4 represents a profile that has constant eddy viscosity in the range
$-\delta _m/2 \lt b \lt 0$
, and a rapid (exponential) decrease in
$\nu$
for
$b\lt -\delta _m/2$
. In particular, the corresponding mass transport velocity has a bump in the region close to
$b=-l$
, then quickly approaches
$U_{Li}$
due to a rapid exponential decay in
$\nu$
. Curves for profiles 6 and 8 are practically coincident with that for profile 1. This is due to the constant value of
$\nu$
in the range
$-3\delta _m \lt b \lt 0$
, and the exponential decay of eddy viscosity occurring outside the boundary layer thickness
$\delta _m$
. We note that although the curves for mass transport velocity are indistinguishable, the corresponding eddy viscosity and thermometric conductivity profiles result in significant differences in the temperature field and heat transfer, as analysed in the next subsection. Figure 1(d) shows the mass transport velocity at the free surface
$U_0$
versus wave frequency and water depth for constant eddy viscosity
$\nu _s=\nu _m=10^{-2}$
m
$^2$
s
$^{-1}$
. The velocity
$U_{0}$
increases with
$\omega$
as also shown by (3.39). Furthermore,
$U_0$
increases or decreases with water depth
$h$
, depending on whether the value of wave frequency is greater or smaller than
$\omega \sim 1.5$
rad s
$^{-1}$
, underlining the importance of the new finite-depth results that we have derived. We highlight that among the important consequences, the greater
$U_{0}$
, the greater the total heat flux
$Q$
, as also shown in (4.23).
Table 1. Parameters describing the eddy viscosity profiles defined by (3.20) and shown in figure 1(b). The corresponding Lagrangian mass transport velocity behaviour is depicted in figure 1(c).
Figure 2. Steady-state normalised temperature field
$T/T_0$
in the plane (
$a/L$
,
$b/\delta _T$
) for
$A=0.5$
m,
$h=5$
m and
$\nu =\chi$
. (a–i) The results for profiles 1–9, respectively, listed in table 1.
5.2. Temperature field and heat flux
We now investigate the behaviour of the steady-state temperature field (4.10)–(4.11) and heat flux generated by a heat source of finite length
$L=20$
m for fixed water depth
$h=5$
m, wave amplitude
$A=0.5$
m, wave frequency
$\omega =1.5$
rad s
$^{-1}$
, wavenumber
$k=0.264$
m
$^{-1}$
, and the eddy viscosity profiles adopted in the previous subsection. Figure 2 shows the normalised temperature fields
$T/T_0$
in the plane
$(b/\delta _T, x/L)$
, where
$\delta _T\sim 3.64\sqrt {\chi _m L/U_0}$
for each of the nine profiles listed in table 1. We observe that the thermal boundary layer thickness is approximately
$b \sim \delta _T$
at
$a \sim L$
, and is larger than the viscous boundary layer thickness
$\delta$
. This can also be predicted by recalling that the order of magnitude of
$\delta _T$
is
$O(3.64\sqrt {\chi }/(\epsilon \sqrt {\omega }))$
, i.e.
$\delta \sim \epsilon \delta _T$
, which implies that
$\delta$
is approximately an order of magnitude smaller than
$\delta _T$
(Michele et al. Reference Michele, Borthwick and van den Bremer2023). Furthermore, we note that in each case, the growth of this thermal boundary layer is very rapid at
$a \sim 0$
and slows down as
$a$
increases. Figure 2(a) shows that the case of constant thermometric conductivity (profile 1) is characterised by the greatest normalised thermal boundary layer thickness because diffusion properties remain constant with depth. Similar behaviour is displayed in figure 2(f) representing profile 6. The temperature fields depicted in figures 2(b) and 2(c) represent the results for profiles 2 and 3, and are similar to each other. However, the sudden increase in
$\nu$
for profile 3 causes greater rightward bending of
$T$
in figure 2(c) due to a greater
$U_L$
. In figures 2(d) and 2(e) we observe a smaller normalised thickness, primarily due to the rapid exponential decay in
$\chi$
for
$b \lt -l$
. The temperature is almost uniform in
$-l\lt b\lt 0$
, and this will result in a small heat flux due to the gentle temperature gradient as
$b$
approaches 0. Figures 2(f) and 2(g) depict the temperature fields corresponding to profiles 6 and 7. We note that the boundary layer thickness is similar to what is shown in figure 2(a). However, significant differences can be observed in figure 2(g) at small values of
$b$
, where
$\chi$
increases linearly. Figures 2(h) and 2(i) show the results for profiles 8 and 9. There are some similarities with the temperature fields depicted in figures 2(f) and 2(g), except that the thermal boundary layer thickness is smaller. This is because the exponential decay of
$\chi$
is much more rapid, preventing the temperature from spreading into the region
$b \lt -l$
. Let us now analyse the behaviour of the averaged steady-state heat flux ratio
$\overline {Q}/\overline {Q}_c$
, where
$\overline {Q}_c$
is the heat flux (4.23) evaluated in inviscid flow conditions, i.e. using (3.37) instead of (3.38), and constant thermometric conductivity throughout the water depth
$\chi =\chi _m$
. Figure 3(a) illustrates the behaviour of
$\overline {Q}/\overline {Q}_c$
with varying wave frequency for the eddy viscosity profiles listed in table 1. We observe that the curves for profiles 1, 6 and 8 reach the largest values. Specifically, the blue curve for profile 1 is approximately 0.9, representing the effect of the decrease in Lagrangian mass transport velocity within the viscous layer relative to the inviscid approximation
$U_{Li}$
. Notably, the heat flux tends to increase with the square root of
$U_0$
, as also shown in (4.23). This highlights that neglecting the effect of the free-surface viscous boundary layer leads to an overestimation of free-surface heat transfer. Profiles 4 and 5 share a common feature: the rapid exponential decrease of thermometric conductivity in
$b \lt -\delta _m/2$
. This implies that the temperature gradient in
$b\sim 0$
, and consequently
$\overline {Q}$
, is smaller than in the case of constant eddy viscosity. As a result,
$T$
remains almost constant and confined to a region very close to
$b = 0$
(see also figure 2
d,e). This explains why the green and purple curves in figure 3(a) attain values much smaller than unity. Results for profiles 3, 7 and 9 are similar to one another, with curves
$\sim 0.5$
. This is associated with a decrease in
$\chi$
towards the free surface, and a steep temperature gradient, resulting in a significant heat transfer ratio despite
$\nu _s\ll \nu _m$
. The reason why this behaviour is not observed for profile 5, which is also characterised by the linear decrease in
$\chi$
, is that
$\chi$
rapidly reaches
$\chi = \chi _m$
at
$b = -\delta _m/2$
before decreasing rapidly at greater depths (large
$\beta$
). This is also evident in figure 2(e), where it can be seen that
$T$
is almost constant and entirely confined within the free-surface boundary layer. Profiles 2 and 3 give similar results; this means that for small
$l$
and
$\beta$
, the impact of
$\chi _s$
being smaller than
$\chi _m$
is not significant. We conclude by analysing
$\overline {Q}$
for the same eddy viscosity profiles. Figure 3(b) presents the behaviour of the averaged heat flux versus wave frequency, where each curve corresponds to a different profile case. These curves can be broadly classified into three groups based on the free-surface thermometric conductivity and
$l$
. Profiles with higher
$\chi _s$
and
$l$
(profiles 1, 6 and 8) are associated with greater free-surface heat flux. Variations within the same group are attributed to differences in the shapes of the thermometric conductivity and eddy viscosity profiles, as discussed previously. Note that
$\overline {Q}$
reaches values of
$O(-10^4){-}O(-10^5)$
W m
$^{-2}$
, which are significantly higher than the measured latent heat fluxes in the ocean, typically of order
$O(10^2)$
W m
$^{-2}$
(Yu & Weller Reference Yu and Weller2007).
Figure 3. (a) Averaged heat flux ratio
$\overline {Q}/\overline {Q}_c$
and (b) averaged heat flux
$\overline {Q}$
versus wave frequency, for
$A=0.5$
m,
$h=5$
m, and each eddy viscosity profile listed in table 1.
6. Conclusions
This paper has described a mathematical model of free-surface heat transfer forced by long-crested progressive waves in the presence of free-surface and seabed viscous boundary layers. The velocity field and steady-state convection–diffusion equation for water temperature are solved by means of Lagrangian coordinates and Fourier transform. We derived a new form of Lagrangian mass transport velocity that is valid throughout the entire water depth, and obtained analytical solutions that account for spatial variations in eddy viscosity. Our investigation into the effects of various eddy viscosity profiles revealed that neglecting the free-surface viscous boundary layer leads to overestimation of the Lagrangian mass transport velocity. This has a crucial impact, particularly for the analysis of mass transport phenomena coupled with wave-induced hydrodynamics, such as the spreading of microplastics or pollutants. We then analysed the effect of heat sources of finite length localised at the free surface, and demonstrated that the total heat flux is significantly influenced by the boundary layer structure and the behaviour of thermometric conductivity within the free-surface thermal layer. Specifically, we found that when the thermometric conductivity decays rapidly just beneath the free surface, there is a decrease in free-surface heat transfer. Conversely, when the eddy viscosity increases slowly from the free surface, the opposite occurs, resulting in greater heat flux values. In addition, we derived an explicit expression for the temperature field that is valid for constant eddy viscosity, and found that heat flux increases with the square root of the Lagrangian mass transport velocity. This implies that neglecting the effects of the free-surface viscous layer leads to overestimation of heat transfer. We also discovered that the free-surface pattern positively enhances heat flux, yielding second-order effects proportional to the square of wave steepness. This result is independent of the behaviour of eddy viscosity or thermometric conductivity, and was previously identified for deep-water waves by other authors. However, here the result is generalised and valid for finite water depth.
The present mathematical model is limited to two-dimensional viscous and thermal boundary layers driven by simple heat source distributions, without accounting for complex effects such as the coupling between momentum and thermodynamic equations. It should be noted that our model can be extended to random waves by superposition of monochromatic components. However, when the wave field varies spatially, as in random waves, an Eulerian mean velocity forms, known as the return flow underneath wave groups (Longuet-Higgins & Stewart Reference Longuet-Higgins and Stewart1962), and the velocity field becomes more complex. This extension would facilitate a more detailed investigation of mass transfer phenomena in real seas. Expanding the model to encompass more complex domains, fluid properties and velocity fields could also provide valuable insights into additional topics of considerable practical interest in fluid mechanics, such as free-surface boundary layer instability influenced by temperature. Additional effects not considered here, such as those due to capillary waves and breaking waves, are obviously important. We have shown that the increase in surface area yields greater heat flux, and these phenomena are potentially capable of making important contributions. In this work, we adopt the simplistic assumption of prescribing the temperature on the free surface (i.e. a Dirichlet boundary condition). Although this approach is suitable for laboratory-controlled environments, in realistic oceanic environments, the heat flux should be prescribed instead (i.e. a Neumann boundary condition), and additional effects such as evaporation and wind should be considered to capture more accurately the complex heat transfer mechanisms on the free surface. Analysis of these airside processes should be addressed in future work, and we expect the results in the present paper to act as a starting point for such an analysis. In addition, an accurate understanding of subsurface temperature distribution, particularly in the presence of heterogeneous fluid layers, is advantageous in other research areas, such as melting sea ice. Furthermore, the model could serve as a basis for extending satellite measurements of ocean surface temperature to subsurface layers, enabling a more comprehensive understanding of boundary layer temperature profiles without relying solely on in situ sensors.
Acknowledgements.
The authors are grateful to the reviewers for their insightful comments.
Declaration of interests.
The authors report no conflict of interest.
Open access
