Coupling of waves to sea surface currents via horizontal density gradients

The mathematical models and numerical simulations reported here are motivated by satellite observations of horizontal gradients of sea surface temperature and salinity that are closely coordinated with the slowly varying envelope of the rapidly oscillating waves. This coordination of gradients of fluid material properties with wave envelopes tends to occur when strong horizontal buoyancy gradients are present. The nonlinear models of this coordinated movement presented here may provide future opportunities for the optimal design of satellite imagery that could simultaneously capture the dynamics of both waves and currents directly. The model derived here appears in two levels of approximation: first for rapidly oscillating waves, and then for their slowly varying envelope (SVE) approximation obtained by using the WKB approach. The WKB wave-current-buoyancy interaction model derived here for a free surface with significant horizontal buoyancy gradients indicates that the mechanism for the emergence of these correlations is the ponderomotive force of the slowly varying envelope of rapidly oscillating waves acting on the surface currents via the horizontal buoyancy gradient. In this model, the buoyancy gradient appears explicitly in the WKB wave momentum, which in turn generates density-weighted potential vorticity whenever the buoyancy gradient is not aligned with the wave-envelope gradient.


Submesoscale sea surface dynamics
Capabilities in sea surface observation have been improving rapidly during the past two decades [1]. In particular, new high-resolution satellite observation capabilities are revealing sea surface features seen for the first time at submesoscale spatial scales of 100 m -10 km and time scales of hours to weeks. Invariably, the new satellite imagery reveals a plethora of coupled dynamical surface phenomena, including currents, spiral filaments, flotsam patterns, jets and fronts, some of which are detected indirectly through gradients of sea surface temperature, salinity or colour, in addition to the imagery [5,10,13,20,26].
The new capabilities in sea surface observation are still developing. For example, the impending Surface Water Ocean Topography (SWOT) mission will map the ocean surface mesoscale sea surface height field, as well as a large fraction of the associated submesoscale field, including buoyancy fronts [17]. A sample of this type of submesoscale data taken from [5] is shown in figures 1 and 2.
The coming new age of higher-resolution upper ocean observations will present a formidable array of challenges for the next generation in data management, computational simulation and mathematical modelling. This paper will offer a mathematical modelling framework that is flexible enough to admit uncertainty quantification through stochastic modelling and analysis, applied in concert with highresolution observations, computational simulations, and stochastic data assimilation for large data sets. This framework involves decomposing the surface motion into a two-dimensional horizontal flow map 2 SUBMESOSCALE THERMAL WAVE-CURRENT DYNAMICS representing transport by the current acting on a one-dimensional vertical flow map representing wavelike motion of the elevation. This composition-of-maps modelling framework is described and applied to model sea-surface dynamics in two deterministic examples in section 2 of the present paper. Figure 1: Wave activity in the submesoscale ocean is dynamically complex, as illustrated in this figure showing the zoomed image of a submesoscale sea surface elevation, seen in Envisar MERIS glitter observations. This image shows the wave elevation tracking a cyclonic eddy visible in the sea surface glitter observations. The pixel resolution is 250m. This glitter image demonstrates the complex, highly-coordinated dynamical forms taken in wave-current interaction on the submesoscale sea surface. In particular, notice the instabilities developing in the eddy's outer boundary. Image courtesy of B. Chapron.
Emergent coherence (EC). Combining high-resolution thermal data (buoyancy) with glitter data for the wave elevation as in figure 2 has recently revealed yet another interesting feature of submesoscale dynamics. Namely, the observed submesoscale data show extremely high correlations of wave, current and thermal properties [5]. This emergent spatial-temporal coherence of dynamic and thermal properties presents a significant challenge for dynamical submesoscale modelling. Accepting this challenge, the aim of this paper is to derive a mathematical model of nonlinear sea surface dynamics whose solutions also demonstrate the emergent coherence observed in combining different types of submesoscale data. This paper derives new two-dimensional equations that show the emergent coherence (EC) seen in the sea surface features appearing in figure 2. The EC behaviour produced by the equations derived here are demonstrated in figure 3 which shows a snapshot of the coherence of buoyancy and wave amplitude distributions in the dynamics of divergence-free two-dimensional flow acting on free surface vertical elevation wave features moving under gravity. In the model equations, the horizontal buoyancy gradients mediate the interactions between the vertical elevation waves and the horizontal currents. The equations of motion represent the current as a time-dependent, area-preserving map of the horizontal plane into itself and the waves as the composition of the horizontal flow map with a time-dependent vertical elevation map. Thus, the model involves a dynamical composition of maps (C˝M).
2 Submesoscale thermal wave-current dynamics on a free surface 2.1 Surface waves as symmetry-breaking features of local force imbalances Waves are propagating symmetry-breaking features that signify the response to a local imbalance of forces. Thus, from the viewpoint of satellite oceanography, observations of waves -defined as propagating sea surface elevation features -signify processes at the surface or below the surface 2.  Figure 2: Comparison of the two images above demonstrates the emergent coherence between sea surface temperature and the glitter patterns visible from satellite imagery. The thermal fronts visible are dynamic, and sea surface roughness is most obvious along the strongest fronts. Discussions of the interpretation of sun glitter measurements are given in [5,20,26]. Images courtesy of B. Chapron. whose presence introduces forces that locally break the symmetry of the surface. The sea surface would otherwise follow the stable global gravitational balance of the geoid, which we regard here as being spherical. Thus, waves arise from a spatially local imbalance of forces in the neighbourhood of a stable equilibrium. The propagating feature of relevance here is the wave elevation, measured as the local departure of the surface level in the direction normal to its equilibrium mean level. The symmetry broken here is the invariance of the sea surface under spatial translations tangent to the equilibrium surface level, also known as the local horizontal direction. Hence, from the viewpoint of satellite oceanography, sea surface waves are observed as local vertical displacements of the otherwise horizontal motion of the ocean currents on the sea surface. From the mathematical modelling viewpoint, sea surface waves are local vertical oscillations of the horizontal surface that are carried along by the horizontal current flow, envisaged as a smooth invertible time-dependent map of the horizontal surface into itself. This is the composition of maps (C˝M) modelling approach for describing the dynamics of horizontal fluid flows (currents) acting on oscillating vertical elevations (waves). Since the surface current velocity, its advected material properties and the wave elevation are all that can be observed in satellite oceanography, the task in three-dimensional ocean modelling for satellite oceanography devolves into determining the dynamical surface features that are wave phase (bottom right) in the numerical simulation of the dynamics of divergence-free flow on a free surface moving under gravity. The simulation began with a spin-up period with zero wave amplitude. After the spin-up period, as explained in section 3, a checker-board pattern of finite wave amplitude with zero phase was introduced and the simulation was resumed. The 'mixing' of these wave patterns eventually brought them into coherence with the spatial distributions of thermal properties and potential vorticity. These features show an emergent coherence in patterns similar to those seen in the corresponding high-resolution satellite data in figure 2.
produced by the three-dimensional flow processes below the surface arising from e.g., bathymetry, stratification, rotation, Langmuir circulation, and thermal effects such as frontogenesis. The dynamics of the surface signatures of these three-dimensional flow processes, as well as effects of air-sea interactions on the surface need to be interpreted, in order to interpret what satellite oceanography observes. interacting at widely separated space-time scales. In this composition of maps (C˝M) approach, the waves are regarded as local vertical disturbances that rapidly oscillate as they are swept along by the broad, slowly changing horizontal currents. Thus, the slow current motion is a Lagrangian coordinate for the rapid wave oscillations. This wide separation in space-time scales invokes the classical WKB description. The standard WKB approach seeks a rapidly oscillating wave packet solution whose phase-averaged amplitude possesses a slowly varying envelope (SVE) spatially. The WKB method is often applied via a variational principle because in a variational setting the phase average naturally leads to an adiabatic invariant known as the wave action density, cf. for example, [3] for a review of the WKB or SVE method in fluid dynamics. Here we will follow the variational approach of [4,11] guided by the classical work of [22,24,25].

A tale of two maps: currents and waves
Submesoscale sea-surface motion: Composition of two time-dependent maps. The position and velocity of fluid parcels in motion under gravity on a 2D free surface embedded in R 3 have both horizontal and vertical components. The corresponding flow maps are denoted as the map φ t : R 2 Ñ R 2 for the horizontal current flow, and as the composite map ζ t φ t for the vertical elevation of the waves as a function of time and position in R 2 . The flow lines of these two components of the flow map of a free surface can be written as r t " φ t r 0 and z t " ζ t pφ t r 0 q ": ζ t pr t q , where r t " px t , y t q P R 2 is the horizontal position along the flow at time t and ζ t pr t q is the vertical elevation at horizontal position r t at time t, starting at position r 0 at time t " 0. Thus, one may say that the initial position of the flow line, r 0 , is a Lagrangian coordinate for the horizontal motion, and the horizontal motion is a Lagrangian coordinate for the vertical motion. That is, the 'footpoint' at time t of the vertical component of the flow map ζ t is located in the horizontal plane along a curve φ t r 0 parameterised by time t. Likewise, one can simply say that the wave dynamics is advected, or swept along, by the current dynamics.
Hence, the corresponding horizontal and vertical components of velocity along a stream line r t in That is, in the dynamics of free surface flow, the vertical velocity p wpr, tq at a given Eulerian point r and time t is related to the wave elevation ζpr, tq and horizontal velocity p vpr, tq at that point by p wpr, tq " B t ζpr, tq`p vpr, tq¨∇ r ζpr, tq .
In terms of these fluid variables, one could propose a Hamilton's principle for wave-current interaction of a free surface by following [8] for the variational modelling framework and applying [24,7] for the potential energy to find 1 To interpret the variational principle proposed in (2.1) we rewrite its Lagrangian as a sum of an Eulerian spatial integral and an integral over material mass elements d 2 r 0 " Dρ d 2 r which follow the paths of the horizontal fluid motion, rpr 0 , tq " φ t r 0 , Variations of the first summand in (2.2) at fixed spatial position prq yield the Euler fluid equations for 2D divergence free flow with advected buoyancy, ρpr, tq " ρpφ t r 0 q " ρ 0 pr 0 q, Variations of the second summand in (2.2) taken at fixed mass element pr 0 q yield equations for vertical harmonic oscillations of the elevation of each material mass element The wave-elevation equation in (2.4) is unrealistic, though, because it implies that fluid mass elements with different labels pr 0 q would be oscillating in phase and all with the same frequency, as they follow the flow of the Euler fluid equations (2.3) for 2D divergence free flow with advected buoyancy. This unrealistic synchronisation and resonance can be removed by including the inertia of each mass element. This can done by including the initial buoyancy of each mass element, as At this point in our reasoning, we have not yet considered the differences in space and time scales between the fluid flow and the wave activity. In what follows, we will use the simple compositionof-maps idea explained here along with estimates of relative space and time scales to investigate the applicability of this class of models. To improve the applicability of the model comprising (2.3) and (2.5) for describing the effects of currents on waves, we will derive a related model in the slowly varying envelope (SVE) approximation. The SVE approximation allows considerations of current and wave dynamics at the same space and time scales. The comparisons of the simulated solutions of these C˝M models with the observations in figures 1-4 above indicate that these models can indeed produce results that match some aspects of observed features. However, these models are not derived from three dimensional fluid equations. Instead, they are derived from the simple solution ansatz in Hamilton's principle that the vertical elevation of the sea surface wave activity is carried by divergence-free horizontal fluid motion. The latter assumption is a weakness of the current approach, because it precludes effects of vertical up-welling and down-welling, which are observed to occur along with convergence and divergence of currents [?]. The equations derived here are also not associated with classical surface wave equations such as the nonlinear Schroedinger (NLS) equation, or other celebrated surface wave equations. This departure from the classical water wave literature may be regarded as another weakness of the current approach.
Estimating parameters σ 2 and F r 2 for satellite observations. The Lagrangian pp v, ζ, D, ρq in (2.1) represents the dimension-free difference of the kinetic and potential energies, augmented by the incompressibility constraint imposed by the Lagrange multiplier p. Two dimension-free parameters (σ 2 and F r 2 ) appear in this Hamilton's principle. The coefficient σ 2 " prHs{rLsq 2 in formula (2.1) is the square of the vertical-to-horizontal aspect ratio. Typically, for satellite observations of submesoscale dynamics one finds rHs « p3ˆ10´4´3ˆ10´3qkm and rLs « p10´1´10qkm , so σ 2 « 10´3´10´6 ! 1 for the squared aspect ratio σ 2 ! 1 of the height of the waves rHs relative to the breadth rLs of the two-dimensional domain. The squared 'Froude number' F r 2 in this regime is estimated by the square of the ratio of horizontal and vertical frequency scales at the sea surface, Here, the horizontal velocity on the sea surface is taken as rV s " p0.1´1qm{sec, rHs " p0.3´3qm. According to [9], the Brunt-Väisälä buoyancy frequency in the sea surface wave regime is given by N « p10´3´10´4q{sec. The ratio of horizontal and vertical frequency scales at the sea surface in (2.6) is selected for use later in applying the slowly varying envelope (SVE) wave approximation in section 2.4. Hence, we estimate that the squared product of the 'Froude number' and aspect ratio for satellite observations of the sea surface can reasonably be estimated over the range Modelling the dynamic effects of surface density variations. As mentioned earlier, the observed oscillations of sea surface waves are by no means simultaneous across the whole domain, although the observations show that they are indeed coordinated spatially with the buoyancy of the fluid. To correct this solution behaviour, the kinetic energy and potential energy need to be de-synchronised from the buoyancy. The dynamic dependence of the wave kinetic energy on the density is physically required. However, to de-synchronise the wave oscillations we can introduce a constant reference density ρ ref into the wave potential energy, by writing The advected quantities Dpr, tqd 2 r and ρpr, tq evolve via push-forward by the horizontal flow map, φ t . For example, D t d 2 r t " φ t˚p D 0 d 2 r 0 q and ρ t " φ t˚ρref denote, respectively, evolution of the determinant of the Lagrange to Euler map and of the local scalar value of the mass density. Conservation of mass is then expressed as the push-forward relation, D t ρ t d 2 r t " φ t˚p D 0 ρ ref d 2 r 0 q. The pressure p in (2.9) acts as a Lagrange multiplier to enforce conservation of area, so that D t " 1 " φ t˚D0 , and the horizontal flow is incompressible, which implies that the horizontal velocity is divergence-free, i.e., div r p vpr, tq " 0. Taking variations of the action integral (2.9) yields the following set of equations, (2.10) From their definitions as advected quantities, one also knows that D and ρ satisfy where L p v denotes the Lie derivative operation along the horizontal velocity vector field, p v, which provides coordinate-free brevity in the notation. Proof. The Euler-Poincaré (EP) theorem in this case yields Here the diamond p˛q operator is defined by (2.14) In addition, X P X is a (smooth) vector field defined on R 2 and a P V , a vector space of advected quantities, which are here the scalar function, ρ, and the areal density D d 2 r. Using the advection relations for D and ρ in (2.11) and the corresponding variational derivatives in (2.10) simplifies the EP equation in (2.14) to Equation (2.10) then yields pB t`Lp v q`p v¨dr`σ 2 p w dζ˘"´ρ´1dr p`d r . Inserting the last relation into the following standard relation for the time derivative of a loop integral then completes the proof of equation (2.12) appearing in the statement of the theorem, dr p`d r . (2.16) Using the advection relations for D and ρ in (2.11) again and combining with the variational relations with respect to ζ in (2.10) simplifies the p w and ζ equations, as follows.
pB t`Lp v qζ " pB t`p v¨∇ r qζ " p w . (2.17) After deriving these equations, one may finally evaluate the constraint D " 1 imposed by the variation in pressure p to obtain further simplifications. , was also observed in [8]. This behaviour is consistent with the Charney-Drazin 'non-acceleration' theorem [6,23]. Namely, in certain circumstances, wave activity does not create circulation in the mean current. A modification that allows exchange of circulation between wave (vertical) and current (horizontal) components of the flow was proposed in [8]. The instabilities observed around the edges of eddies in the satellite imagery shown in figure 1 suggests that a coupling of this sort may exist at high wave number.

Thermal potential vorticity (TPV) dynamics on a free surface
The momentum map arising from the variations in (2.10) is given by As expected from the well-known non-acceleration theorem [6,23], the dynamics of the Euler-Poincaré equations separate (2.15) gives the dynamics of the fluid and wave components of the momentum oneform (2.20). (2.21) The mass-weighted thermal potential vorticity (TPV) also separates into fluid and wave components Q " Q F`QW with following definitions where buoyancy weighted vertical velocity is defined as r w :" ρ p w. The dynamics of of Q F d 2 r and Q W d 2 r can be computed from (2.21) as (2.23) From the two relations in (2.23), one sees that the buoyancy gradient ∇ρ couples the PV dynamics of the waves pQ W q and currents pQ F q, each to their corresponding kinetic energy. In the case of constant buoyancy, dρ " 0 in (2.23); so, the PVs of the waves and currents would be separately advected. The operator divpρ∇q is invertible, so long as ρ is a differentiable positive function, which can be ensured by requiring that this condition holds initially. Consequently, the stream function ψ is related to the other fluid variables by ψ :" pdivρ∇q´1Q F . (2.24) The potential vorticity dynamics can then be written in coordinate form as  (2.26) The energy Hamiltonian hpQ, ρ, p w, ζq associated with this system is given by hpQ, ρ, r w, ζq " Proof. The Casimirs C Φ,Ψ for the direct sum of the Lie-Poisson brackets for Q and ρ and canonical Poisson brackets for r w and ζ follows by direct verification that the C Φ,Ψ are conserved for any differentiable functions, pΦ, Ψq.

C˝M equations in the slowly varying envelope (SVE) approximation
The SVE solutions apply to satellite observations of sea surface waves. From the viewpoint of satellite observations, the vertical motion on the sea surface typically oscillates much more quickly than the rate of change of features in the horizontal motion of the ocean surface currents. In this situation, the standard WKB approximation introduces a solution Ansatz for the slowly varying envelope (SVE) of the rapidly oscillating vertical wave elevation in the standard form [2,11], ζpr, tq " ˆa pr, tq exp´i θpr, tq ¯˙w ith ! 1 . (2.29) The SVE solution Ansatz (2.29) comprises the product of a slowly varying complex amplitude apr, tq P C multiplied by a rapidly oscillating phase θpr, tq{ P R with ! 1 in which the phase factor θpr, tq may also vary slowly as a function of the space and time variables, pr, tq. Following [11], let us substitute the SVE solution Ansatz (2.29) into Hamilton's principle in (2.9) and find the condition on the parameter ! 1 that will allow higher order wave terms to be neglected. For this, one computes (2.30) 12

C˝M EQUATIONS IN THE SVE APPROXIMATION
The leading order wave term Op ´2 q with ! 1 in Hamilton's principle will dominate the solution and the remaining wave terms in the second line of equation (2.31) may be neglected, when 2 ! 1 , 2 σ 2 F r 2 " Op1q, and σ 2 F r 2 ! 1 .

(2.31)
According to the estimates in (2.7) there is a range of physical parameters relevant to satellite observations in which the SVE approximation applies, for σ 2 F r 2 ! 1.
To continue the investigation of the SVE description of wave-current interactions on the sea surface, we take variations of the action integral (2.31) to find the following set of equations, with ωpr, tq "´B t θ and kpr, tq " ∇ r θ , (2.32) In the second line of (2.32) we see that stationarity of the action integral with respect to variations in |a| 2 acts as a Lagrange multiplier to impose a constraint which relates the dynamics of the wave phase θ to the buoyancy. This constraint relation involves the Doppler-shifted frequency of the waves, as shown in the third line of (2.32). In combination with conservation of the wave action density and the divergence free condition on the fluid flow velocity p v, this constraint relation implies in the last line of (2.32) that the wave magnitude |a| 2 is advected by the fluid flow. Because of the oscillatory nature of the solution Ansatz (2.29), the sign of the wave phase in dθ{dt " B t θ`p v¨∇ r θ in the second line above is immaterial. Hence, hereafter, we will choose the positive root for dθ{dt " ? ρρ ref {ρ. From the conservation of wave action density A in (2.32) and the definitions of the advected fluid variables, one finds that |a| 2 , D and ρ satisfy the following advection relations pB t`Lp v qpD d 2 rq " 0 ùñ B t D`div r pDp vq " 0 with D " 1 , where L p v denotes the Lie derivative operation along the horizontal velocity vector field, p v. The Lie derivative notation L p v provides coordinate-free brevity in proving the following Kelvin circulation theorem for thermal wave-current theory. Proof. The Euler-Poincaré (EP) theorem [16] in this case yields Here, the diamond p˛q operator is defined for a fluid advected quantity f by (2.36) In (2.36), X P XpR 2 q is a (smooth) vector field defined on R 2 and f P V is a vector space of advected quantities. These advected quantities are the scalar function, ρ, and the areal density, D d 2 r. Upon using the advection relations for D and ρ in (2.33) and the corresponding variational derivatives in (2.32), the EP equation in (2.35) simplifies to Equation (2.32) then yields pB t`Lp v q´p v¨dr`N d dθ dt¯"´ρ´1 dp`d´1 2 |p v| 2¯. (2.37) Inserting the last relation into the following standard relation for the time derivative of a loop integral then completes the proof of equation (2.34) appearing in the statement of the theorem, Note, however, that equations (2.32) imply the following combination of advected quantities, Consequently, the wave-momentum 1-form N dp dθ dt q is advected by the fluid flow and the Kelvin circulation theorem in equation (2.38) reduces to the standard circulation theorem for the 2D Euler fluid equations. Remark 2.3 (Separation of wave and current motion in the SVE approximation). The decoupling of the Kelvin-Noether circulation theorem into its wave and current components for the SVE approximation is inherited from the un-approximated model. When modifications to the un-approximated model which removes this property are added, one would expect the new SVE approximation to lose the non-acceleration result.
Remark 2.4. Equation (2.39) implies advection of the 1-form |a| 2 dρ, which in turn implies advection of the Jacobian Jp|a| 2 , ρq. Since the fluid flow is area preserving, divp v " 0, the following 2-form will also be advected,`B t`p v¨∇ r˘´d |a| 2^d ρ¯" 0 . (2.40) Thus, the divergence-free flow of p v preserves the area element d|a| 2^d ρ. This means that if the gradients ∇|a| 2 and ∇ρ are not aligned initially, then they will remain so. It also means that equilibrium solutions of (2.40) will be symplectic manifolds [14].
After deriving these equations, one may finally evaluate the constraint D " 1 imposed by the variation in pressure p to obtain further simplifications.

Thermal potential vorticity dynamics with SVE on a free surface
The momentum map arising from the variations of the action in (2.32) is given by This means that the mass-weighted thermal potential vorticity (TPV) dynamics also separates into the following fluid and wave components, Q " Q F`QW , given by with Q F :" divpρ∇ψq and Q W :" ΓJ`?ρρ ref , |a| 2˘. (2.43) Then, again, the differentials of the separate equations in (2.42) yield the 'non-acceleration' result, Equivalently, in coordinates one has with Q F :" divpρ∇ψq and Q W :" ΓJ`?ρρ ref , |a| 2˘, B t ρ`p v¨∇ r ρ " 0 and Γ " σ 2 4F r 2 " Op1q , B t |a| 2`p v¨∇ r |a| 2 " 0 , (2.45) The operator pdivρ∇q is invertible, so long as ρ is a differentiable positive function, which can be ensured by requiring that this condition holds initially, since ρ is advected. Consequently, the stream function ψ is related to the other fluid variables by ψ :" pdivρ∇q´1Q F .

Numerical implementation
Our implementation of the C˝M equations (2.25) and the C˝M equations in the SVE approximation (2.45) used the finite element method (FEM) for the spatial variables. The FEM algorithm we used is based on the algorithm formulated in [15] and is implemented using the Firedrake 3 software. In particular, for (2.25) we approximated the fluid potential vorticity Q F , buoyancy ρ, wave elevation ζ and bouyancy weighted wave vertical velocityw using a first order discrete Galerkin finite element space. Similarly, for (2.45), we approximated Q F , ρ, square of the wave amplitude |a| 2 and wave phase θ using a first order discrete Galerkin finite element space. The stream function ψ for both models was approximated by using a first order continuous Galerkin finite element space. For the time integration, we used the third order strong stability preserving Runge Kutta method [12]. Figures 3 and 4 present snapshots of high resolution runs of the C˝M equations and the C˝M equations in the SVE approximation. These simulations were run with the following parameters. The domain is r0, 1s 2 at a resolution of 512 2 . The boundary conditions are periodic in the x direction, and homogeneous Dirichlet for ψ in the y direction. To see the effects of the waves on the currents, the procedure was divided into two stages for both set of equations. The first stage was performed without wave activity for T spin " 100 time units starting from the following initial conditions Q F px, y, 0q " sinp8πxq sinp8πyq`0.4 cosp6πxq cosp6πyq`0.3 cosp10πxq cosp4πyq0 .02 sinp2πyq`0.02 sinp2πxq , ρpx, y, 0q " 1`0.2 sinp2πxq sinp2πyq and ρ ref " 1 . The purpose of the first stage was to allow the system to spin up to a statistically steady state without any wave activity. The PV and buoyancy variables at the end of the initial spin-up period are denoted as Q spin px, yq " Q F px, y, T spin q and ρ spin px, yq " ρpx, y, T spin q. initial conditions for the flow variables being the state achieved at the end of the first stage. To start the second stage for (2.25), wave variables were introduced with the following initial conditions ζpx, y, 0q " sinp8πxq sinp8πyq`0.4 cosp6πxq cosp6πyq`0.3 cosp10πxq cosp4πyq0 .02 sinp2πyq`0.02 sinp2πxq , wpx, y, 0q " 0 , Q F px, y, 0q " Q spin px, yq , ρpx, y, 0q " ρ spin px, yq , To start the second stage for (2.45), wave variables were introduced with the following initial conditions |a| 2 px, y, 0q "`sinp8πxq sinp8πyq`0.4 cosp6πxq cosp6πyq`0.3 cosp10πxq cosp4πyq0 .02 sinp2πyq`0.02 sinp2πxq˘2 , θpx, y, 0q " 0 , Q F px, y, 0q " Q spin px, yq , ρpx, y, 0q " ρ spin px, yq .

Conclusion and Outlook
This paper models the effects of thermal fronts on the dynamics of the ocean's waves and currents. It introduces and simulates two models of thermal wave-current dynamics on a free surface. The original C˝M model is derived from Hamilton's principle via the composition of two maps which represent the horizontal and vertical motion respectively. The second, a slowly varying encelope (SVE) model, is introduced via the standard WKB approximation which takes advantage of large separation of the space-time scales between the slow horizontal currents and fast vertical oscillations. In particular, the second model introduces the WKB solution Ansatz into Hamilton's principle, whereupon the time integral averages over the phases of the rapid oscillations that are out of resonance with the slowly varying envelope. Model runs of both models are presented in which the buoyancy mediates the dynamics of the currents and waves, as seen in Figures 3 and 4. These simulations also validate the use of the WKB approximation for two reasons. First, the resolved small scale wave features of the original C˝M model lie primarily within the envelope defined by the SVE approximate model. This means that the dynamics of the spatial features of the SVE approximate model are consistent with those of the original C˝M model, although the resolved space and time scales differ. Secondly, requiring that 2 {pF r 2 σ 2 q " Op1q ensures that the time scale for the wave envelope dynamics matches that for the fluid motion.
Nonetheless, the two models introduced here merit further study in several directions. For example, it remains to: (1) quantify the correlations observed visually; (2) determine their rate of formation; and (3) parameterise the model for comparison and analysis of the satellite data on which their derivations were based. Furthermore, the models discussed here involve only variables that are evaluated on the free surface and therefore they neglect bathymetry. A scientific challenge persists in understanding regions of the ocean where bathymetry has profound effects on the observable surface dynamics, such as in the Lofoten vortex [21]. This is a multiscale issue that might be addressed by including mesoscale modulations of the sub-mesoscale models derived here. One candidate for providing the mesoscale modulations would the thermal quasi-geostrophic (TQG) model in which bathymetry has recently been included [15].
The currents are modelled here by the two dimensional incompressible Euler equations, as seen in equations (2.2) and (2.3). Incompressibility is a reasonable assumption in some regions of the ocean, for example when the quasigeostrophic approximation is valid. There are regions in the upper ocean where other equations are more suitable for modelling currents, and the development and investigation of such two dimensional models is an open problem which warrants further consideration.
As mentioned in Remark 2.1, the wave component of the model presented here does not create circulation in the currents. The instabilities present in satellite simulations indicate that additional modelling is needed to fully capture this effect. Future work will investigate approaches for modelling these instabilities.
Many other questions remain about wave-current interaction. The full extent of submesoscale ocean dynamics is by no means adequately described by existing models. For example, we have little understanding of the formation and dynamics of various sea-surface phenomena, including the socalled 'spirals on the sea' [18]. Other questions are emerging because the ocean has absorbed in excess of 90% of the heat present in the earth system as a result of human activity during the post-industrial era [19]. The absorption of heat from the warming atmosphere is ongoing and it is forecast to become more dramatic. This absorption has resulted in 'marine heat waves', which are predicted to increase in frequency and severity. These changes to the upper ocean, where most of this heat is stored, could have a profound effect on the dynamical landscape of our oceans. These effects may, in turn, influence our weather and climate systems. Over the millennia, the ocean has approached statistical equilibrium under its current forcing conditions. Using modelling terminology, one says the ocean is well 'spunup'. However, the continued warming of the ocean is likely to influence the number and intensity of thermal fronts. One hopes that mathematical models will provide a useful framework for estimating some of the potential impacts of these thermal fronts on atmospheric effects, as well.