1. Introduction
The cross-shore exchange of heat, mass and momentum is a key problem in coastal oceanography and limnology. In general, cross-shore transport in the inner shelf region can be forced by an array of mechanisms including tides, surface waves and internal waves and wind stress, over a wide range of time scales (cf. Lentz & Fewings Reference Lentz and Fewings2012).
Buoyancy-driven flows can also provide an important mechanism for cross-shore exchange over non-uniform bathymetry. Considering an idealized wedge domain, relevant to many near-shore lake and coastal ocean settings, surface buoyancy fluxes across variable depth can set up cross-shore thermal variations with associated baroclinic pressure gradients. The resulting baroclinic exchange can affect flushing times for near-shore waters and deep water renewal (Ivanov et al. Reference Ivanov, Shapiro, Huthnance, Aleynik and Golovin2004; Molina et al. Reference Molina, Pawlak, Wells, Monismith and Merrifield2014), with implications for coastal ecology and near-shore water quality (Hatcher Reference Hatcher1997; Monismith et al. Reference Monismith, Genin, Reidenbach, Yahel and Koseff2006).
On seasonal time scales, winter cooling can drive cascading flows that rapidly flush shelf regions and produce significant mixing (Cooper & Vaux Reference Cooper and Vaux1949; Hill et al. Reference Hill, Souza, Jones, Simpson, Shapiro, McCandliss, Wilson and Leftley1998; Wang & Symonds Reference Wang and Symonds1999; Fer, Lemmin & Thorpe Reference Fer, Lemmin and Thorpe2002; Wang Reference Wang2005). On shorter time scales, diurnal cycles in surface heat fluxes can create a ‘thermal siphon’, as described by Monismith, Imberger & Morison (Reference Monismith, Imberger and Morison1990), where shallower near-shore waters heat and cool more rapidly than deeper regions, driving diurnally reversing flows. Thermally driven diurnal exchanges have been observed in lakes (Monismith et al. Reference Monismith, Imberger and Morison1990) and tropical reefs (Neiman et al. Reference Neiman, Richter, Jonkers and Badran2004; Monismith et al. Reference Monismith, Genin, Reidenbach, Yahel and Koseff2006; Molina et al. Reference Molina, Pawlak, Wells, Monismith and Merrifield2014; Herdman, Hench & Monismith Reference Herdman, Hench and Monismith2015), suggesting that they may be a persistent mechanism for cross-shore exchange in these environments. For tropical reefs, this mechanism for daily cross-shore heat flux may be important in moderating thermal stress that can lead to coral bleaching (Hoegh-Guldberg Reference Hoegh-Guldberg1999; Monismith Reference Monismith2007; Safaie et al. Reference Safaie, Silbiger, McClanahan, Pawlak, Barshis, Hench, Rogers, Williams and Davis2018).
Cross-shore exchange is an essential physical process controlling larval connectivity on reefs and driving the shoreward transport of plankton that make up a large portion of the nutrient supply required to support reef metabolism (Hatcher Reference Hatcher1997; Genin et al. Reference Genin, Monismith, Reidenbach, Yahel and Koseff2009; Monismith et al. Reference Monismith2010). Diurnal thermally driven cross-shore exchange may interact with similar time scale biological mechanisms, such as the diel vertical swimming patterns of some larvae, to produce observed cross-shore distributions of larvae and patterns of connectivity among populations (Garland, Zimmer & Lentz Reference Garland, Zimmer and Lentz2002; Tapia et al. Reference Tapia, DiBacco, Jarrett and Pineda2010). Cross-shore residual flow patterns may also influence the variability of pH and dissolved oxygen concentrations experienced by near-shore organisms, which may differ significantly from the offshore environment due to spatial gradients in rates of primary production and respiration (Hofmann et al. Reference Hofmann2011).
Coastal ocean observations of thermally driven exchange suggest that the nature of the flow response to cyclical buoyancy forcing varies significantly in character depending on the dynamical balances that apply to a given site. For example, Molina et al. (Reference Molina, Pawlak, Wells, Monismith and Merrifield2014) described observations on a fore reef slope on Oahu, Hawai'i where strong along-shore flows provided high effective vertical mixing. For this case, high turbulent transport distributes surface buoyancy fluxes vertically and an unsteady thermal response is observed, where the temperature is in quadrature with the surface buoyancy forcing. The resulting cross-shore baroclinic pressure gradient is in phase with the thermal response and drives a momentum response that is similarly in phase (i.e. in quadrature with the surface buoyancy forcing), reflecting a diffusive balance. In contrast, Monismith et al. (Reference Monismith, Genin, Reidenbach, Yahel and Koseff2006) described thermally driven cross-shore flow, observed at 8–16 m depth on a steep reef slope in the Gulf of Aqaba, where advective acceleration was dominant, leading to thermal and momentum responses that are in phase with the surface forcing.
A theoretical framework for buoyancy driven flow on a slope has been established for the steady-state flow response (Sturman, Oldham & Ivey Reference Sturman, Oldham and Ivey1999), the transient response to steady thermal forcing (Lei & Patterson Reference Lei and Patterson2002, Reference Lei and Patterson2005) and the periodic response to harmonic forcing (Farrow & Patterson Reference Farrow and Patterson1993; Farrow Reference Farrow2004). Of particular relevance for our work is the seminal investigation of Farrow & Patterson (Reference Farrow and Patterson1993), who considered the motion induced in a fluid wedge by diurnal heating and cooling, with the thermal forcing introduced as a vertically uniform heating/cooling term in the energy equation. A solution was sought for small values of the bottom slope. In this limit, advection terms were found to be negligible at leading order, thereby enabling a closed-form solution in which the temperature is independent of the depth and the velocity is given by a linear balance between the local acceleration, the viscous force, and the buoyancy-induced horizontal pressure gradient. Higher-order terms were described in a follow-up study by Farrow (Reference Farrow2004), who compared the results of the asymptotic analysis for the infinite wedge with those of numerical integrations in a triangular fluid domain limited by a vertical no-slip boundary. As expected, good agreement was found away from the vertical boundary for moderately large values of the relevant Grashof number, when the solution remains weakly nonlinear, with larger differences emerging when nonlinear effects become prevalent at large Grashof numbers. A triangular fluid domain was also considered by Mao, Lei & Patterson (Reference Mao, Lei and Patterson2009) in a scaling analysis of natural convection induced by absorption of radiation. Also of interest is the investigation of Farrow (Reference Farrow2016), who extended the early work by including interactions of thermal forcing with the simultaneous presence of flow forcing associated with a sea-breeze/gully wind system.
The heating/cooling model used by Farrow & Patterson (Reference Farrow and Patterson1993) neglects the attenuation of solar radiation, instead imposing a depth uniform heat flux that, in turn, results in a vertically uniform temperature structure. A more elaborate model was employed by Lei & Patterson (Reference Lei and Patterson2006) in a numerical solution of the near-shore flow response to periodic thermal forcing in a trapezoidal flow domain, including a wedge-like section of linearly increasing depth followed by a constant-depth section ending at a vertical adiabatic wall. The model uses Beer's law for describing the decrease of the diurnal radiation intensity with depth along with a prescribed heat flux at the water surface for modelling the nocturnal heat loss. The same trapezoidal flow domain was employed in a subsequent numerical investigation (Bednarz, Lei & Patterson Reference Bednarz, Lei and Patterson2009), with the periodic heating/cooling introduced by harmonically varying the surface temperature. Both numerical studies (Lei & Patterson Reference Lei and Patterson2006; Bednarz et al. Reference Bednarz, Lei and Patterson2009) revealed a time lag of the flow response to the switch between day-time heating and night-time surface cooling as well as the existence of offshore thermal instabilities, more prominent in shallow reservoirs.
One important outcome of the previous theoretical works on thermally driven exchange is an understanding that the resulting dynamic regimes can vary with distance from shore and variable depth. For example, Farrow & Patterson (Reference Farrow and Patterson1993) examined the linearized case and found a phase shift between the shallow diffusive region and an unsteady response in deeper offshore waters. The cross-shore phase shift introduced transient convergence/divergence regions with associated down/upwelling regions. Additionally, thermally driven cross-shore exchange can be modulated by other persistent forcing mechanisms such as diurnal winds (Farrow Reference Farrow2013, Reference Farrow2016) or tidal variations in along-shore flows (Molina et al. Reference Molina, Pawlak, Wells, Monismith and Merrifield2014; Ulloa et al. Reference Ulloa, Davis, Monismith and Pawlak2018). Monismith et al. (Reference Monismith, Genin, Reidenbach, Yahel and Koseff2006) examined the dynamics of the flow in the absence of rotation, using scaling of the momentum and thermal energy balance equations to examine different possible regimes of flow behaviour.
Harmonic forcing can give rise to residual components where a steady flow appears due to nonlinear contributions. Residual flows in the coastal zone are most commonly associated with wind-driven flow (Lentz & Fewings Reference Lentz and Fewings2012) or with transport due to surface waves via Stokes drift (Xu & Bowen Reference Xu and Bowen1994; Lentz et al. Reference Lentz, Fewings, Howd, Fredericks and Hathaway2008) and in the bottom boundary layer (Longuet-Higgins Reference Longuet-Higgins1953; Davies & Villaret Reference Davies and Villaret1999). Exchanges at longer time scales can also be driven by low frequency seasonal buoyancy fluxes. The potential for residual flow components associated with thermally driven exchange flows has not been examined previously.
For the flows driven by harmonic thermal forcing analysed by Farrow & Patterson (Reference Farrow and Patterson1993) and Farrow (Reference Farrow2004), at leading order in the limit of small slopes, the velocity and temperature exhibit a harmonic temporal variation with zero time average. The first-order correction to the temperature, computed by Farrow (Reference Farrow2004), displays a steady component driven by the nonlinear heat advection associated with the leading-order solution. Steady-streaming motion (Riley Reference Riley2001), partly driven by the steady buoyancy force associated with the temperature correction and partly driven by nonlinear advective terms, can be expected to appear in the first-order correction to the velocity field. The steady velocities for this residual motion, which were not examined by Farrow (Reference Farrow2004), scale with the square of the slope, and therefore would provide a negligibly small contribution to the cross-shore transport rate. Clearly, other types of harmonic thermal forcing might give rise to a more pronounced residual motion, as could be tested by time-averaging results of numerical computations (e.g. Lei & Patterson Reference Lei and Patterson2006; Bednarz et al. Reference Bednarz, Lei and Patterson2009), but such an analysis has not been attempted.
Here, we explore a theoretical formulation for buoyancy-driven flow in a fluid wedge to show how cyclical thermal forcing may drive not only a harmonic flow response but also a steady residual motion, with both components appearing in the leading-order solution for small bottom slopes. Thermal forcing is introduced through a harmonically varying surface flux, an appropriate representation of heating/cooling in turbid waters, where the radiation penetration distance is much smaller than the local depth. Our analytical approach builds upon previous work by considering vertically varying heat input over sloping bathymetry – a situation common in natural near-shore environments. The solutions presented here, although limited in parameter space, show that a periodic thermal forcing can drive a net advective heat flux and persistent environmental gradients, an important result in connection with natural flows. The occurrence of a steady residual flow is reasoned to be a more general feature for any harmonic wedge flow with a vertically non-uniform thermal response.
2. Formulation
The two-dimensional problem considered here is schematically represented in figure 1. Cartesian coordinates are used in the description, with 
$x$ measuring the horizontal distance to the shore and 
$z$ pointing in the vertical direction opposite to the gravitational acceleration 
$\bar {g}=-g \bar {e}_z$, with 
$\bar {e}_z$ denoting the vertical unity vector. The water layer extends for 
$x>0$ between the surface 
$z=0$ and the bottom 
$z=z_s(x)<0$. The wedge-like near-shore region described here is assumed to be very slender, in that
\begin{equation} -z_s/x \sim \beta \ll 1, \end{equation}
with the average slope 
$\beta$ used as an asymptotically small parameter in the following description.
Figure 1. Schematic view of the model flow configuration used here to investigate buoyancy-induced cross-shore flow.
To model diurnal heating by solar radiation and corresponding nocturnal cooling, a periodic heat flux 
$q=-q_o \cos (\omega t)$ is applied at the water surface, heat being transferred across the water layer by turbulent transport associated with a superposed along-shore motion, with corresponding turbulent thermal diffusivity 
$\kappa _t$ and associated eddy viscosity 
$\nu _t$, both assumed to be constant in the following analysis. As follows from a balance between local accumulation and heat conduction, away from the shore this heat flux results in the appearance of a near-surface Stokes thermal layer of characteristic thickness 
$\delta =(\kappa _t/\omega )^{1/2}$ where we find temperature variations 
$T-T_o$ from the bulk temperature 
$T_o$ of order 
${\rm \Delta} T_c = [q_o/(\rho _o c)]/(\kappa _t \omega )^{1/2}$, with 
$\rho _o$ and 
$c$ representing the bulk water density and the water specific heat, respectively. The associated variations of the density due to heating are given by 
$(\rho _o-\rho )/\rho _o = \alpha (T-T_o)$, with 
$\alpha \ll 1/T_o$ denoting the thermal-expansion coefficient, resulting in changes of the pressure from the unperturbed hydrostatic distribution 
$p'=p-(p_a-\rho _o g z)$ of order 
$p'\sim g \alpha q/(c \omega )$ where 
$p_a$ is the ambient reference pressure. In the near-shore region, where 
$-z_s \sim \delta$, the interactions of the heated layer with the bottom surface produce a horizontal pressure gradient. Considering an unsteady thermal balance with a corresponding balance between local acceleration and the pressure gradient yields a scale for the resulting cross-shore velocity 
$u$ of order 
$u_c=\beta g \alpha [q/(\rho _o c)]/(\kappa _t \omega ^3)^{1/2}$. Since the fluid domain is slender, the accompanying vertical velocity component 
$w$ is much smaller, of order 
$\beta u_c \ll u_c$, as follows from continuity.
The description of the periodic heating and associated induced motion in the near-shore region requires integration of the continuity, momentum and energy equations written in the Boussinesq approximation. The scales identified above are used to define dimensionless variables
\begin{equation} \left.\begin{gathered} x^*=\frac{x}{\delta/\beta},\quad z^*=\frac{z}{\delta},\quad t^*=\omega t,\quad u^*=\frac{u}{u_c},\quad w^*=\frac{w}{\beta u_c},\\ T^*=\frac{T-T_o}{{\rm \Delta} T_c},\quad p^*=\frac{p'}{g \alpha q/(c \omega)}. \end{gathered}\right\} \end{equation}In what follows, the asterisks will be removed to simplify the notation. The problem reduces to that of integrating
\begin{gather} \frac{{\partial} u}{{\partial} x}+\frac{{\partial} w}{{\partial} z} =0 \end{gather}
\begin{gather}\frac{{\partial} u}{{\partial} t}+\beta^2 {{Ra}} \left(\frac{{\partial} (u^2)}{{\partial} x}+ \frac{{\partial} (u w)}{{\partial} z} \right)={-}\frac{{\partial} p}{{\partial} x} +{Pr} \left(\beta^2 \frac{{\partial}^2 u}{{\partial} x^2}+\frac{{\partial}^2 u}{{\partial} z^2}\right) \end{gather}
\begin{gather}\frac{{\partial} w}{{\partial} t}+\beta^2 {{Ra}} \left(\frac{{\partial} (u w)}{{\partial} x}+ \frac{{\partial} (w^2)}{{\partial} z} \right)= \frac{1}{\beta^2}\left(-\frac{{\partial} p}{{\partial} z}+T\right) +{Pr}\left(\beta^2 \frac{{\partial}^2 w}{{\partial} x^2}+\frac{{\partial}^2 w}{{\partial} z^2}\right) \end{gather}
\begin{gather}\frac{{\partial} T}{{\partial} t}+\beta^2 {{Ra}} \left(\frac{{\partial} (u T)}{{\partial} x}+ \frac{{\partial} (w T)}{{\partial} z} \right)= \beta^2 \frac{{\partial}^2 T}{{\partial} x^2}+\frac{{\partial}^2 T}{{\partial} z^2}, \end{gather}
where 
${Pr} = \nu_t/\kappa$ is the turbulent Prandtl number and
\begin{equation} {{Ra}}=\frac{g \alpha [q/(\rho_o c)]}{\omega^2 \kappa_t}, \end{equation}is the relevant Rayleigh number. The necessary boundary conditions include
\begin{equation} \frac{{\partial} u}{{\partial} z}=w=\frac{{\partial} T}{{\partial} z}+\cos t=0 \quad \textrm{at}\ z=0, \end{equation}corresponding to stress-free flow with prescribed heat flux at the free water surface, which is modelled with the rigid-lid approximation, and
\begin{equation} u=w=\frac{{\partial} T}{{\partial} z}-\beta^2 \frac{\textrm{d}z_s}{\textrm{d}x} \frac{{\partial} T}{{\partial} x}=0\quad \textrm{at} \ z=z_s(x), \end{equation}
corresponding to an adiabatic bottom with no-slip flow. As 
$x \rightarrow \infty$ the fluid is at rest with a boundary temperature distribution
\begin{equation} T=\textrm{e}^{z/\sqrt{2}} \cos (t+z/\sqrt{2}+3{\rm \pi}/4) \quad \textrm{as}\ x \rightarrow \infty, \end{equation}
which can be obtained from integration of the associated one-dimensional heat problem arising far from the shore. Besides the specific topography, defined by the negative function 
$z_s(x)$ with 
$|z_s| \sim 1$ for 
$x \sim 1$ and 
$z_s \rightarrow -\beta x$ as 
$x \rightarrow \infty$, and the turbulent Prandtl number 
$ {Pr}=\nu _t/\kappa _t$, the solution depends on two parameters, namely, the Rayleigh number defined in (2.7) and the square of the average bottom slope 
$\beta^2$.
The oscillatory nature of the heat flux applied at the surface implies the existence of unstable thermal stratification over some portion of the cycle. Rayleigh–Bénard instabilities can be expected to develop in flow configurations for which the period of the heat-flux cycle is sufficiently large compared with the characteristic time for instability growth. The tendency toward instability development is measured by the Rayleigh number, the relevant bifurcation parameter governing convective instabilities. Periodic solutions to (2.3)–(2.10) can be expected to become unstable when the Rayleigh number increases above a critical value, leading to modifications in the flow that enhance heat transport rates (see, e.g. Yang et al. (Reference Yang, Chong, Wang, Verzicco, Shishkina and Lohse2020) for a recent analysis of effects of time periodic modulation on Rayleigh–Bénard convection). This aspect of the problem is to be discussed further below in connection with the results of the direct numerical simulations.
An asymptotic solution is obtained below for small values of 
$\beta$ in the distinguished limit 
$\text {Pr} \sim 1$ and 
${Ra} \sim 1$. At leading order in the limit 
$\beta \ll 1$ the problem becomes linear, as is clear from inspection of (2.3)–(2.9). An interesting feature of the associated periodic solution is that, besides harmonic temperature and velocity fields resulting from the harmonic forcing through the surface heat flux, the leading-order solution exhibits a steady component, including a residual temperature increase and associated steady-streaming velocity. The existence of this non-zero steady solution, additional to the expected harmonic response, can be anticipated by taking the time average
\begin{equation} \langle \cdot \rangle =\frac{1}{2 {\rm \pi}} \int_0^{2 {\rm \pi}} \cdot \textrm{d}t \end{equation}of (2.6) to yield
\begin{equation} \beta^2 {{Ra}} \left(\frac{{\partial} \langle u T \rangle}{{\partial} x}+ \frac{{\partial} \langle w T \rangle}{{\partial} z} \right)= \beta^2 \frac{{\partial}^2 \langle T \rangle}{{\partial} x^2}+\frac{{\partial}^2 \langle T \rangle}{{\partial} z^2}. \end{equation}
Integrating in 
$z$ across the fluid layer with account taken of the boundary conditions given in (2.8) and (2.9) provides
\begin{equation} \frac{{\partial}}{{\partial} x} \left({{Ra}} \int_{z_s}^0 \langle u T \rangle \textrm{d} z- \int_{z_s}^0 \frac{{\partial} \langle T \rangle}{{\partial} x} \textrm{d} z \right)=0, \end{equation}which further reduces to
\begin{equation} \int_{z_s}^0 \left({{Ra}} \langle u T \rangle- \frac{{\partial} \langle T \rangle}{{\partial} x} \right)\textrm{d} z=0, \end{equation}
when use is made of the solution as 
$x \rightarrow \infty$. The above exact result, independent of 
$\beta$, indicates that, since the product of two harmonic functions in general does not average to zero, the average temperature gradient must also be non-zero, which in turn induces an average motion with 
$\langle u \rangle \ne 0$ and 
$\langle w \rangle \ne 0$, as seen below.
3. Analytical solution for small slopes
3.1. Leading-order problems
An asymptotic solution to (2.3)–(2.10) in the limit 
$\beta \ll 1$ can be constructed by introducing expansions for the different variables in even powers of 
$\beta$ (e.g. 
$T=T_0+\beta ^2 T_1+\beta ^4 T_2+ \cdots$). The leading-order terms, denoted by a zero subscript, satisfy linear problems, obtained by collecting the dominant terms in (2.3)–(2.9). The solution simplifies by expressing the velocity in terms of the streamfunction 
$\psi _0$, defined such that
\begin{equation} u_0=\frac{{\partial} \psi_0}{{\partial} z} \quad \textrm{and} \quad w_0={-}\frac{{\partial} \psi_0}{{\partial} x}. \end{equation}
Furthermore, guided by the previous discussion, we anticipate that each dependent variable can be expressed at leading order as the sum of a harmonic component and a steady component, the latter linearly proportional to 
${Ra}$, as can be expected from (2.14). Thus, we choose to write the streamfunction and the temperature in the form
\begin{equation} \psi_0= \varPsi_h +{{Ra}} \, \mathcal{F} \quad \textrm{and} \quad T_0=T_h +{{Ra}} \, \mathcal{T}, \end{equation}
where the terms 
$\varPsi _h = \mbox {Re}(\mathrm {e}^{\mathrm {i} t} \varPsi )$ and 
$T_h = \mbox {Re}(\mathrm {e}^{\mathrm {i} t} \varTheta )$, carrying the temporal variation of the solution, are written in terms of the complex functions 
$\varPsi (x,z)$ and 
$\varTheta (x,z)$, which satisfy
\begin{gather} {Pr} \frac{{\partial}^4 \varPsi}{{\partial} z^4}-\mathrm{i} \frac{{\partial}^2 \varPsi}{{\partial} z^2}=\frac{{\partial} \varTheta}{{\partial} x} \end{gather}
\begin{gather}\frac{{\partial}^2 \varTheta}{{\partial} z^2}-\mathrm{i} \varTheta =0 , \end{gather}with boundary conditions
\begin{equation} \varPsi=\frac{{\partial}^2 \varPsi}{{\partial} z^2}=\frac{{\partial} \varTheta}{{\partial} z}+1=0 \quad \textrm{at}\ z=0 \quad \textrm{and} \quad \varPsi=\frac{{\partial} \varPsi}{{\partial} z}=\frac{{\partial} \varTheta}{{\partial} z}=0\quad \textrm{at}\ z=z_s(x). \end{equation}
The second terms on the right-hand sides of (3.2a,b) correspond to the time averages 
$\langle \psi _0 \rangle ={Ra} \, \mathcal {F}$ and 
$\langle T_0 \rangle ={Ra} \, \mathcal {T}$, where the functions 
$\mathcal {F}$ and 
$\mathcal {T}$, satisfy
\begin{gather} {Pr} \frac{{\partial}^4 \mathcal{F}}{{\partial} z^4} = \frac{{\partial} \mathcal{T}}{{\partial} x} \end{gather}
\begin{gather}\frac{{\partial}^2 \mathcal{T}}{{\partial} z^2} = 0 , \end{gather}with boundary conditions
\begin{equation} \mathcal{F}=\frac{{\partial}^2 \mathcal{F}}{{\partial} z^2}= \frac{{\partial} \mathcal{T}}{{\partial} z}=0 \quad \textrm{at}\ z=0 \quad \textrm{and}\quad \mathcal{F}=\frac{{\partial} \mathcal{F}}{{\partial} z}=\frac{{\partial} \mathcal{T}}{{\partial} z}=0\quad \textrm{at}\ z=z_s(x). \end{equation}3.2. Harmonic response
The harmonic component can be expressed in closed form. The development begins by integrating (3.4) with the boundary conditions given in (3.5a,b) to give
\begin{equation} \varTheta=\left(\frac{1-\mathrm{i}}{\sqrt{2}} \right) \frac{\cosh(Z-Z_s)}{\sinh (Z_s)}, \end{equation}which has been expressed in terms of the auxiliary variables
\begin{equation} Z=\frac{1+\mathrm{i}}{\sqrt{2}} z \quad \textrm{and} \quad Z_s=\frac{1+\mathrm{i}}{\sqrt{2}} z_s. \end{equation}Introducing this result into (3.3) and integrating yields
\begin{equation} \frac{\varPsi}{\sinh^{{-}2}(Z_s)\frac{\textrm{d} z_s}{\textrm{d}x}}=A_1 \cosh(Z/\sqrt{{Pr}})+A_2 \sinh(Z/\sqrt{{Pr}})+A_3+A_4 Z +\varPsi_p, \end{equation}
involving the particular solution 
$\varPsi _p(Z)$ and the four functions 
$A_1(Z_s)$, 
$A_2(Z_s)$, 
$A_3(Z_s)$ and 
$A_4(Z_s)$, to be obtained by application of the four homogeneous boundary conditions given for 
$\varPsi$ in (3.5a,b). The solution takes different forms depending on the value of 
$\text {Pr}$. Thus, we find
\begin{equation} \varPsi_p=\frac{\cosh(Z)}{{Pr}-1} \quad \textrm{if}\ {Pr} \ne 1 \quad \textrm{and} \quad \varPsi_p=\frac{Z}{2} \sinh(Z) \quad \textrm{if}\ {Pr} = 1, \end{equation}with corresponding integration constants given by
\begin{gather} A_1=\frac{{Pr}}{1-{Pr}}, \end{gather}
\begin{gather}A_2={-}\frac{\cosh(Z_s)-1-Z_s \sinh(Z_s) -{Pr}\left[\cosh\left(\dfrac{Z_s}{\sqrt{{Pr}}}\right)-1- \left(\dfrac{Z_s}{\sqrt{{Pr}}}\right) \sinh\left(\dfrac{Z_s}{\sqrt{{Pr}}}\right)\right]}{(1-{Pr}) \left[\left(\dfrac{Z_s}{\sqrt{{Pr}}}\right) \cosh\left(\dfrac{Z_s}{\sqrt{{Pr}}}\right)-\sinh \left(\dfrac{Z_s}{\sqrt{{Pr}}}\right)\right]}, \end{gather}
\begin{gather}A_3=1, \end{gather}
\begin{gather}A_4={-}\frac{\cosh\left(\dfrac{Z_s}{\sqrt{{Pr}}}\right) [\cosh(Z_s)-1+{Pr}]/\sqrt{{Pr}}-\sinh(Z_s) \sinh\left(\dfrac{Z_s}{\sqrt{{Pr}}}\right) -\sqrt{{Pr}}} {(1-{Pr})\left[\left(\dfrac{Z_s}{\sqrt{{Pr}}}\right) \cosh\left(\dfrac{Z_s}{\sqrt{{Pr}}}\right)-\sinh\left(\dfrac{Z_s}{\sqrt{{Pr}}}\right)\right]}, \end{gather}
if 
$\text {Pr} \ne 1$ and by
\begin{gather} A_1={-}1, \end{gather}
\begin{gather}A_2={-}\frac{\cosh\left(Z_s\right)-1-Z_s \sinh\left(Z_s\right)+(Z_s/2) \cosh(Z_s)}{Z_s \cosh\left(Z_s\right)-\sinh\left(Z_s\right)}, \end{gather}
\begin{gather}A_3=1, \end{gather}
\begin{gather}A_4={-}\frac{\sinh(Z_s)}{2}-\frac{ \cosh\left(Z_s\right)-1-\left(Z_s/2\right) \sinh\left(Z_s\right) \cosh(Z_s)}{Z_s \cosh\left(Z_s\right)-\sinh\left(Z_s\right)}, \end{gather}
if 
$\text {Pr} = 1$.
3.3. Steady components
In addition to the harmonic velocity and temperature fields described above, the solution at leading order has steady components, including a time-averaged temperature distribution 
$\langle T_0 \rangle (x)={Ra} \mathcal {T}$ that is independent of the depth, as can be seen by integrating (3.7) subject to 
$ {\partial } \mathcal {T}/ {\partial } z=0$ at 
$z=0$ and 
$z=z_s$. The function 
$\mathcal {T}(x)$ can be obtained by extending the analysis given here to a higher order or, alternatively, from the solvability condition
\begin{equation} \frac{\textrm{d} \mathcal{T}}{\textrm{d}x} ={-}\frac{1}{z_s} \int_{z_s}^0 \left\langle \mbox{Re}\left(\mathrm{e}^{\mathrm{i} t} \frac{{\partial} \varPsi}{{\partial} z} \right) \mbox{Re}\left(\mathrm{e}^{\mathrm{i} t} \varTheta \right) \right\rangle \textrm{d}z, \end{equation}
stemming from (2.14). Integrating (3.21) with the boundary condition 
$\mathcal {T}(\infty )=0$, consistent with the far-field distribution (2.10), provides
\begin{equation} \mathcal{T}=\int_x^\infty \frac{1}{z_s} \left[\int_{z_s}^0 \left\langle \mbox{Re}\left(\mathrm{e}^{\mathrm{i} t} \frac{{\partial} \varPsi}{{\partial} z} \right) \mbox{Re}\left(\mathrm{e}^{\mathrm{i} t} \varTheta \right) \right\rangle \textrm{d}z \right] \textrm{d} \tilde{x}, \end{equation}
where 
$\tilde {x}$ is a dummy integration variable. The accompanying streamfunction for the steady motion,
\begin{equation} \mathcal{F}=\frac{z_s^4}{24 {Pr}} \frac{\textrm{d} \mathcal{T}}{\textrm{d}x} \left(\frac{z}{z_s}\right) \left[\left(\frac{z}{z_s}\right)^3 - \frac{3}{2} \left(\frac{z}{z_s}\right)^2 + \frac{1}{2} \right], \end{equation}is obtained in terms of the cross-shore temperature gradient (3.21) by straightforward integration of (3.6) subject to the boundary conditions stated in (3.8a,b).
4. Discussion of the analytical results
The expressions developed above can be used to investigate the resulting flow. Results are given below for the case of a plane bottom surface (i.e. 
$z_s=-x$). Besides analysing the temperature and velocity fields, different aspects of cross-shore exchange will be quantified. A convenient measure of the instantaneous net transport rate associated with the thermally driven flow at a given offshore location is given by the depth-integrated speed
\begin{equation} \varGamma_{ex}=\tfrac{1}{2} \int^0_{z_s} \left| u \right| \, \textrm{d}z, \end{equation}
where the 
$1/2$ factor allows this quantity to be interpreted as an exchange flux, consistent with field studies (e.g. Molina et al. Reference Molina, Pawlak, Wells, Monismith and Merrifield2014). In addition, transport of heat is conveniently measured with use of the cross-shore advective flux
\begin{equation} G_{ex}= \int^0_{z_s} u T \, \textrm{d}z. \end{equation}
The above two quantities and their corresponding time-averaged values 
$\langle \varGamma _{ex} \rangle (x)$ and 
$\langle G_{ex}\rangle (x)$ can be evaluated for the leading-order solution with use of
\begin{equation} u_0=\frac{{\partial} \varPsi_h}{{\partial} z} +{{Ra}} \, \frac{{\partial} \mathcal{F}}{{\partial} z} \end{equation}
and 
$T_0= T_h +{Ra} \, \mathcal {T}$. It is worth mentioning that, since 
$\mathcal {T}$ is independent of 
$z$, its contribution to the advective flux 
$G_{ex}$ vanishes, with the result that the corresponding time-averaged value reduces to
\begin{equation} \langle G_{ex}\rangle(x)= \int^0_{z_s} \left\langle \frac{{\partial} \varPsi_h}{{\partial} z} \, T_h \right\rangle \, \textrm{d}z, \end{equation}independent of the Rayleigh number.
4.1. Harmonic temperature and velocity variations
Streamlines and temperature distributions for the harmonic solution with 
$\text {Pr}=1$ are shown in figure 2 for varying phases of the surface heat-flux cycle. The heat flux is defined so that midnight coincides with the beginning of the cycle, with 
$t={\rm \pi}$ correspondingly representing conditions at noon, when the heating rate is maximum. The solution includes a near-shore region 
$x \sim 1$ where the temperature exhibits cross-shore variation, with the vertical temperature distribution eventually approaching the Stokes profile (2.10) for 
$x \gg 1$. The associated motion is characterized by a series of counterrotating vortices in the 
$x-z$ plane. As can be seen, a pair of vortices is created in the near-shore region during each heat-flux cycle, with the generation from midnight to noon of a counterclockwise rotating vortex followed by the generation of a clockwise rotating vortex from noon to midnight. These vortices weaken as they propagate offshore.
Figure 2. Cross-sections for the harmonic solution with 
$\text {Pr} = 1$ at selected phases of the heating/cooling cycle. Colour indicates instantaneous temperature distribution 
$T_h$ with streamlines 
$\varPsi _h =$ constant overlaid. Streamlines are at intervals of 0.003 with grey contours indicating positive (counterclockwise) rotation and white contours indicating negative (clockwise) rotation.
The resulting phase relationships are more clearly evident in the surface velocity and temperature contours shown in figure 3. For 
$x\ll 1$, both temperature and velocity are in quadrature with the surface heat flux so that the cooling response, with onshore flow at the surface, extends from midnight to noon. For 
$x\gg 1$, the lag in the surface thermal response is diminished, reflecting a balance between unsteady cooling/heating and downward thermal diffusion. Temperature variations approach a constant amplitude offshore, consistent with the boundary condition in (2.10), while the vortical structures weaken exponentially although their propagating response persists as 
$x\rightarrow \infty$.
Figure 3. Temporal variation of the surface distributions of temperature (a) and velocity (b) corresponding to the harmonic solution with 
$\text {Pr} = 1$.
To understand the propagating cell response, we examine the harmonic temperature component, 
$T_h = \mbox {Re}(\mathrm {e}^{\mathrm {i} t} \varTheta )$ using (3.9). For a plane bottom surface (
$z_s = -x$), far from shore (
$x \gg 1$), this reduces to
\begin{align} T_{h} &\simeq \textrm{e}^{z/\sqrt{2}} \cos (t+z/\sqrt{2}+3{\rm \pi}/4) \nonumber\\ &\quad +\exp({-(2x+z)/\sqrt{2}}) \cos [t-\left(2x+z\right)/\sqrt{2}+3{\rm \pi}/4] \end{align}
when exponentially small terms are neglected. The first term on the right-hand side is the solution to the one-dimensional heat problem, given in (2.10). The second term is a perturbation to this solution that accounts for the presence of the lower boundary. This latter component is characterized by propagating temperature bands that follow lines with 
$2x+z = \textrm {const}$. These perturbations to the one-dimensional solution are amplified near the bottom boundary and propagate upward via vertical diffusion. The horizontal gradient that results from the sloping bands thus propagates offshore, driving the translating momentum response apparent in the cells in figures 2 and 3(b).
The persistent cell structure in the harmonic velocity solution is in contrast with that observed when the surface heat flux is assumed to be distributed uniformly throughout the depth, as in the analyses of Farrow & Patterson (Reference Farrow and Patterson1993) and Farrow (Reference Farrow2013). For that case, no propagation is observed for 
$x \rightarrow \infty$. There, the homogeneous vertical distribution of the cyclic surface heat flux ensures that a harmonic cross-shore pressure gradient persists, driving a flow that is in quadrature with the forcing, reflecting an unsteady dynamic balance.
4.2. Time-averaged temperature distribution and associated residual streaming
As indicated in (3.2a,b), the periodic solution includes a steady component proportional to the Rayleigh number. The spatial distributions of the time-averaged temperature and velocity are described by the functions 
$\mathcal {T}(x)$ and 
$\mathcal {F}(x,z)$, which can be evaluated from (3.22) and (3.23), respectively. Illustrative results are given in figure 4, including temperature profiles 
$\mathcal {T}(x)$ for 
$\text {Pr} = 0.1$, 1 and 7 (figure 4a) and streamlines 
$\mathcal {F}=$constant for 
$\text {Pr}=1$ (figure 4b). The dashed curves in figure 4(a) correspond to results of direct numerical integrations, to be discussed later. While results for 
$\text {Pr} = 0.1$ are not environmentally relevant, these are included here to illustrate the trends associated with varying 
$\text {Pr}$.
Figure 4. (a) The solid curves represent the function 
$\mathcal {T}$ evaluated from (3.22) for selected values of 
$\text {Pr}$, while the dashed curves represent distributions of depth-averaged mean temperature 
$\langle T \rangle /{Ra}$ evaluated from numerical integrations of (2.3)–(2.10) with 
$\text {Pr}=1$ and 
${Ra}=5$ for two small values of 
$\beta$. (b) Colour indicates temperature 
$\mathcal {T}$ with corresponding contours of 
$\mathcal {F}$ for 
$\text {Pr} = 1$. Contours are at intervals of 0.001 with grey contours indicating positive (counterclockwise) rotation and white contours indicating negative (clockwise) rotation.
The plots in figure 4(a) indicate that heat tends to accumulate near shore, with the mean temperature reaching its maximum value at 
$x=0$. The magnitude of the normalized temperature 
$\mathcal {T}$, larger for smaller values of 
$\text {Pr}$, remains relatively small. The temperature variation with 
$x$ is non-monotonic. The initial decay, eventually leading to negative values of 
$\mathcal {T}$, is followed by a gentle increase. As can be inferred from (3.23), this non-monotonic variation of 
$\mathcal {T}$ leads to the formation of a sequence of cells with alternating circulatory flow of decreasing magnitude, with the boundaries of each cell defined by the locations where 
$\textrm {d} \mathcal {T}/\textrm {d}x$ vanishes. The steady vortex emerging near shore corresponds to clockwise rotation (i.e. negative values of the normalized streamfunction 
$\mathcal {F}$). The next vortex, with counterclockwise rotation, has a similar magnitude for 
$\mathcal {F}$, (
$\mathcal {F}_{peak}=-1.1\times 10^{-2}$ for the first vortex and 
$\mathcal {F}_{peak}=9.5 \times 10^{-3}$ for the second vortex) although with notably weaker velocities due to the greater depth. The third vortex, whose streamlines do not appear in the figure, is much weaker, and plays a negligible role in the resulting near-shore dynamics.
4.3. Cross-shore exchange for the harmonic solution
The contribution of the harmonic motion to the net cross-shore exchange rate 
$\varGamma _{ex}$ can be investigated by evaluating (4.1) using 
$u = {\partial } \varPsi _h/ {\partial } z$. Corresponding time-averaged values 
$\langle \varGamma _{ex} \rangle$ are represented in figure 5 for the three Prandtl numbers considered earlier in figure 4. Near shore, the exchange increases with increasing total depth, peaking at a distance 
$x \sim 1$. For 
$x\gg 1$, 
$\langle \varGamma _{ex} \rangle$ decreases as the bottom depth exceeds the extent of the surface thermal layer (
$z_s>1$) and the associated cross-shore baroclinic gradient diminishes.
Figure 5. Time-averaged cross-shore exchange flux 
$\langle \varGamma _{ex} \rangle$ as obtained from (4.1) with the harmonic expression 
$u= {\partial } \varPsi _h / {\partial } z$ for 
$\text {Pr}=(0.1,1,7)$ (dotted curves) and with the combined solution (4.3) for 
$\text {Pr}=1$ and 
${Ra} =(5,20)$ (solid curves); the dashed curves represent results of direct numerical simulations with 
$\text {Pr}=1$, 
${Ra} =5$ and 
$\beta =(0.1,0.2)$.
Effects of variations in the Prandtl number on 
$\langle \varGamma _{ex} \rangle$ are most easily interpreted in terms of changes in viscosity, since length scales for the problem are defined in terms of a thermal diffusion length. A reduction in 
$\text {Pr}$ leads to an increase in stream function magnitudes as viscous effects are diminished, as evident in the first two terms in (3.11), with a corresponding increase in velocity magnitude and in 
$\langle \varGamma _{ex} \rangle$. Length scales for the harmonic flow response are primarily determined by the structure of the harmonic thermal response, which is independent of 
$\text {Pr}$. Cross-shore scales for the velocity field decrease slightly for smaller values of 
$\text {Pr}$, however, moving the peak in 
$\langle \varGamma _{ex} \rangle$ shoreward.
We can similarly use the asymptotic results to examine the cross-shore advective flux of heat 
$G_{ex}$. Its time-averaged value 
$\langle G_{ex} \rangle$, independent of the Rayleigh number, can be evaluated for different values of 
$\text {Pr}$ using (4.4), with results given in figure 6. The observed trends, including the effects of variations in the Prandtl number, are consistent with those discussed earlier in connection with figure 5. The time-averaged advective heat flux associated with the harmonic solution is negative close to shore, switching to positive farther offshore (i.e. at 
$x \simeq 5$ for 
$\text {Pr} = 1$). The resulting divergence leads to a heat deficit at 
$x\simeq 5$ and an accumulation of heat toward shore. The shoreward advective heat transport is offset by an offshore-directed horizontal diffusive transport as reflected in the balance expressed in (2.14). The magnitude of the resulting time-averaged cross-shore thermal gradient is thus determined by the relative strength of thermal diffusion, which accounts for the dependence on 
${Ra}$. The pressure difference associated with this thermal gradient then results in a steady component of the streamfunction that is also proportional to 
${Ra}$. The resulting steady flow, shown in figure 4(b) for 
$\text {Pr} = 1$, is characterized by counterrotating cells with boundaries at 
$x$ locations at which the cross-shore thermal gradient vanishes. The cross-shore scales for these cells are modulated by 
$\text {Pr}$, shifting shoreward for decreasing 
$\text {Pr}$.
Figure 6. Time-averaged advective heat flux 
$\langle G_{ex} \rangle$ as obtained from (4.4) for 
$\text {Pr}=(0.1,1,7)$ (solid curves) and as obtained from (4.2) using velocity and temperature fields from direct numerical simulations for 
$\text {Pr}=1$, 
${Ra} =5$ and 
$\beta =(0.1,0.2)$ (dashed curves).
The temporal variation in the harmonic flow contribution to 
$G_{ex}$ is shown in figure 7. Here, 
$G_{ex}$ is evaluated with 
$u= {\partial } \varPsi _h/ {\partial } z$ and 
$T = T_h$, for 
$\text {Pr} = 1$. The observed variations in depth-integrated advective heat flux arise as a result of the phase relationship between the velocity and temperature. Very near to shore (
$x\ll 1$), the temperature response is largely uniform with depth so that the harmonic velocity field cannot drive a depth-averaged flux. For 
$x\gtrsim 1$, the thermal response diminishes and eventually shifts in phase as function of depth, while the velocity exchange structure spans the full depth, as shown in figure 2. The resulting differences in vertical scales lead to a phase relationship that can drive a net flux over the full heat cycle while varying with cross-shore distance. In figure 2, the cell closest to shore at 
$t = 0$ has a negative rotation and spans depths that exceed the thermal layer thickness (
$z>1$). As the cell propagates offshore, it drives cool water offshore at the surface and warm water onshore at depth over the late cooling phase, through 
$t\simeq 3{\rm \pi} /4$ extending offshore to 
$x\simeq 4$. The flow at the surface and at depth therefore both contribute to the negative (shoreward) heat flux evident in the early part of the cycle in figure 7. A new cell, with positive rotation, forms close to shore by 
$t=3{\rm \pi} /4$ (see figure 2) advecting warm water shoreward at the surface and cool water offshore at depth over the late heating phase through 
$t\sim 3{\rm \pi} /2$ and 
$x\sim 4$, resulting in a shoreward heat flux that appears as a second peak at 
$t\sim 5{\rm \pi} /4$, as shown in figure 7. Together, these features result in the peak in cycle-averaged shoreward flux evident in figure 6 at 
$x\simeq 3$. Similar velocity–temperature phasing leads to a weaker seaward flux that peaks at 
$x\simeq 6$ (for 
$\text {Pr} = 1$). As for the cross-shore exchange, the cross-shore advective heat flux associated with the harmonic flow is also function of 
$\text {Pr}$, increasing and shifting shorewards at low 
$\text {Pr}$ as the harmonic velocity magnitude increases.
Figure 7. Depth-integrated advective heat flux 
$G_{ex}$ evaluated from (4.2) using the harmonic components 
$u= {\partial } \varPsi / {\partial } z$ and 
$T= T_h$ for 
$\text {Pr} = 1$.
4.4. Rayleigh number dependence
It is of interest to use the results of the analysis to address effects of increasing Rayleigh numbers, whose influence is carried by the time-independent terms in the asymptotic expressions (3.2a,b). Since the magnitude of the functions 
$\mathcal {F}$ and 
$\mathcal {T}$ is fairly small for all Prandtl numbers of interest, as revealed by the results in figure 4, significant quantitative effects can be expected to emerge only for moderately large values of 
${Ra}$. In this respect, it is worth mentioning that, although the analysis assumes values of 
${Ra}$ of order unity in establishing the asymptotic ordering, the resulting leading-order solution (3.2a,b) remains valid as long as 
${Ra} \ll \beta ^{-2}$, so that effects of advective transport can be neglected in writing the reduced equations (3.3)–(3.8a,b). For slopes with 
$\beta \sim 0.1$ (comparable to steeper coastal settings), the asymptotic ordering leading to (3.2a,b) can thus be expected to hold up to 
${Ra} \sim 10^2$ (additional limitations arising from the development of hydrodynamic instabilities are to be addressed separately in the next section).
As previously mentioned, the analysis indicates that the time-averaged advective flux of heat 
$\langle G_{ex} \rangle$ is independent of the Rayleigh number, so that the predictions given in figure 6 should remain valid irrespective of the value of 
${Ra}$, a result to be tested later by means of direct numerical simulations. By way of contrast, the residual motion yields a non-zero contribution to the time-averaged net cross-shore exchange rate 
$\langle \varGamma _{ex} \rangle$, larger for larger values of 
${Ra}$. This is illustrated in figure 5, where the solid curves represent results corresponding to 
${Ra}=5$ and 
${Ra}=20$ for 
$\text {Pr}=1$. These results are to be compared with the purely harmonic case, depicted with a dotted curve. As expected, all three curves coincide at points where the time-averaged reduced velocity 
$ {\partial } \mathcal {F}/ {\partial } z$ vanishes, corresponding to the boundaries separating the different circulation cells in figure 4(b). It is of interest that, since the definition of 
$\varGamma _{ex}$ involves the instantaneous absolute value of the velocity, the departure of the time-averaged cross-shore exchange rate 
$\langle \varGamma _{ex} \rangle$ from its purely harmonic value is not linearly proportional to 
${Ra}$, as might have been expected from (4.3).
The asymptotic expressions (3.2a,b) are used to generate the cross-sectional plots of temperature and streamlines shown in figure 8 for 
$\text {Pr}=1$ and 
${Ra}=5$, to be compared with those shown in figure 4 for the harmonic solution with the same value of 
$\text {Pr}=1$. As in the harmonic case, the plots reveal the alternating shedding of travelling counterrotating vortices from the near-shore region. The propagating cells interact with the stationary residual vortices (depicted in figure 4b) introducing asymmetry in the near-shore flow response (
$x<5$) with strengthening of the negative vortex for 
$0<t<3{\rm \pi} /4$ and a corresponding weakening of the positive vortex for 
${\rm \pi} <t<7{\rm \pi} /4$. For 
$x>5$, the pattern is reversed, evidenced in the breakdown of the negative cell, apparent offshore beginning around 
$t = 5{\rm \pi} /4$.
Figure 8. Cross-sections for the asymptotic solution (3.2a,b) with 
$\text {Pr} = 1$ and 
${Ra} = 5$ at selected phases of the heat-flux cycle. Colour contours are used for the instantaneous temperature distribution with streamlines 
$\psi _0=$ constant overlaid. Streamlines are at intervals of 0.003 with grey contours indicating positive (counterclockwise) rotation and white contours indicating negative (clockwise) rotation.
The effects of increasing 
${Ra}$ on the velocity field and the advective heat flux are illustrated in figure 9 which can be compared with the associated plots for the harmonic solution in figures 3(b) and 7. The development of the diurnal asymmetry is apparent in the surface-velocity plot, which exhibits strengthened offshore velocity close to shore, associated with the warming phase response, and a corresponding increase in onshore flow farther from shore. As 
${Ra}$ increases, propagating cells offshore are replaced by quasi-steady cross-shore cells, reflecting the steady flow pattern in figure 4(b). Notably, a persistent surface divergence is established near 
$x\simeq 5$.
Figure 9. Combined harmonic and steady solution: surface velocity and depth-integrated advective heat flux for 
$\text {Pr} = 1$ and (a,b) 
${Ra} = 1$, (c,d) 
${Ra} = 5$, (e,f) 
${Ra} = 20$.
For the heat flux, 
$G_{ex}$, although the time-averaged advective flux is unaffected by changes in 
${Ra}$ (figure 6), the product of the steady flow with the harmonic temperature introduces a temporal variation in 
$G_{ex}$ that modulates the diurnal patterns, as evident in figure 9. The superposition of the negative near-shore cell associated with the steady solution, with nightly cooling leads to an increase in onshore heat flux for 
$0<t<{\rm \pi}$, which is offset by an offshore flux during the heating response. The amplitude of the time-variable advective fluxes increase concurrently with increasing 
${Ra}$, with the resulting patterns eventually reflecting the influence of the steady flow cells.
5. Numerical simulations for finite values of the bottom slope
5.1. Computational model
To test the accuracy of the asymptotic results, full numerical solutions of the conservation equations (2.3)–(2.6) were computed for finite values of the slope 
$\beta$. The numerical method employs a second-order Lagrange–Galerkin temporal discretization scheme in combination with P2/P1 Taylor–Hood finite elements for velocities and pressure, which are solved in a coupled fashion (for detailed information about the numerical method see, for example, Carpio, Prieto & Vera Reference Carpio, Prieto and Vera2016). The method was implemented in the open source solver FreeFEM (Hecht Reference Hecht2012), using an interface with the multifrontal massively parallel sparse direct solver (MUMPS) to solve the discretized systems.
The computational domain, shown in figure 10, has a finite extent 
$x_{max}$ in the offshore direction. A clipped corner of size 
$x_{min}$ is introduced near the origin to avoid singularities. The values 
$x_{min} = 0.1$ and 
$x_{max} = 24$ were found to be sufficiently small/large to give well-converged results. The computational mesh was constructed so that the elements decrease in size towards the top surface and towards the near-shore corner, as shown in figure 10. Convergence tests were carried out to ensure the results are converged with respect to mesh refinement and temporal resolution. In particular, the instantaneous temperature and velocity at a central point in the domain, 
$(x=10,\ z=-5)$, was monitored over the course of ten cycles. Maximal changes in this temperature and velocity evolution were tracked between subsequent mesh refinements, each refinement consisting of a doubling of the number of triangular elements in the mesh. Mesh convergence was established when subsequent refinements produced maximal changes in the aforementioned velocity and temperature evolutions of less than 1 %. In the final mesh, which was employed for the computations, the smallest element measures 0.027 and the largest element has size 0.61. With respect to the temporal discretization, an analogous convergence test, in which the time step was decreased in 25 % decrements, showed that a fixed time step of 
$2{\rm \pi} /192$ was sufficiently small to ensure stability and convergence (i.e. a 25 % smaller time step produced changes in the results smaller than 1 %). In the numerical model, the bottom domain boundary 
$\varSigma _b$ and the small vertical near-shore boundary 
$\varSigma _n$ are adiabatic, no-slip surfaces, as dictated by (2.9), whereas the top boundary 
$\varSigma _t$ is a slip boundary subject to the periodically fluctuating heat flux expressed in (2.8). To facilitate convergence, advective boundary conditions 
$\textrm {D} u/\textrm {D} t = \textrm {D}w/\textrm {D}t = \textrm {D}T/\textrm {D}t = 0$ with 
$\textrm {D}()/\textrm {D}t= {\partial }()/ {\partial } t+u {\partial }()/ {\partial } x$ are imposed on the offshore boundary 
$\varSigma _f$.
Figure 10. The computational domain and mesh employed in the numerical simulations.
The numerical computations were started at time 
$t = {\rm \pi}/2$ – such that the heat flux at the surface is zero – using as initial conditions the approximate solution (3.2a,b) for 
$\beta \ll 1$ that contains both the harmonic response at 
$t={\rm \pi} /2$ and the steady residual motion. The computations were carried out until a 
$2{\rm \pi}$ periodic solution was reached, the criterion for convergence being that the root-mean-squared difference between average velocity and temperature fields (integrated over the entire fluid domain) be smaller than 
$10^{-2}$ over subsequent cycles. It was found that the duration of the transient stage leading to the establishment of a periodic solution is very sensitive to the initial conditions. For instance, initializing the computations with an incorrect energetic balance sets up persistent convective currents across the domain that take a prohibitively long time to decay.
5.2. Validation of the analytical results
Numerical results corresponding to 
${Ra}=5$ and 
$\text {Pr}=1$ are shown in figures 4–6. Time averages of the converged periodic solution confirmed the existence of near-shore heat accumulation with an accompanying residual motion. Vertical variations of the time-averaged temperature 
$\langle T \rangle$ at any cross-shore location 
$x$ were verified to be smaller than 
$10^{-3}$, in agreement with the theoretical findings. The cross-shore variation of the depth-averaged mean temperature 
$-\int _{z_s}^0 \langle T \rangle \textrm {d}z/z_s$ obtained numerically for 
$\beta =0.1$ and 
$\beta =0.2$ is compared in figure 4 with the leading-order predictions of the asymptotic analysis for 
$\beta \ll 1$. The resulting agreement is excellent, with observed discrepancies being consistent with the expected errors of order 
$\beta ^{2}$ associated with the asymptotic development.
The velocity field obtained numerically was employed to compute the net transport rate 
$\varGamma _{ex}$ and the cross-shore advective heat flux 
$G_{ex}$, defined in (4.1) and (4.2), respectively. Time-averaged values 
$\langle \varGamma _{ex} \rangle (x)$ and 
$\langle G_{ex} \rangle (x)$ computed for 
$\beta =0.1$ and 
$\beta =0.2$ are compared in figures 5 and 6 with the corresponding asymptotic predictions. The degree of accord is excellent in both cases. Computations for both larger and smaller values of 
${Ra}$ were seen to give distributions of 
$\langle G_{ex} \rangle (x)$ that are virtually indistinguishable from those shown in figure 6 for 
${Ra}=5$, as expected from the analytical findings. Clearly, the consistency displayed by the comparisons in figures 4–6 provide increased confidence in the validity of the analytical results.
5.3. Development of thermal instabilities
As mentioned earlier, the application of a harmonic heat flux at the surface necessarily leads to the appearance of unstable thermal stratification over some portion of the cycle. Instabilities may develop if the fluid remains in an unstable configuration for a sufficiently long time, i.e. if the frequency of the heating cycle 
$\omega$ decreases below a threshold value. Correspondingly, the periodic solution described above can be expected to undergo a bifurcation to a different solution when the Rayleigh number 
${Ra}$ exceeds a critical value 
${Ra}_c$. Since the unperturbed solution is periodic, the linear stability analysis needed to determine 
${Ra}_c$ involves Floquet theory (Davis Reference Davis1976). Results are available for the related case of a horizontal fluid layer bounded by parallel surfaces subject to zero-mean boundary temperature modulation (Gershuni & Zhukhovitskii Reference Gershuni and Zhukhovitskii1963; Or & Kelly Reference Or and Kelly1999). As expected from classical Rayleigh–Bénard analyses (Drazin & Reid Reference Drazin and Reid2004), the instability tendency is seen to augment for increasing layer depths. The stability of a semi-infinite fluid domain subject to a surface heat flux with zero mean, relevant for the conditions found offshore in our flow configuration, has been addressed recently (T. Christopher, personal communication). The associated linear stability analysis reveals that the flow undergoes a period-doubling bifurcation for 
${Ra}_c=12.44$.
To further explore this aspect of the problem, additional numerical computations were performed with 
$\beta =0.1$ and 
$\text {Pr}=1$ by sequentially increasing the value of 
${Ra}$ in unity increments. No instability development was observed up to 
${Ra}=13$, for which the resulting velocity and temperature fields were verified to exhibit excellent agreement with the analytical predictions across the domain. As shown in the temporal evolution of surface velocity in figure 11(a), convergence to a periodic solution with vanishing offshore motion was no longer achieved for 
${Ra}=14$. Instead, thermal-instability cells appeared offshore after a few tens of cycles and were seen to spread rapidly, eventually occupying the entire water surface beyond 
$x \simeq 6$. It is of interest that the offshore convective cells are 
$4{\rm \pi}$-periodic, consistent with the expected double-period bifurcation, while the near-shore solution remains stable owing to the stabilizing effect of the bottom wall, its 
$2{\rm \pi}$-periodic dynamics being still described by the analytical solution given above. Similar transient dynamics can be observed at higher values of 
${Ra}$, for which the instability grows faster. As seen in figure 11(b), compared with the near-critical solution for 
${Ra}=14$, the cells of the bifurcated solution for 
${Ra}=20$ become larger and develop at smaller distances from shore, so that the cross-shore extent of the region of 
$2{\rm \pi}$-periodic stable flow shrinks for increasing 
${Ra}$. These findings are consistent with those of previous numerical simulations (Mao, Lei & Patterson Reference Mao, Lei and Patterson2010).
Figure 11. Contours of surface velocity obtained in numerical computations with 
$\text {Pr} = 1$ and 
$\beta =0.1$ for 
${Ra} = 14$ (a) and 
${Ra} = 20$ (b).
The cross-sectional structure of the double-period thermal instability away from the bifurcation (i.e. for 
${Ra}=20$) is depicted in the snapshots of figure 12, corresponding to 16 evenly spaced phases spanning two periods of the heating/cooling cycle. The size of the non-slender convective cells is seen to be much larger than the characteristic thickness of the Stokes thermal layer 
$\delta$ used in scaling the coordinates 
$x$ and 
$z$. The induced circulatory motion penetrates to depths exceeding 10
$\delta$, creating a non-uniform temperature field, markedly different from that of the stable case 
${Ra}=5$, depicted in figure 8.
Figure 12. Cross-sections for the numerical solution with 
$\text {Pr} = 1$, 
${Ra} = 20$, and 
$\beta =0.1$ at selected instants of time covering two complete heating/cooling cycles. Colour contours are used for the instantaneous temperature distribution with arrows indicating velocities overlaid.
6. Discussion
The theoretical analysis presented here describes the cyclical buoyancy-driven cross-shore exchange that arises within an idealized coastal wedge-shaped domain due to a harmonically varying surface heat flux. An important consequence of applying the heat flux at the surface is that the resulting harmonic temperature distribution is a function of both cross-shore distance and depth, which, coupled with the harmonic velocity field, results in a net advective cross-shore heat flux associated with the leading-order solution.
The surface heat-flux boundary condition used in (2.8) is thus critical to the steady temperature distribution and to the associated steady flow component. If we consider, instead, a heat flux applied uniformly in depth, through the introduction of a source term in (2.6), for example, then, at leading order, we would obtain a depth uniform temperature. The continuity equation (2.3), requires that the depth integral of the cross-shore velocity be zero, so that the horizontal advective heat flux would similarly be zero for all values of 
$x$. Any boundary condition for heat input that leads to a vertically inhomogeneous temperature distribution can introduce a net horizontal advective heat flux, however. For realistic parameters characteristic of environmental flows, this situation might arise from depth-variable radiative input, or from limited vertical mixing. Phase shifts between the vertical structure of temperature fields and cross-shore momentum consequently can lead to net heat fluxes and time-averaged temperature structure in the cross-shore.
For the heating phase, the heat-flux boundary condition used here is consistent with turbid conditions, for which the short-wave radiation input is confined to a near-surface layer. For more general water conditions, an exponentially decaying heat input (i.e. Beer's law, cf. Lei & Patterson Reference Lei and Patterson2006) would be more representative. This would alter the details of the theoretical solution, although we can expect that the resulting vertical heat distribution would similarly result in a non-uniform cross-shore gradient, as suggested by (2.14). Steady streaming components can thus be expected for numerical solutions like that presented by Lei & Patterson (Reference Lei and Patterson2006). For very clear waters, bottom heating can further complicate the vertical distribution of heat, leading to thermal inversions and daytime convective mixing (Mao et al. Reference Mao, Lei and Patterson2009; Wells, Fram & Pawlak Reference Wells, Fram and Pawlak2012); see also Hattori, Patterson & Lei (Reference Hattori, Patterson and Lei2015) for an analysis of the interplay between direct radiative heat absorption in the main water body, which leads to stable thermal stratification, and local heating from the bottom surface, which results in a near-boundary layer with unstable thermal stratification. Radiative input introduces a new length scale associated with the optical path length (or optical depth). In systems where this length scale is comparable to or larger than the heat-conduction thickness 
$\delta$, the cross-shore extent of the thermally driven flow would be larger than that described above. Future work should also consider a depth-uniform heat flux as a more appropriate representation of the cooling phase in systems where vertical convection distributes heat rapidly, thereby homogenizing the existing thermal stratification.
Here we have presented a solution for a class of flows in which periodic thermal forcing over non-uniform bathymetry drives a net advective heat flux, resulting in persistent environmental gradients. While the parameter regime considered theoretically here cannot generally represent realistic environmental parameters, the generation of persistent gradients arising from the mechanisms highlighted in the analytical solution can be expected to be generic for natural flows where temperature and scalar fields will have vertical structures that differ from that of the diurnally varying velocity field. Spatial gradients in the physical environment, such as water temperature, are important in determining the community structure of near-shore coastal ecosystems (Blanchette et al. Reference Blanchette, Miner, Raimondi, Lohse, Heady and Broitman2008; Tapia et al. Reference Tapia, Navarrete, Castillo, Menge, Castilla, Largier, Wieters, Broitman and Barth2009). For example, in coral reef ecosystems, persistent spatial differences in mean temperature and daily temperature variability have been shown to be among the most important factors determining resilience to thermal stress and bleaching (McClanahan et al. Reference McClanahan, Maina, Moothien-Pillay and Baker2005; Carilli, Donner & Hartmann Reference Carilli, Donner and Hartmann2012; Safaie et al. Reference Safaie, Silbiger, McClanahan, Pawlak, Barshis, Hench, Rogers, Williams and Davis2018).
The asymptotic analysis presented here assumes small values of the bottom slope 
$\beta \ll 1$. Effects of advective transport are found to be negligible at leading order, provided that 
$\beta ^2 {Ra} \ll 1$, so that our analysis is strictly valid for flow configurations with associated Rayleigh numbers 
${Ra} \ll \beta ^{-2}$. For environmental scales, estimates for 
${Ra}$ using typical values (cf. Monismith et al. Reference Monismith, Genin, Reidenbach, Yahel and Koseff2006; Molina et al. Reference Molina, Pawlak, Wells, Monismith and Merrifield2014) for cooling heat fluxes (200–500 W m
$^{-2}$), along-shore flow velocities (
${\sim }0.1$ m s
$^{-1}$), drag coefficients (
$C_D \sim 10^{-2}$) and using a depth-averaged eddy viscosity (cf. Fischer et al. Reference Fischer, List, Koh, Imberger and Brooks1979), yield values of 
${Ra}$ of a few thousand at depths of 10 m. For these flows, therefore, the advective-free flow predicted here applies only in situations where the slope is of order 
$\beta \sim 0.01$, whereas for larger values of 
$\beta$ one should expect significant nonlinear effects. Given this, however, it is worth mentioning that observations and analysis by Molina et al. (Reference Molina, Pawlak, Wells, Monismith and Merrifield2014) for a moderately shallow slope 
$\beta \sim 0.03\text {--}0.04$ at depths of 10 m show that a linear diffusive regime persists in configurations with 
$\beta ^2 {Ra} = O(1)$, and that the emergence of significant nonlinear advective effects requires larger values of 
$\beta ^2 {Ra}$, associated with steeper slopes (Monismith et al. Reference Monismith, Genin, Reidenbach, Yahel and Koseff2006). The persistence of the linear regime in environmental settings with relatively large 
${Ra}$ illustrates the relevance of the analysis presented here.
The applicability of the theoretical analysis developed here to realistic flows is further limited by the use of constant eddy viscosities and diffusivities. For environmental scale flows, both can be expected to vary in space and time. Because purely cross-shore flows are generally weak, directly forced turbulence and the associated turbulent transport can also generally be expected to be weak. In most coastal settings, however, turbulence will be dominated by wind- and tide-driven processes that operate on different time scales so that eddy diffusivities will be notably larger and can thus modulate the resulting cross-shore flow (Molina et al. Reference Molina, Pawlak, Wells, Monismith and Merrifield2014; Ulloa et al. Reference Ulloa, Davis, Monismith and Pawlak2018).
The analytical solution outlined here does not address any potential instability arising from inverted thermal profiles that occur during the cooling phase, consistent with the assumption of 
${Ra} \sim 1$. Numerical results for 
$\text {Pr}=1$ confirm that the periodic solutions are stable for 
${Ra}\le 13$ but become unstable for 
${Ra} \ge 14$. The ensuing dynamics is consistent with the onset of a double-period bifurcation at a critical value of the Rayleigh number 
$13 < {Ra}_c < 14$, slightly larger than the critical value 
${Ra}_c \simeq 12.44$ predicted by Floquet linear stability analyses in a semi-infinite domain (T. Christopher, personal communication), with the observed difference being attributable to the stabilizing effect of the wedge bottom surface. For the supercritical values of 
${Ra}$ considered in the simulations, the near-shore flow structure remains largely unaffected by the instability offshore, although we can anticipate unstable flow closer to shore with increasing 
${Ra}$, in agreement with previous findings by Mao et al. (Reference Mao, Lei and Patterson2010). For the observed bifurcation mode, unstable flow can be expected to effectively homogenize the temperature field, thus reducing the net advective heat flux that drives the residual circulation.
At the large values for 
${Ra}$ associated with environmental scales we can further anticipate that nightly surface cooling will drive convective instability over time scales that are much shorter than the diurnal cycle with vertical mixing reducing thermal stratification (Sevadjian, McManus & Pawlak Reference Sevadjian, McManus and Pawlak2010; Molina et al. Reference Molina, Pawlak, Wells, Monismith and Merrifield2014). Advective thermal fluxes will thus be weaker during this phase. However, a cooling phase response that persists into the stabilizing heating phase, as observed here and in field observations (Molina et al. Reference Molina, Pawlak, Wells, Monismith and Merrifield2014; Ulloa et al. Reference Ulloa, Davis, Monismith and Pawlak2018) will result in an associated onshore heat flux. The stabilizing heating phase will tend towards increasing vertical stratification with a vertical phasing that will vary relative to the velocity structure in a similar manner to that described by the theory. Cross-shore variations in phasing can then lead to divergences in heat flux with associated cross-shore thermal gradients.
The numerical simulations described in § 5 required initialization using a solution with the correct energetic balance in order to achieve convergence within a practical time frame. This is consistent with the magnitudes for the function 
$\mathcal {T}$ and the advective heat flux 
$\langle G_{ex}\rangle (x)$ shown in figures 4(a) and 6 which indicate that the time scale for establishing the steady temperature distribution will be 
$t \sim {Ra}$. The long time scale evolution of the perturbed temperature field associated with the unstable cases was not examined, although we can anticipate that the mean temperature field will be modified accordingly with an associated long transient for adjustment. Though the direct analogy with realistic, turbulent environmental scale flows is tenuous, the long adjustment time scale suggests that high variability observed in real, coastal situations (Molina et al. Reference Molina, Pawlak, Wells, Monismith and Merrifield2014; Ulloa et al. Reference Ulloa, Davis, Monismith and Pawlak2018) would result in a slowly varying residual.
The Eulerian residual flow that arises here at leading order has further consequences for cross-shore exchange and dispersion. In particular, at higher 
${Ra}$, the closed streamlines in the cross-shore cells associated with the streaming motion dominate the structure (e.g. figure 9) potentially introducing barriers to cross-shore flow. Convergences in surface flow (e.g. figure 8f) can also be conducive to frontogenesis with important implications for larval transport (Scotti & Pineda Reference Scotti and Pineda2007). Lagrangian transport associated with the non-uniform harmonic flow was not examined here, but can also be expected to play a role in dispersion.
Acknowledgements
Comments from three anonymous reviewers were very helpful in improving the paper.
Funding
The work here received support from the National Science Foundation via awards OCE-1436254 and OCE-1436522.
Declaration of interests
The authors report no conflict of interest.
Open access
