Tidal Forcing in Icy‐Satellite Oceans Drives Mean Circulation and Ice‐Shell Torques

Tidal forces generate time‐varying currents in bodies with fluid layers, such as the icy ocean moons of the outer solar system. The expectation has been that tidal currents are periodic—they average to zero over a forcing period—so that they are not associated with a mean flow. This expectation arises from the assumption of linearity. Here, we relax this assumption and develop a theory that predicts the emergence of mean currents driven by any periodic forcing. The theory, derived in the context of a global, uniform, shallow ocean, constitutes a set of mean flow equations forced by non‐linear eddy fluctuations. The latter are the canonical, periodic tidal currents predicted by the Laplace Tidal equations. We show that the degree‐2 tide‐raising potential due to obliquity and/or orbital eccentricity can drive time‐averaged currents with zonal wavenumbers from 0 to 4. The most prominent of these is a retrograde zonal jet driven by the obliquity‐forcing potential. Assuming Cassini state obliquities, this jet has speeds ranging from 0.01 to 1 mm s−1, which can exert torques up to roughly 1015 N m at the ice–ocean interfaces of Europa, Callisto, Titan, and Triton. Depending on the viscosity of the ice shell, these torques could drive ice shell drift rates of tens to potentially hundreds of meters a year. Thinner or stably stratified global oceans can experience much faster mean currents.


Introduction
Global oceans of liquid water are known to exist beneath the frozen exteriors of the outer solar system moons Europa, Ganymede, Callisto, Titan, Enceladus, potentially Dione, and even Mimas (Beuthe et al., 2016;Kivelson et al., 2000Kivelson et al., , 2002;;Lainey et al., 2024;Rhoden & Walker, 2022;Thomas et al., 2016;Zimmer et al., 2000).These bodies have attracted wide interest due to their potential to provide a habitat for life.On Earth, much of our climate and biology is controlled by time-mean ocean currents that transport tracers such as heat, salt, and nutrients.For icy moons, time-mean currents may be critical for the thermal evolution of the ice shell, and the chemical transport of oxidants and other astrobiologically useful products.A recently suggested, profound geophysical consequence of time-mean zonal (east-west) flows is their net frictional torque at the ice-ocean interface, which may alter the rotation rate of icy shells (Ashkenazy et al., 2023;Hay et al., 2023;Kang, 2024).This torque may drive non-synchronous rotation of the ice shell.Such surface reorientation has been suggested to play a role in the formation of surface tectonic features on Europa (e.g., Hurford et al., 2007;Rhoden & Hurford, 2013;Rhoden et al., 2010) and in altering the hemispherical distribution of its impact craters (Zahnle et al., 2001).
Large-scale flows are theorized to develop within icy-moon oceans due to thermal convection (e.g., Soderlund et al., 2014), tidal forcing (e.g., Tyler, 2008), and magnetic pumping (Gissinger & Petitdemange, 2019).For synchronously rotating moons, ocean tides can become important when the body has an elliptical orbit and/or a tilted rotation axis, both of which induce large-scale oscillating currents driven by periodic variations in the gravitational potential across the body.In this manuscript we focus on tidally driven flow, and the capacity for tidal forcing to drive time-averaged ocean circulation within the internal oceans of outer solar system moons.A fundamental question is whether a perfectly periodic tidal forcing can generate a time-averaged flow field.We hypothesize that this is indeed the case due to nonlinear interactions within the oscillating portion of the tidal flow.We present a theoretical framework that describes this process, and we demonstrate the plausibility and potential geophysical consequences of this hypothesis.
Here we address this shortcoming of linearized tidal flow models by developing a weakly nonlinear theory of time-mean currents driven by tidal forcing.We do so in the context of classic shallow-water theory, where we vertically integrate the nonlinear momentum and mass conservation equations to solve for the depth-averaged velocity of the ocean, and the deformation of its surface.Using a Reynolds decomposition, we split the momentum and mass conservation equations into a time-mean system and an oscillating system.By assuming that the time-mean velocity and surface deformation are much smaller than their oscillating counterparts, we are able to partially decouple the mean flow from the oscillating tide.
Focusing on Europa, we use our theory to determine the spatial structure of time-mean flow fields that arise from eccentricity-and obliquity-forced tides.We find that obliquity forcing can drive retrograde jets with flow speeds of less than a mm s 1 , which is capable of creating ice-ocean torques of up to 10 15 N m.Our method is general and readily applied to any global fluid layer in the shallow-water limit.
The manuscript is organized as follows.In the next section we develop a theoretical framework to calculate tidally driven time-mean flows.The dominant time-mean tidal currents on Europa are presented in Section 3 for endmember ocean thicknesses, alongside approximate analytical solutions for the strongest circulation and its sensitivity to drag.Implications of these time-mean currents, including the potential retrograde ice-ocean torques predicted by our theory, are discussed in Section 4. The paper is concluded in Section 5 with a summary of the main findings.

Theory
Here we give an overview of the LTE, the manipulation of the LTE into semi-coupled systems for the mean and periodic flows, the tidal forcing potentials, and numerical benchmarks of our theory.

The Laplace Tidal Equations
If the ocean is sufficiently thin in comparison to the radius at the top of the ocean, its barotropic dynamics can be described by the nonlinear shallow water equations.With the inclusion of a gravitational tidal forcing potential, U T , these are known as the LTE.The momentum of an ocean with homogenous density ρ, depth-averaged velocity u and thickness h is conserved via, Mass is correspondingly conserved by adjusting the thickness of the incompressible ocean column to compensate for divergence in lateral thickness flux, Here, the satellite is rotating with angular velocity Ω and has surface gravity g.The ocean surface is perturbed, relative to the initial, spherically uniform thickness h 0 , by a radial distance η, which is the tide.The total column thickness is therefore h = h 0 + η.Dissipation, the conversion of kinetic and/or potential energy to heat via friction, is accounted for via the operator D. The velocity field is tangential to the surface along the colatitude, e θ (positive southward), and longitude, e ϕ (positive eastward), unit vectors, u = ve θ + ue ϕ .We refer to v and u as the meridional/southward and zonal/eastward velocity components, respectively.The divergence and gradient operators are horizontal operators only.All differential operators used in this manuscript are defined in Appendix A.
The momentum equation is sometimes written in the vector-invariant form (Vallis, 2017), which introduces the relative vorticity of the ocean column, where e r is the radial unit vector, and f = 2 Ω cos θ is the Coriolis parameter that varies with co-latitude θ.We use both forms of the momentum equation in this manuscript.
Drag is typically modeled using one of the following three terms; (5) The first term is a simple linear drag approximation with proportionality constant α.It does not necessarily describe any specific dissipation process, though it has been used on Earth to parameterize the conversion from barotropic to baroclinic tides and the associated wave breaking (e.g., Jayne & St. Laurent, 2001).If icy satellite oceans are stratified, then this processes may be important (Abdulah et al., 2022;Idini & Nimmo, 2024;Rovira-Navarro et al., 2023).Prescribing drag in this linear way is common in tidal dissipation studies in icy satellite oceans (e.g., Matsuyama et al., 2018;Rovira-Navarro et al., 2020, 2023;Tyler, 2011).As α has units of inverse seconds, it can be thought of the inverse of the decay timescale of the momentum of an ocean.The second term accounts for viscous dissipation and is sometimes called Navier-Stokes drag (Chen et al., 2014).This can account for dissipation at the molecular level, or large-scale dissipation due to eddy diffusion.The eddy diffusivity for icy satellites is unknown, though it has been recently estimated to be ν ∼ 300 m 2 s 1 in Europa's ocean using eddypermitting numerical simulations by Ashkenazy and Tziperman (2021).The last term, quadratic or bottom friction, accounts for dissipation due to turbulence in the ocean top and bottom boundary layers.As it is nonlinear, numerical methods are typically required to include its effect (Hay & Matsuyama, 2017, 2019;Vincent et al., 2022).Chen et al. (2014) derived scaling laws to account for bottom drag dissipation, which were later extended to subsurface oceans by Hay and Matsuyama (2019).This is the preferred drag mechanism when attempting to make quantitative predictions of oceanic tidal heating because the dimensionless drag coefficient c D is well constrained on Earth.However, to make analytical progress, we consider only linear and Navier-Stokes drag in this manuscript.
We note that the momentum equations above have not been modified to include gravitational and mechanical iceocean coupling (Beuthe, 2016;Matsuyama et al., 2018), nor the corrections from deformation of the seafloor due to oceanic mass loading and gravitational self-attraction by the ocean (Matsuyama, 2014).These corrections can be important for calculating energy dissipation rates, particularly in smaller icy satellites with relatively thick ice shells such as Enceladus (Hay & Matsuyama, 2019), or when near resonance (Beuthe, 2016;Matsuyama, 2014).
However, the focus of this manuscript is on identifying the circulation patterns that can be driven by tidal forcing, so we neglect these effects.We discuss this simplification in Section 4.3.

Forcing-Generated Mean Flow
We consider tidal forcing that is perfectly periodic.In particular, the tidal forcing potential integrated over a forcing period is zero.Here we show, by making use of a Reynold's decomposition, that periodic forcing can generate a mean flow through nonlinear effects.
We start by splitting u into a forcing period-average u to describe the time-mean tidal currents, and a timefluctuation u′ to describe periodic tidal currents.Substituting this into the vector-invariant form of the momentum Equation 3, taking a time-average over the forcing period, and assuming that the growth timescale of the mean flow is much longer than the forcing period, gives Taking the radial component of the curl gives the time-mean vorticity balance, and taking the divergence gives, Both Equations 7 and 8 are written so that the left-hand terms depend only on the mean flow, while the right-hand side contains non-linear terms from the fluctuating flow.We call the terms on the right-hand side eddy forcing terms.
The continuity equation can be manipulated in a similar way, where we decompose the total ocean thickness as h = h 0 + η′ + η, where h 0 is the spatially uniform ocean thickness, η is the period-averaged tide, and η′ is the time-fluctuating tide.Note that η is a deformation of the ocean surface driven by ocean currents, and is not related to the time-independent portion of the tide-raising potential (which we anyway exclude).Once again taking a time-mean over the forcing period, and assuming that the growth timescale of the mean component is much greater than the tidal forcing period, Equation 2 can be written as, which again has an eddy forcing term on the right-hand side.
The time derivatives in Equations 7-9 are useful for interpreting our results later in the manuscript, but in the following we will now assume that the mean flow is in steady state.Next, we introduce a Helmholtz decomposition for the mean flow that splits it into spheroidal Φ) and toroidal Ψ) scalar components, where R is the mean radius of the body.Sometimes Φ is referred to as the poloidal or consoidal component.The Helmholtz decomposition is useful because, and Substituting Equation 10 into Equations 7-9 at steady state gives, Here we have made the critical assumption that the mean flow and mean tide are much smaller than their fluctuating counterparts, |u| ≪ |u′| and η ≪ η′, so that quadratic terms in mean quantities from Equations 7-9 can be neglected.Additionally, we assume that the ocean is much thicker than the time-mean tide, η ≪ h 0 .This yields a coupled linear system (Equation 13) of three equations for the three time-mean unknowns, the spheroidal velocity potential Φ, toroidal stream function Ψ, and time-mean ocean tide η.The system is forced through eddy transport by the fluctuation terms on the right-hand side.The last equation comes from continuity, and essentially tells us that the velocity potential exists solely to balance fluxes of ocean thickness η′u′) due to correlations in the fluctuating parts of the flow.Unlike in the classic LTE (i.e., Equation 20), the relationship between the time-mean tide and the velocity field is not straightforward and can be found only through force balance (Equation 13b), rather than through mass conservation.Equations 13a and 13b reduce to geostrophic balance if drag and the righthand sides are neglected.
Using the periodicity of the sphere, the unknown, time-mean quantities can be expressed in the form, where Y m n = P m n ( cos θ) e imϕ is the complex, unnormalized, degree n and order m spherical harmonic, P m n ( cos θ) is the unnormalized Associated Legendre Polynomial, Ψ m n is the corresponding complex spherical harmonic expansion coefficient, c.c. is the complex conjugate, and the superscripts s and t refer to spheroidal and toroidal vector harmonic components, respectively.The time-mean eddy transport vector fields, ζ′u′ and η′u′, are expanded in terms of the vector spherical harmonics These are the spheroidal and toroidal spherical harmonic vectors, respectively.We refer to all spherical harmonic expansion coefficients of the eddy terms (Z m,s n , Z m,t n , H m,s n , H m,t n , and J m n ) as eddy forcing coefficients.
Substituting Equations 14 and 15 into Equation 13 and following Longuet-Higgins (1968) (i.e., collecting coefficients of terms in P m n , cos(mϕ), and sin(mϕ)), gives a system of linear equations that is coupled across spherical harmonic degree n, and partially uncoupled in order m (the zonal wavenumber), where, Given some mean eddy fluxes of vorticity Z m,s n and Z m,t n , surface displacement H m,s n , and kinetic energy per unit mass J m n , the resulting secondary mean flow can be found by solving system Equation 16.This system is conveniently decoupled in m.However, the full system is not decoupled since the eddy forcing components on the right hand side depend on coupling between different m components of the fluctuating tide, as discussed below in Section 2.4.Indeed, the task now is to determine what Z m,s n , Z m,t n , H m,s n , and J m n are.Crucially, because we have assumed the mean flow/deformation to be much smaller than its fluctuating counterpart, the fluctuating tidal flow field remains accurately predicted through linear theory, provided that the forcing and linear response are small.

The Linearized Laplace Tidal Equations' Full Analytical Solution
Linearizing about a state of rest and uniform ocean thickness h 0 , the momentum Equation 1 for the fluctuating, depth-averaged velocity field u′, and continuity Equation 2 for the fluctuating ocean tide η′, become ∂u′ ∂t As with the time-mean system derived above, we decompose the fluctuating velocity field into solenoidal (Φ) and toroidal (Ψ) components using a Helmholtz decomposition, where we omit primes to ease the notation.This decomposition simplifies the linearized mass conservation in Equation 19 to, Taking the radial component of the curl of Equation 18gives, Equations 21 and 22 are equivalent to Equations 16 and 17 in Chen et al. (2014).
As we are dealing with periodic forcing in both space and time we assume solutions of the form, where Ψ m,X n is the corresponding spherical harmonic expansion coefficient of the fluctuating stream function, c.c. is the complex conjugate, and X is the direction of propagation.For eastward traveling waves (X = E) we have ω E = Ω, whereas for westward traveling waves (X = W) we have ω W = Ω.Equivalent expansions are taken for Φ, η′, and U T .Substituting this form into Equations 21 and 22, alongside the substantial algebra outlined in Longuet-Higgins (1968), gives an algebraic system of equations coupled in spherical harmonic degree n, but uncoupled in order m and propagation direction where, and p m n , q m n , and b n are defined in Equation 17.With the exceptions that we neglect ocean loading/self-attraction, deformation of the solid regions, ice-ocean coupling, and we have included Navier-Stokes drag, the linear algebraic system in Equation 24 is equivalent to Equations 21 and 22 in Chen et al. (2014), Equation 50in Beuthe (2016), Equation C.8 in Matsuyama et al. (2018), and Equations A1 and A2 in Hay et al. (2022).Equation 24 must be solved for each spherical harmonic degree n for each component of the perturbing potential at frequency ω, degree n, and order m, so as to recover the linearized, periodic, barotropic, ocean flow response to tidal forcing.These, together, can be used to find the eddy forcing coefficients in the mean flow system, Equation 16.To find the full semi-analytical solution, up to an arbitrary Nth degree of truncation, we use the code developed and made available in Hay et al. (2022).In what we refer to as the approximate analytical solution, we truncate the solution to only a few degrees such that it is analytically tractable by hand.

Calculating the Eddy Forcing
Now we have two sets of algebraic equations, one for the steady state mean flow (Equation 16) and one for the periodic flow (Equation 24).In order to solve for the mean flow, we need to first determine the eddy coefficients Z m,s n , Z m,t n , H m,s n , and J m n , that act as the forcing in the mean flow system.To do this, we first note that the fluctuating velocity field can be written in a form involving vector spherical harmonics, u′ = u s + u t , where which is functionally identical to substituting Equation 23 into the Helmholtz decomposition.
Due to the orthogonality of the spherical harmonics, the spheroidal and toroidal eddy forcing coefficients and kinetic energy coefficient can be found with, where dΩ = sin θdθdϕ is the solid angle and tildes represent a complex conjugate.The normalization factor N m n is defined as, We can calculate the eddy forcing coefficients numerically on a grid using Equation 29, though it is useful to find analytical expressions for them to understand which modes in the oscillating tide couple together to produce the eddy transport; these we provide in Appendix B. Note that we use the term mode to refer to an individual component of the tidal solutions, for example, Ψ m,X n or Φ m,X n for a specific n, m, and X.The key to determining the eddy forcing coefficients lies in calculating what is essentially a surface integral of the product of three spherical harmonics.Many combinations of different degrees and orders are trivially zero for such an integral (Arfken & Weber, 1999, p. 803).One particularly important rule is that the sum of all the ms in the integral must be zero so that the integrand is not periodic in longitude.By comparing how many complex conjugates are in each integral, and whether the coupling coefficient couples waves traveling in the same or opposite directions, we find that the zonal wavenumber m of the eddy forcing coefficients obeys, for waves traveling in the same direction, where m 1 and m 2 are the zonal wavenumbers of the interacting modes.A mode that interacts with only itself (|m 1 | = |m 2 |) will therefore only ever produce a zonally invariant flow (m = 0).If the waves travel in opposite directions, These rules can help to provide some intuition into the results that follow.

Tidal Forcing
We first consider periodic tidal forcing due to orbital eccentricity.This is followed in Section 2.5.2 by obliquity tidal forcing.

The Eccentricity Tide
If the satellite is in an elliptical orbit but has no obliquity, an observer on the surface will see the tide-raiser librate in longitude as well as change apparent size due to changes in planet-satellite distance.This generates the eccentricity tide, a triaxial distortion of the moon that is symmetric about the equator.In a reference frame rotating with the satellite and aligned along the permanent tidal bulge, the diurnal (|ω| = Ω) eccentricity tide-raising potential is (Tyler, 2011), which is accurate to first order in eccentricity e.In the second line, we have split this into contributions from an eastward E, westward W, and zonal Z (i.e., independent of longitude) component, where, ). (34)

The Obliquity Tide
An observer standing on the surface of a synchronously locked satellite with a nonzero obliquity and circular orbit will see the tide-raising object appear to librate in latitude.This generates the obliquity tide, an asymmetric tidal distortion of the moon that oscillates across the equator.In the same reference frame as before, the tidal forcing potential that generates this distortion is, where only the first-order terms in obliquity angle θ o have been retained (Tyler, 2011).Similarly to eccentricity forcing, the second line splits the potential into a westward and eastward traveling component, Only the westward component of the obliquity tide-raising potential can excite Rossby-Haurwitz waves in a global ocean (Tyler, 2008).

Summary of the Procedure
To summarize, the mean flow can be calculated using the following steps.First, the LTE system in Equation 24 is solved to find the fluctuating tide, η′, and its currents, u′, for each spherical harmonic order m that is present in the forcing potential (Equation 33 and/or Equation 35).All components of the fluctuating tide and currents are then used to determine the eddy transport vectors, ζ′u′ and η′u′, either numerically or by using Equation B1.The timeaveraged currents and tide, u and η, are then obtained by solving system Equation 16 with the forcing prescribed through the eddy transport coefficients.This last step must be repeated for every m component present in the eddy forcing.

Model Verification
We use the numerical model outlined in Hay and Matsuyama (2017) and Hay and Matsuyama (2019), Ocean Dissipation in Icy Satellites (ODIS), to verify our analytical results.ODIS solves the flux form of the shallow water equations discretized over the sphere using an icosahedral geodesic grid (e.g., Heikes & Randall, 1995).The model has been modified to include the nonlinear advective terms that are required to produce time-averaged flow.These additions follow the implementation described in Thuburn et al. (2009) and Ringler et al. (2010) for the Coriolis, vorticity, and momentum-advection operators.Advection of ocean thickness is added using the third-order accurate transport scheme outlined in Skamarock and Gassmann (2011).These additions to the code have been verified by comparison with test cases 1, 2, and 5 of the classic shallow water test suite outlined in Williamson et al. (1992), as well as the Gaussian hills test proposed by Nair and Lauritzen (2010).
In order to find the mean flow predicted by the theory presented above, we first run ODIS to convergence with the tidal forcing of interest.We define convergence as when the global-and period-averaged kinetic energy varies by <0.1% between successive orbits.Once converged, we take a high-resolution time series of the velocity and height fields and average these over several forcing periods.The averaging removes the fluctuating part of the solutions, and thus yields the steady velocity and height fields that can be compared to both our full and approximate analytical solutions.

Results
The theory outlined above enables us to consider any tidal forcing (or combination thereof) and predict its consequent mean flow.It is general in its capacity to model any weak flow that is generated by periodic or timeinvariant eddy transport in a global ocean.In the present manuscript, our focus is the most prominent mean flow that can be generated by a combination of eccentricity and obliquity forcing.In particular, we consider the combination relevant for the Jovian moon, Europa.Model parameters relevant to Europa are given in Table 1.The tidal response of the ocean is particularly sensitive to its thickness.Because this thickness is only loosely constrained (e.g., Kivelson et al., 2000;Zimmer et al., 2000), we investigate two end-members; a 100 km thick ocean and a near-resonant 220 m thick ocean.46, corresponding drag coefficient, and ice shell drift rate Equation 47 due to the mean flow induced by the westward obliquity tide are also given.The torques are starred because they are upper limits, taken from the maximum of Equation 45.Drift rates Equation 47 are calculated using an ice viscosity of 10 16 Pa s and shear modulus of 3.5 GPa.Obliquities are theoretical values based on the Cassini state calculations from Chen et al. ( 2014), except for Titan, which was constrained by Stiles et al. (2008) and Meriggiola et al. (2016) using radar imagery from the Cassini spacecraft.Uniform ice shell thicknesses are estimates following Hay and Matsuyama (2019), except for Enceladus and Dione, which are based on the global averages from Beuthe et al. (2016) and Hemingway and Mittal (2019).We include Mimas, Titania, and Oberon as potential ocean worlds due to the recent work by Rhoden et al. (2017), Rhoden and Walker (2022), Rhoden (2023), Bierson and Nimmo (2022)

Tidal Mean Flow
We first focus on circulation that arises from forcing due to eccentricity, then obliquity, and finally the case of simultaneous eccentricity-and obliquity-forcing.(d and e) are the m = 2 and 4 components, respectively, and (f) is the sum of all components.There are no odd harmonic components.For m = 0 currents, east, west, and zonal eccentricity-forcing need only be considered separately, whereas all forcing components are required simultaneously to find the m = 2 and 4 components of the mean flow.The solutions are truncated at degree N = 8.The maximum velocities are roughly (d) 4 × 10 5 mm s 1 , (e) 3 × 10 5 mm s 1 , and (f) 8 × 10 5 mm s 1 .that the time-averaged flow contains multiple zonal "jets" (panel a), mostly due to the m = 0 component of the forcing potential.There are two eastward jets at high latitude, and two westward jets at low latitudes.These are the strongest features in this particular scenario, as is evident from the total time-averaged flow shown in panel f.The m = 2 circulation (panel d) contains high-latitude gyres, which are also evident in the total flow.As outlined in Section 2.4, the different zonal wavenumbers of each time-averaged flow component arise due to coupling between various modes in the fluctuating tide.The m = 0 component of the time-averaged flow is a result of a single fluctuating mode interacting with itself, which is why we can split up the flow in the first row of Figure 1 into contributions from westward, eastward, and zonal forcing.This is not possible for m > 0 component of the timeaveraged flow, which are a result of coupling between different modes of the fluctuating tide.In Figure 1, panel d is a result of coupling between the zonal and westward fluctuating eccentricity tide, whereas the m = 4 flow in panel e is due to east-west eccentricity tide coupling.

Eccentricity Forcing
Figure 1c shows that the ocean surface has a time-mean deformation, which we refer to as the time-mean tide.Like the mean flow, this is dominated by the m = 0 component, which leads to a (small) increase in the planetary oblateness.The meridional velocity component v due to eastward eccentricity-forcing, shown in panel b, indicates that the time-mean flow is directed down gradient of the time-mean tide.At first glance, this seems to imply mass imbalance, with a net advection of ocean thickness toward the poles.However, in this case (and the rest of the results in this paper) the time-mean meridional velocity is actually balanced by an eddy transport of ocean thickness, η′v′.This eddy transport is the underlying physical mechanism generating the time-mean tide.That is, η′v′ deforms the ocean surface by advecting ocean thickness, and h 0 v opposes this to conserve mass.This can be understood further by noting that the eddy transport divided by the ocean thickness, η′v′ / h, has units of velocity and is sometimes referred to as the bolus velocity.For linear waves, the bolus velocity is equivalent to the Stokes velocity (Lee et al., 1997).Manipulating Equation 9 for linearized, zonally invariant time-mean flows shows that v ≈ η′v′ / h 0 .The time-mean meridional velocity thus exactly opposes the meridional Stokes velocity.
Finally, in the Europan ocean case considered here, it is important to note that the predicted time-averaged flow speeds are very small.This is because (a) the tidal-forcing is relatively small, and (b) the fluctuating velocity u′ is weak due to the ocean being so thick (thick oceans, when forced via orbital eccentricity, respond predominantly via rotational-gravity wave propagation (Hay & Matsuyama, 2019) which results in slower current speeds than in a thin layer).
Figure 2 shows that stronger mean flows can be generated by a more energetic fluctuating tide.It uses a much thinner, nearly resonant, 220 m thick ocean.Now the zonal jets and m = 0 time-averaged tide are much larger: on the order of centimeters per second and tens of centimeters, respectively.In contrast to the thick-ocean case in Figure 1, the dominant jets are exclusively westward and driven by the eastward-propagating forcing potential.This is because the eastward-propagating fluctuating tide is nearly resonant at this ocean thickness, and thus has the most energetic currents available to drive the mean flow.Another difference is the shape of the dominant timemean tide shown in panel (f), which has a much stronger n = 4 and m = 0 component compared to the dominant n = 2, m = 0 shape in the thick-ocean case in Figure 1f.

Obliquity Forcing
We now consider time-averaged flows driven by obliquity forcing, again assuming a 100 km thick ocean.Europa's obliquity has not been measured, so we use the Cassini state value of 0.053°from Chen et al. (2014).Obliquity forcing is thus rather small.Despite this, obliquity forcing drives the most energetic of all the thickocean time-averaged currents explored here.
Figure 3 shows that there is only one, entirely westward, zonal jet (panel a), which has a maximum westward velocity of 0.05 mm s 1 at the equator.This jet is driven by the Rossby-Haurwitz waves caused by the westward propagating portion of the obliquity forcing (Tyler, 2008); we discuss this in much more detail in Section 4.
The accompanying time-mean tide, shown in Figure 3c, has a small, n = 2 and m = 0, anti-oblate figure.This is an interesting result because even though the ocean is forced by a purely order m = 1 potential (Equation 35), the ocean response includes an m = 0 surface deformation, a consequence of the non-linear effects included in the theory.Obliquity forcing can thus lessen the planetary oblateness driven by rotational deformation.This timemean tide is again caused by the meridional eddy transport of tidal height, η′v′, which in this case forces mass to flow toward the poles.The time-averaged meridional flow in Figure 3b opposes this eddy transport, maintaining steady state.The m = 0 time-averaged current forced by the eastward obliquity potential is completely negligible when compared to that driven by the westward potential.The m = 2 circulation is caused by coupling between the eastward and westward fluctuating tides, which results in a series of gyre-like flows.The total obliquity-driven, time-averaged flow shown in panel f is completed dominated by the westward jet from panel a.  (d and e) are the m = 2 and 4 components, respectively, and (f) is the sum of all components.There are no odd harmonic components.For m = 0 currents, east, west, and zonal eccentricity-forcing need only be considered separately, whereas all forcing components are required simultaneously to find the m = 2 and 4 components of the mean flow.The solutions are truncated at degree N = 16.The maximum velocities are roughly (d) 1 cm s 1 , (e) 2 cm s 1 , and (f) 10 cm s 1 .
Due to the strength of the m = 0 time-averaged circulation in Figure 3a, much of the later discussion in this manuscript is devoted to this component of the time-averaged tidal response.This includes the derivation of analytical solutions for the mean tide and currents, which are plotted as the dashed lines in Figures 3a-3c.

Combined Eccentricity and Obliquity Forcing
The last set of time-averaged flows that we find are shown in Figure 4 for a 100 km thick ocean, which arises only when there is simultaneous forcing from orbital eccentricity and obliquity.This simultaneous forcing leads to odd zonal-wavenumber mean flows, with m = 1 and m = 3 components shown in panels a and b, respectively.
The m = 1 mean flow is especially surprising, and comprises what is effectively a single gyre flowing north-south across the equator and poles.In this respect, it is somewhat similar to the Rossby-Haurwitz wave generated by the obliquity tide (Tyler, 2008), except that in this case it is a zero-frequency wave.This m = 1 component of the mean flow is a result of coupling between the westward obliquity and zonal/westward eccentricity fluctuating tides.
The m = 3 component in panel b is slightly weaker than the m = 1 mean flow, and is characterized by six gyres that straddle the equator and do not extend to the poles.This mean flow component is almost entirely due to interactions between the westward obliquity and eastward eccentricity fluctuating tides.The sum of all the odd-m flow components is shown in panel c of Figure 4, which comprises a rather complex flow pattern featuring various time-mean tidal highs/lows in the low to high-latitudes.
Given that all forcing components are included in Figure 4, it is useful to remind ourselves what the total, timeaveraged flow looks like; that is shown in panel d.It is dominated by the m = 0 westward jet from Figure 3a, which, for our 100 km thick ocean, is much stronger than all other time-mean flow components.There are, however, small deflections of the tidal amplitude contours across lines of latitude, hinting at the contribution from the smaller m = 1 and 3 mean flows.These introduce some zonal shear.

Approximate Analytical Solutions to the Westward Obliquity Tide Mean Flow
We have shown that the strongest mean flow that develops in thick oceans is the westward jet shown in Figure 3a.This jet is driven by the eddy forcing that results from non-linear interactions of the fluctuating tidal flow response to westward obliquity forcing (discussed in Section 4.1).Here, we derive analytical solutions that accurately predict this time-mean circulation.
When the ocean layer is thick, the time-averaged flow that results from westward obliquity forcing is dominated by only three modes: the time-mean tide η 0 2 , velocity potential Φ 0 2 , and stream function Ψ 0 1 .Manipulating the mean-flow system in Equation 16 with only these three modes yields, To solve this system we must find the four eddy-forcing coefficients, H 0,s 2 , Z 0,s 1 , Z 0,t 2 , and J 0 2 , which are a result of nonlinear interactions of the fluctuating tide.
As shown by Chen et al. (2014), the fluctuating ocean tide due to westward obliquity-forcing can be approximated by two westward-traveling modes when the ocean is thick: the stream function Ψ 1,W 1 and velocity potential Φ 1,W 2 .The stream function is the biggest of the fluctuating modes and can be approximated in the thick ocean limit as where ϵ = 4Ω 2 R 2 /gh 0 is Lamb's parameter.This equation is equivalent to Equation 58in Beuthe (2016) and Equation 16in Hay and Matsuyama (2019) if the ice shell is not taken into account.The factor outside the brackets is the inviscid solution first determined by Tyler (2008).Chen et al. (2014, Equation 32) showed that Ψ 1,W 1 and Φ 1,W 2 are straightforwardly related to each other through drag (Equation B7).The eddy-forcing and its spherical harmonic coefficients H 0,s 2 , Z 0,s 1 , Z 0,t 2 , and J 0 2 , can therefore be written analytically using only Ψ 1,W 1 , as shown in Equation B2.Substituting these analytical expressions for the eddy coefficients into Equation 37 eventually yields a solution for η 0 2 , Φ 0 2 and Ψ 0 1 .Then, plugging the latter two modes into Equation 14and then Equation 10 yields an analytical solution for the time-mean velocity field.Together, The analytical expressions show that all of the mean flow components depend on the obliquity angle squared; in the inviscid limit, only the north-south velocity tends to zero.The above solutions are shown as the dashed lines in Figures 3a-3c and D1, where we see that the match with the full semi-analytical solution is excellent.A detailed derivation of the solutions in Equation 39 is given in Appendix C.

Sensitivity to Drag
The analytical solutions derived above are useful to illustrate the role of drag in damping the magnitude of the westward obliquity-forced mean flow.The actual drag coefficient of these ocean worlds, as well as the coefficient's dependence on spatial and temporal scales, is unknown.The exception to this is quadratic drag, where c D is constrained to be on the order of 10 3 -10 2 on Earth (e.g., Egbert & Ray, 2001); there is little reason to suspect that it would be significantly different at the boundaries of an internal water ocean.Our analysis above has focused on linear drag mechanisms as these are more analytically tractable and the results so far have only used a Rayleigh drag term αu.Below, we further consider the sensitivity of our results to α, as well as to the Navier-Stokes-type viscosity, ν.
It is helpful to recast the effective drag, regardless of the mechanism, in terms of a tidal dissipation factor Q. We follow Beuthe (2016) and define the dissipation factor at spherical harmonic degree n as Here, Q n depends on spherical harmonic degree because Navier-Stokes drag depends on the spatial scale of the flow and is more efficient at smaller wavelengths (larger n).As defined here, Q n does not necessarily inform us about the efficiency of tidal heating, which is how it is classically used (Goldreich & Soter, 1966), because it does not take into account the decrease of the maximum oceanic kinetic energy as drag becomes increasingly large.This point is discussed in detail in Hay and Matsuyama (2017) and Appendix D of Matsuyama et al. (2018).Q n does, however, provide a convenient way to contrast the effect of either α or ν at a fixed n.We therefore find it helpful to view Q n , at least in the following, as the inverse of an effective, non-dimensional, drag coefficient.In panels a and c we see that the obliquity-forced mean flow is not greatly altered by our choice of drag mechanism.The solution asymptotes to a maximum jet speed in the limit of large Q 1 .As Q 1 is decreased, the maximum jet speed also decreases.This behavior can be anticipated from our analytical solution in Equation 39a.
In contrast to the obliquity-forced mean flow, the drag mechanism has a profound effect on the structure of the zonal flow in the near-resonant case in Figures 5b and 5d.Linear drag (b) is far more damping than viscous drag (d), causing the mean flow to be roughly three times faster in the viscous case, despite using similar values of Q.The spatial form of the jets is also different.Linear drag creates two distinct mid-latitude jets, while viscous drag permits, roughly, one single jet.This difference is due to the fact that linear drag preferentially damps the fastest parts of the flow, whereas viscous drag acts to smooth the sharpest gradients in the flow.These findings are in accordance with Huthnance (1981), who showed that the Lagrangian mass transport associated with long waves in the ocean would be controlled by the drag mechanism, but be independent of the drag magnitude, provided that drag was small enough.We conclude from Figure 5 that eccentricity-forced mean flow strongly depends on the type of drag mechanism, whereas westward obliquity-forced mean flow does not.The obliquity-forced mean flow is unaffected by the drag mechanism because the αu and (ν⁄ R 2 )∂ θθ ū terms in the time-mean zonal velocity balance (Equation 42) have the same spatial form for the east-west velocity solution in Equation 39a-a mathematical coincidence.In all cases, when the magnitude of drag is weak, the time-mean zonal jets approach a maximum (but finite) speed.

Discussion
We have shown that our weakly nonlinear theory predicts several types of time-averaged currents driven by periodic ocean tides, the fastest of which is a single westward jet (Figure 3) that originates from the westwardpropagating obliquity-forcing potential (Equation 35).In Section 4.1 we examine the physical mechanism that drives this jet, and in Section 4.2 we consider the effect that this flow may have on the rotation of floating ice shells.Finally, in Section 4.3 we re-summarize the assumptions that have lead to our results, and caveats and limitations that arise from them.

Driving Mechanism for Westward Zonal Jets
The analytical solutions presented in Section 3.2 do not provide accessible physical insight into the mechanism through which westward obliquity tides drive retrograde mean flow.For that, we turn to the forcing periodaveraged vorticity balance.Neglecting second-order terms in time-mean quantities in Equation 7, the periodaveraged vorticity balance becomes, Here, vorticity accumulates/depletes where there is any imbalance in the fluxes of planetary ( f ) and relative vorticity (ζ′), except where it is dissipated by drag.Because zonal jets are zonally invariant, we can take a zonal average of the above which eliminates the derivatives of all quantities with respect to ϕ. Multiplying the result by R sin θ and integrating with respect to latitude gives a momentum balance for the zonal-mean, east-west velocity, This equation is equivalent to Equation 15.21 in Vallis (2017), except that our system can be horizontally divergent so v is non-zero.Term I is the meridional transport of eddy vorticity by the fluctuating flow, II is the meridional transport of planetary vorticity by the mean flow, and term III is damping of the time-mean zonal velocity field by drag.Thus, Equation 42 tells us that, initially, the east-west flows develop through an imbalance in the meridional transports of vorticity.At steady state, this imbalance is countered by damping due to drag.
Term I is non-zero due to the tilting of the Rossby waves generated by the obliquity tidal forcing.The tilting is illustrated in Figure 6a, which shows the relative vorticity (ζ′, black contours) and meridional velocity (v′, solid contours) associated with these waves using α = 10 6 s 1 .Drag tends to distort the Rossby waves into a "bow" shape, so that lines of constant vorticity tilt toward the north-east and south-east in the northern and southern  40).The values used are Q 1 = 1, 5, 10, 50, 100, 500 and 1,000, which spans α ≈ 10 5 s 1 to 10 8 s 1 , and ν ≈ 1.25 × 10 7 m 2 s 1 to 1.25 × 10 4 m 2 s 1 .
hemispheres, respectively (this shape can also be seen in obliquity-forced ocean tidal heating maps (Matsuyama et al., 2018, Figure 8b)).Figure 6b shows term I by taking a zonal-mean of the product of ζ′ and v′ from the snapshot in panel a.Panel b thus shows that correlation between ζ′ and v′ due to the tilting of Rossby waves causes a northward transport of eddy vorticity at high latitudes, while southwards transport occurs near the equator.Figure 6b also shows term II, vf , which is slightly larger than term I and transports planetary vorticity southwards, except near the equator.The sum of terms I and II is therefore not in balance, meaning that overall vorticity is transported from the northern to southern hemisphere.This net transport of vorticity causes an equatorially symmetric east-west deceleration of the ocean and thus the generation of the westward jet seen in Figure 3a.
Each of the three labeled terms in Equation 42 inherently depends on drag.This is obvious for term III, but less so for terms I and II.The Rossby waves driving the mean flow only develop a tilt (and therefore correlation of v′ and ζ′) when drag is present, meaning that term I tends to zero in the inviscid limit.The mean meridional velocity, v, arises to balance the eddy transport of tidal height, η′v′.Because η′ tends toward zero in the absence of drag (Tyler, 2008), this means v and therefore term II also tends to zero in the inviscid limit.Thus, in the (physically impossible) limit of no drag, terms I, II, and III in Equation 42 all vanish identically.Nevertheless, the jet speed remains finite and at its maximum in this limit, as predicted by our analytical solution in Equation 39a and shown in Figure 5.The reason for this is that terms I and II vary linearly with α and ν, provided that those parameters are small.At steady state, the dependence on α or ν cancels out in Equation 42and the solution for u becomes independent of the magnitude of drag/viscosity (but not necessarily the drag mechanism, as shown in Figure 5).

Ice-Ocean Torque
Persistent circulation can exert a torque on the ice shell because momentum is exchanged between the ocean and ice due to viscous drag in the mechanical boundary layer.These torques may play a role in ice shell reorientation over geological time, which has recently been proposed for buoyancy-driven zonal jets (Ashkenazy et al., 2023;Hay et al., 2023;Kang, 2024).This effect is relevant for the m = 0 and m = 1 jets that emerge from westward obliquity-forcing and coupled obliquity-eccentricity forcing, respectively (Figures 3a and 4a).Below, we consider the torque due to the strongest m = 0 mean flow.
We take the simplest possible method to estimate the axial torque that is aligned with the rotation axis T z .The torque is given by integrating the boundary layer stresses multiplied by the distance to the rotation axis, over the surface area of the ice-ocean interface, which we assume to be uniform.Hence, The ocean thickness is held at h 0 = 10 km, the drag coefficient is α = 10 6 s 1 , and the solution is truncated at N = 3.
Considering only linear drag (ν = 0, b 1 = α/2 Ω), and assuming that all friction occurs only within the boundary layers, the zonal stress field may be approximated as, The factor of 1/2 appears because, in our depth-averaged model, stresses are equally partitioned between the iceocean and core-ocean interfaces.Substituting our zonal-mean zonal velocity from Equation 39 into the above and evaluating the integral in Equation 43gives This function peaks when α = Ωϵ/40 = Ω 3 R 2 /10gh 0 .Plugging this value of α into the above gives us the maximum possible retrograde torque that the westward obliquity mean flow can exert through linear drag, where ρ is the bulk density of the satellite and G is the universal gravitational constant.The maximum torque is thus independent of ocean thickness, and is proportional to the obliquity angle squared.The peak in the retrograde torque is similar to the peak in tidal heating expected from obliquity tide oceanic dissipation at particular values of drag coefficient (e.g., Downey et al., 2020;Hay & Matsuyama, 2017;Tyler, 2011).When α is small, drag has a negligible effect on the fluctuating tide and the magnitude of the mean flow, so the surface stresses can grow as drag is increased.Once α is large enough, however, drag damps the fluctuating tide which lowers the ocean's kinetic energy and therefore the mean flow and surface stresses.The maximum retrograde torque occurs between these two states.
Table 1 shows the torque calculated from Equation 46for a set of ocean worlds in the solar system.The obliquity values in Table 1 are from Chen et al. (2014), where it is assumed that the rotation of each satellite lies in a low Cassini state.We find that because of the R 5 dependence of the maximum torque in Equation 46, the larger satellites can generally experience more significant torques.Satellites that are further away have slower rotation rates due to being synchronously locked, which by itself decreases the torque experienced.This effect is countered by an increase in the obliquity angle (i.e., forcing amplitude), which is generally larger for bigger, further out satellites as they can have more inclined orbits.Larger inclinations result in more significant obliquities, according to the Cassini state relationship (Ward, 1975).
The maximum retrograde torques in Table 1 are ∼10 15 N m, which occurs for Europa, Titan, and Triton.These torques are surprisingly large given how weak the jet speeds are in Figure 3a and Table 1, but the fact that there is only one westward jet means that all surface stresses contribute in the same sense to the torque (i.e., there is no cancellation due to opposing jets).Our calculated torques can approach, though are typically less than, the tidal torque that the central body exerts on the viscoelastic ice shell.
Also shown in Table 1 is the drift rate of the surface, calculated using the viscoelastic ice drift model from Ashkenazy et al. (2023).This model can be rewritten as where the ice shell's moment of inertia I = 8πρ s R 5 (R h s ) 5 )/ 15, ρ s = 10 3 kg m 3 is the density of the shell, k i is the elastic parameter defined in Equation 10 of Ashkenazy et al. (2023), and τ M is the Maxwell time of the ice shell, given by the ratio of the shell's viscosity to shear modulus.Following Kang (2024), we neglect the relative angular speed of the ice shell as it is one to two orders of magnitude smaller than the retrograde jet speed calculated above.For a relatively high ice viscosity of 10 16 Pa s and typical shear modulus of 3.5 GPa, retrograde drift rates of a few tens of m/yr may occur on Europa, Callisto, and Triton.Titan's surface could potentially experience a drift rate magnitude of over 300 m/yr, though this will be complicated by angular momentum exchange between Titan's surface and atmosphere.The obliquity-tide driven drift for Europa may be detectable by the Europa Clipper spacecraft if the ice shell has an average viscosity lower than the 10 16 Pa s ice viscosity assumed in Table 1.Callisto's ice shell drift may also be detectable using imagery from the JUICE spacecraft.Titan's maximum possible drift rate corresponds to an angular velocity of ∼ 0.0079°/yr, which falls within the 1 σ non-synchronous rotation rate of 0.024 ± 0.018°/yr measured from Cassini radar imagery in the period 2004-2009 (Coyette et al., 2018;Meriggiola et al., 2016).

Caveats and Limitations
The theory presented in this manuscript is based on the following assumptions.First, we assume that the tidal flow of an ocean is well described by the LTE.Second, we assume that the ocean is homogeneous and uniform in depth.Third, we do not account for the mechanical suppression of and gravitational effect of the ice shell on the ocean (Beuthe, 2016;Matsuyama et al., 2018), so we assume the effect of the ice shell can be neglected.Fourth, the mean flow is assumed to be much weaker than its fluctuating counterpart, and both are weak enough that their governing equations can be linearized.Fifth, we assume that the underlying physical drag mechanism that damps the flow field is identical for both the fluctuating tide and the mean flow.Finally, we neglect other types of oceanic currents in our theory.Below we discuss the validity of each of these assumptions.
The LTE are roughly valid provided that the thickness of the ocean layer is less than 5% to 10% of the body's radius.This is an acceptable assumption for many of the ocean worlds considered here, with the notable exception of Enceladus (Aygün & Čadek, 2023) and probably Dione.The shallow water assumption precludes the emergence of internal, inertial tidal waves (Rekier et al., 2019;Rovira-Navarro et al., 2019), which would be another manifestation of time-averaged currents, though this is unlikely to affect the results here.A more significant limitation is the second assumption that the background ocean thickness h 0 is spatially uniform.Rovira-Navarro et al. (2020) and Rovira-Navarro et al. (2023) demonstrated that small meridional variations in ocean thickness can reduce the fluctuating oceanic velocities induced by the westward obliquity tide.Because the largest mean flows predicted here depend on the magnitude of the fluctuating tidal flow-speed squared Equation C3, meridional ocean thickness variations may reduce our calculated maximum zonal flow speeds.Future work can take this effect into account.
The mechanical effect of the ice shell decreases the response time of the ocean to tidal forcing, while the additional self-gravitation due to ice deformation enhances the effective tidal force the ocean experiences.The former tends to decrease the flow speeds associated with gravity waves, while the latter increases the flow speeds of all waves (Hay & Matsuyama, 2019).Self-gravity is the dominant effect for obliquity tides, which means that our theory will slightly underestimate the speed of the persistent westward jet driven by the westward obliquity forcing.The mean flow speeds due to all other tidal forcing components will be slightly overestimated because the mechanical suppression of the ocean tide is the more prominent effect.For most of the ocean worlds considered in this work, mechanical suppression of the tide is negligible, as shown in Hay and Matsuyama (2019).Again, Enceladus and Dione are the exceptions to this, where the mechanical effect of the ice-ocean coupling can be significant.
Our fourth assumption of weak fluctuating tides and even weaker mean flows is not violated in any of the results presented in this manuscript, though it may be possible to find instances where the mean flow exceeds the fluctuating tidal currents (e.g., Moore, 1970).For larger tidal forcing or more resonant ocean thicknesses than considered here, it may be that nonlinear terms in the momentum equations cannot be neglected.This would preclude the linearization of both the fluctuating and mean flow LTE systems.Partial decoupling of the mean and fluctuating flow will still be possible as long as the mean flow is much weaker than the fluctuating currents.
We assume that the drag coefficient and mechanism acting on the fluctuating tide and mean flow are identical.In principle, this does not have to be the case.For example, perhaps the barotropic tide is damped primarily through conversion to and breaking of internal tidal waves, while the mean flow is damped through nonlinear top/bottom drag.Without any constraint on the dominant physical drag mechanism affecting these two flows, it is not appropriate to further decompose them, though it may be interesting to explore this further.
The eddy transports are key to determining the mean flow because they act as the forcing in system Equation 16.
In this work we consider only the long-wavelength tidal contributions to the eddy transports, but there are certainly many others.Convection and libration, for example, may alter the eddy transports of momentum that we calculate and, consequently, the mean flow.Our theory has the advantage that the mean flow system is decoupled in zonal wavenumber m.This implies, for example, that only when other mechanisms produce an eddy transport at m = 0 will they alter the m = 0 mean flow.Convection in rotating systems yields zonally symmetric momentum fluxes (Busse, 2002;Kaspi, 2008), as does magnetic pumping (Gissinger & Petitdemange, 2019), so it is an intriguing possibility that our m = 0 results may change if these are taken into account.It is not clear, however, if any mechanism other than tides can produce time-mean eddy transports at zonal wavenumbers m ≥ 1, unless there are inherent lateral inhomogeneities.A possible instance of this is seafloor heatflux variations due to solid-body tidal heating (e.g., Beuthe, 2013;Lemasquerier et al., 2023).

Conclusion
In this manuscript we have developed a theory for the generation of time-average currents in global liquid layers driven by periodic forcing.The theory is applied to the subsurface ocean moons of the outer solar system, with a focus on Jupiter's moon Europa.We force these oceans with the tide-raising gravitational potential due to orbital eccentricity and the body's rotational obliquity to find the residual circulation in both thick (∼100 km) and thin (∼0.1 km) oceans.We investigate the theory in three ways: a fully numerical method, a spectral semi-analytical method, and an approximate analytical method.All three produce consistent results.
We find that periodic tidal forcing can generate a variety of time-averaged currents, with zonal wavenumbers of 0 through 4. The dominant time-averaged current for a thick ocean on Europa is a single, retrograde jet that manifests when the ocean is forced by the obliquity tidal potential.Using analytical solutions for the retrograde jet, we find that the jet speed generally falls between 0.01 to 1 mm s 1 for a range of icy ocean worlds, and can impart a significant torque on the ice shell.Using the ice-drift model of Ashkenazy et al. (2023), we find that the maximum possible torques can induce retrograde ice drift rates of 10s m yr 1 for Europa, Callisto, and Triton, and ∼300 m yr 1 in the case of Titan, though these values depend on the assumed ice shell viscosity.We find that thinner oceans or stably stratified layers can produce much more energetic time-average currents.There is therefore potential for stably stratified portions of icy satellite oceans (e.g., Ashkenazy & Tziperman, 2021;Kang, 2023) to experience non-trivial mean currents induced by tidal forcing, which should be a focus of future work.
This study highlights and overcomes some of the limitations of linear tidal theory.Doing so predicts the emergence of mean flows generated by tides, thus highlighting a path to induce large-scale transport within planetary interiors.We focus on internal oceans, but the general approach here could be applied to, for example, time-mean circulation within partially molten systems such as the core of Enceladus, Io's asthenosphere, or even convective ice layers.Given the importance of cycling chemistry across planetary interiors, tidally induced mean flows may play a role in helping maintain habitable subsurface environments.
After substituting the spherical harmonic and Fourier expansions from Equations 23 to 27 into the above, along with some substantial algebra, we get, H s nm = N m n 4n(n + 1) ) To ease the notation, we have introduced the nonlinear coupling coefficients, ), (B4a) ), (B4e) which can all be calculated entirely by solving the linearized LTE algebraic system (Equation 24).The A, C, and E coefficients capture coupling between modes traveling in the same direction, while B, D, and F account for coupling between modes traveling in opposite directions.The remaining unknowns in Equation B3 that are needed in order to find the eddy coefficients are the integrals of the form, While analytical solutions for the integral of the product of three scalar spherical harmonics exist (Arfken & Weber, 1999, Equation 12.189), we are unaware of similar solutions for the above integrals that involve a combination of scalar and vector harmonics.Fortunately, these can all be calculated numerically for each relevant combination of n and m.Like the triple spherical harmonic integral, many of these combinations are trivially zero (Arfken & Weber, 1999, p. 803).

B2. Analytical Eddy Forcing Coefficients for the Westward Obliquity Tide
A power spectrum of the m = 0 solutions to the westward obliquity-forced mean flow shows that degrees 1 and 2 dominate, restricting the majority of the unknowns to Ψ 0 1 , Φ 0 2 , and η 0 2 .From inspection of the mean flow system (Equation 16), this means that the relevant eddy-forcing coefficients are restricted to these degrees as well, Z 0,s 1 , Z 0,t 2 , J 0 2 , and H 0,s 2 .These four coefficients can be well approximated by only the largest two modes of the fluctuating oceanic response to the westward obliquity forcing, Ψ 1,W 1 and Φ 1,W 2 (Chen et al., 2014;Hay & Matsuyama, 2019).The only possible couplings between these modes yields the non-linear coupling coefficients A + 1121 , A 1121 , C + 1111 , and E + 2121 from Equation B4.Working with only these modes, we can approximate the eddyforcing coefficients using Equation B3, where only the largest coefficients are retained.The coupling coefficients themselves can be simplified using the truncated LTE in Equation 24 at degrees 1 and 2 following Chen et al. (2014), which relates the stream function and velocity potential via, Substituting this into the coupling coefficients in Equation B4 and noting that L 1,W 1 = 0, gives, After substituting these into Equation B6 above, we get approximate solutions for the four eddy-forcing coefficients, with the westward propagating obliquity potential.The numerical results correspond to Grid Levels 5, 6, 7, each of which doubles the resolution of the previous.GL7 discretizes the sphere over 40,962 finite volumes, giving an average horizontal resolution of about 30 km on Europa.Here we use a large obliquity of 1°and drag coefficient of α = 10 7 s 1 to make the signal of the mean flow easier to detect in the numerical results.The ocean thickness is 10 km.The features in the numerical solutions at latitudes of 25°are produced by noise associated with the locations of pentagons in the numerical grid.The approximate analytical solution (dashed line) and full analytical solution (bold red line) are also shown.

Figure 1
Figure 1 shows each component of the mean flow due to eccentricity forcing in a 100 km thick ocean with e = 0.01.The first row (panels a-c) show the zonally invariant component; the second row are the m = 2 and m = 4 components (panels d-e); the total flow field is in the last row (panel e).The major point of interest here is

Figure 1 .
Figure 1.Components of the time-mean currents driven by eccentricity-forcing at Europa, using α = 10 8 s 1 , h 0 = 100 km, and e = 0.01.Panels (a-c) are the m = 0 zonal-mean components,(d and e) are the m = 2 and 4 components, respectively, and (f) is the sum of all components.There are no odd harmonic components.For m = 0 currents, east, west, and zonal eccentricity-forcing need only be considered separately, whereas all forcing components are required simultaneously to find the m = 2 and 4 components of the mean flow.The solutions are truncated at degree N = 8.The maximum velocities are roughly (d) 4 × 10 5 mm s 1 , (e) 3 × 10 5 mm s 1 , and (f) 8 × 10 5 mm s 1 .

Figure 2 .
Figure 2. Components of the time-mean currents driven by eccentricity-forcing at Europa, using α = 10 8 s 1 , h 0 = 220 m, and e = 0.01.Panels (a-c) are the m = 0 zonal-mean components,(d and e) are the m = 2 and 4 components, respectively, and (f) is the sum of all components.There are no odd harmonic components.For m = 0 currents, east, west, and zonal eccentricity-forcing need only be considered separately, whereas all forcing components are required simultaneously to find the m = 2 and 4 components of the mean flow.The solutions are truncated at degree N = 16.The maximum velocities are roughly (d) 1 cm s 1 , (e) 2 cm s 1 , and (f) 10 cm s 1 .

Figure 3 .
Figure 3.All components of the time-mean currents driven by obliquity-forcing at Europa, using α = 10 8 s 1 and h 0 = 100 km.Panels (a-c) are the m = 0 zonal-mean components, (d) is the m = 2 component, and (e) is the sum of all components.For m = 0 currents, east and west obliquity-forcing need only be considered separately, whereas both east and west forcing components are required to find the m = 2 flow.The m = 0 analytical solutions derived in Section 3.2 are shown as the dashed lines in the top row.The maximum velocities are (d) 0.0002 mm s 1 and (e) 0.05 mm s 1 .

Figure 4 .
Figure 4.Additional components of the time-mean currents that arise when driven by simultaneous eccentricity-and obliquity-forcing at Europa, using α = 10 8 s 1 and h 0 = 100 km.Panels (a and b) are for the only odd harmonics, m = 1 and 3, respectively, (c) is for the sum of the odd harmonic components, and (d) is the complete timemean current that includes all components.The maximum velocities are (a) 0.002 mm s 1 , (b) 0.002 mm s 1 , (c) 0.003 mm s 1 , and (d) 0.05 mm s 1 .

Figure 5
Figure5shows how the drag magnitude and mechanism alters the dominant mean flows.The top row uses linear drag (a, b); the bottom row uses viscous drag (c, d).The left column (a, c) is for westward obliquity-forcing and a thick ocean, while the right column (b, d) is for eastward eccentricity-forcing and a thin ocean.Each curve represents a different Q 1 , where we have chosen n = 1 because the biggest mode in the obliquity-driven mean flow is Ψ 0 1 .In panels a and c we see that the obliquity-forced mean flow is not greatly altered by our choice of drag mechanism.The solution asymptotes to a maximum jet speed in the limit of large Q 1 .As Q 1 is decreased, the maximum jet speed also decreases.This behavior can be anticipated from our analytical solution in Equation 39a.

Figure 5 .
Figure 5.Effect of drag on the time-averaged m = 0 zonal currents.Panels (a and b) use a linear drag term, while (c and d) use a Navier-Stokes viscosity.The left column uses westward obliquity forcing and a thick ocean, while the right column uses eastward eccentricity forcing and a thin ocean.Each colored curve represents a different value of Q 1 (Equation40).The values used are Q 1 = 1, 5, 10, 50, 100, 500 and 1,000, which spans α ≈ 10 5 s 1 to 10 8 s 1 , and ν ≈ 1.25 × 10 7 m 2 s 1 to 1.25 × 10 4 m 2 s 1 .

Figure 6 .
Figure 6.Mechanism for the westward obliquity north-south vorticity transport that drives the retrograde jet.(a) Fluctuating vorticity (black) and meridional velocity (colored) contours at a snapshot of the ocean's response to westward obliquity forcing.(b) Corresponding zonal-mean meridional transport of vorticity, calculated by taking the zonal mean of the product of v′ and ζ′ from panel (a).Positive and negative values indicate southward and northward transport of vorticity, respectively.The ocean thickness is held at h 0 = 10 km, the drag coefficient is α = 10 6 s 1 , and the solution is truncated at N = 3.

Figure D1 .
Figure D1.Numerical benchmarking.Ocean Dissipation in Icy Satellites produces a westward jet when the ocean is forced with the westward propagating obliquity potential.The numerical results correspond to Grid Levels 5, 6, 7, each of which doubles the resolution of the previous.GL7 discretizes the sphere over 40,962 finite volumes, giving an average horizontal resolution of about 30 km on Europa.Here we use a large obliquity of 1°and drag coefficient of α = 10 7 s 1 to make the signal of the mean flow easier to detect in the numerical results.The ocean thickness is 10 km.The features in the numerical solutions at latitudes of 25°are produced by noise associated with the locations of pentagons in the numerical grid.The approximate analytical solution (dashed line) and full analytical solution (bold red line) are also shown.

Table 1
Satellite Parameters and Westward Obliquity Mean Flow Torque Results for Select Ocean Worlds Note.Model parameters are the surface radius R, bulk density ρ b , angular rotation rate Ω, obliquity, θ o , and ice shell thickness, h s .The maximum retrograde zonal velocity from Equation 39, (9/4)ΩRθ 2 o , maximum torque from Equation