Elsevier

Computers & Geosciences

Volume 155, October 2021, 104865
Computers & Geosciences

Solving the Laplace Tidal Equations using Freely Available, Easily Extensible Finite Element Software

https://doi.org/10.1016/j.cageo.2021.104865Get rights and content

Highlights

  • We study the use of General-purpose Finite Element Program (PDE2D) to solve Laplace Tidal Equations.

  • We calculate ocean dissipation in three icy satellites, Enceladus, Europa and Titan.

  • We show that modeling results using PDE2D can be a useful tool for researchers studying tidal heating related problems.

Abstract

The Laplace Tidal Equations (LTEs) play a basic role for investigating the dissipation of long-term orbital energy as heat in surface oceans of icy satellites. Here we address solving LTEs by means of finite elements using PDE2D, a general-purpose finite element program for solving multidimensional PDEs. We use PDE2D to solve the equations for eccentricity and obliquity gravitational potentials for a shallow water approximation, and calculate ocean dissipation in three icy satellites, Enceladus, Europa and Titan. Compared with previous modeling procedures, PDE2D solves for unknowns that are approximated by piecewise cubic polynomials, thus attaining a higher O(4) + O(4) numerical precision. Our modeling results show, within few percent errors, a good agreement with earlier semi-analytical models and numerical solutions. We provide an open access to our PDE2D codes which are intended to be simple and easy to adjust for a wide range of tidal heating related problems. Our modeling results show that PDE2D can be a useful tool for researchers studying tidal heating related problems.

Introduction

Thermal energy within the interiors of rocky planets, including Earth, and outer Solar System dwarf planets and satellites is maintained over periods of billions of years primarily by radiogenic decay of naturally radioactive isotopes with significant contributions from the decay series of 238U, 235U and 232Th, and from the isotope 40K (Birch, 1954). The radioactive isotope content of the icy moon interiors is not well known. However, based on chondritic meteorite samples and compared with the Earth's crust, their concentration into the silicate mantle is quite low (Heier, 1965). Whereas the radiogenic decay is important for bodies with a massive mantle, small satellites dissipate efficiently the radiogenic heat over time scales less than the age of our Solar System (Roberts and Nimmo, 2008). Therefore, small celestial bodies, with the exception of Europa (Nimmo and Pappalardo, 2016), should be presently geologically inert. Nevertheless, past and recent space missions (i.e. Voyager, Galileo and Cassini) revealed a totally different picture. For example, Io, a moon of Jupiter, is a prime example where volcanic eruptions are visible from space. Unexpectedly, the small Saturnian icy moon Enceladus and the smallest of the four Galilean moons orbiting Jupiter, the icy Europa, have the youngest, smoothest and brightest surfaces in our Solar System (Porco et al., 2006; Smith et al., 1979). This clearly suggests the existence of a different long-term internal heat source that is maintaining high temperature conditions inside these icy satellites and driving geological resurfacing. Some key observations, as the 2:1 eccentricity-type resonance of Enceladus with the moon Dione, or 4:2:1 orbital resonance for the Io, Europa, and Ganymede, yield a forced nonzero eccentricity for both Enceladus (e = 0.0047) and Europa (e = 0.0094) (Peale, 1999; Prockter and Pappalardo, 2007; Meyer and Wisdom, 2008). Therefore, flexing arising from the eccentric orbits periodically deform the icy moons, producing tidal heating (Matsuyama et al., 2018). As long as the nonzero orbital eccentricities are preserved, tidal heating probably represents the main long-term internal heat source. Another remarkable example of a young icy surface satellite is Saturn's largest moon, Titan, whose young surface is controlled by chemical atmospheric weathering (Neish et al., 2015). Compared with Enceladus and Europa, Titan has the largest eccentricity (e = 0.029), but its orbital eccentricity is not forced by other satellites, and the evidence that Titan has a subsurface ocean is based on its k_2 Love number (Iess et al., 2012; Durante et al., 2019). Other lines of observations observed at Enceladus, as spewing, possible at supersonic rates, plumes containing water vapor (Postberg et al., 2011; Nimmo et al., 2014). Additionally, magnetic field deflections also indicate the presence of a subsurface liquid ocean on Europa (Zimmer et al., 2000; Kivelson et al., 2000). Such liquid layer will decouple the icy shell from an inner rocky porous core (McKinnon 2015; Thomas et al., 2016). Although the presence and extent of a subsurface liquid ocean comes from different indirect observations, the best evidence for a global ocean beneath the ice crust on Enceladus is provided by the presence of harmonic oscillations or physical librations (Thomas et al., 2016). These observations corroborated with gravity and topography data (Iess et al., 2014; Beuthe et al., 2016) provide convincing evidence of a global subsurface water ocean on Enceladus.

Evidences for water plumes are also proposed for Europa (Roth et al., 2014; Sparks et al., 2016, 2017), although recent studies suggest that the amount of water vapor on Europa released by plumes is less than previously estimated (Paganini et al., 2019).

Following the observation that most of Earth's tidal dissipation takes place within the large water bodies as oceans (Lambeck, 1977), in the last decades or so, a major effort has been made to better understand the tidal ocean dissipation for icy satellites. Tidal dissipation inside the porous rocky cores (Choblet et al., 2017) represents a major source of heat, but thermal energy released due to the motion of subsurface liquid oceans might also contribute to preventing the subsurface water from freezing (Tyler, 2011, 2014). Beuthe (2016) shows that the icy crust controls the amount of heat dissipated by the subsurface oceans, and for small bodies as Enceladus crustal dissipation becomes significant, although crustal dissipation is usually smaller for small bodies because their crust deforms less, compared to a large body such as Europa. Although the ocean tidal dissipation is quite reduced when compared to solid-body tidal heating (Chen et al., 2014), the possibility of tidal heating due to turbulent dissipation within subsurface oceans specific for icy satellites is still a matter of debate (Ross and Schubert, 1989; Tyler, 2011).

Mathematical solutions to solve the equations of continuity and motion used to describe ocean tidal flow for a shallow water layer, known also as the Laplace Tidal Equations or LTEs (Lamb, 1932), range from analytical (Sagan and Dermott, 1982), semi-analytical models (Tyler, 2011, 2014; Chen et al., 2014; Beuthe, 2016; Matsuyama et al., 2018), to more advanced numerical solutions (Sears, 1995; Hay and Matsuyama, 2017; Matsuyama et al., 2018), including recent studies based on finite-elements and finite-volume methods (Rovira-Navarro et al., 2020; Hay and Matsuyama, 2019).

In this paper we employ finite elements for solving the LTEs using a shallow water approximation for icy moons. More specifically, we use a general-purpose finite element program (PDE2D – www.pde2d.com) for calculating ocean dissipation in three icy satellites, Enceladus, Europa and Titan, considering the shallow water approximation. The ocean thickness in these icy moons is considered to be in the range from a few to tens of kilometers, and on the order of 100 km thick for Europa or even up to 400 km for Titan (Moore et al., 2001; Nimmo et al., 2003; Bills and Nimmo, 2011; Iess et al., 2012; Cadek et al., 2019). The shallow water equations are appropriate when the water depth is just few percent of their radius. However, ocean thicknesses for Europa, Enceladus and Titan are considered about 5–15% of their radius and the shallow water approximation might not be the best choice (Rovira-Navarro et al., 2019; Rekier et al., 2019). However, the overall advantage of the shallow water approximation is that we can employ a system of two-dimensional partial differential equations that is relatively faster to solve compared with three-dimensional models. Applying the shallow water approximation on which the LTEs depend, we limit the water flow through the subsurface ocean to the horizontal direction. Since there are no observations of flows in subsurface oceans, the assumption that ocean flow is restricted to lateral, or horizontal movement, is not well constrained. However, some studies suggest that hydrothermal plumes can generate geostrophically balanced small-scale eddies, that can interact with others forming larger eddies (Vance and Goodman, 2009; Soderlund et al., 2020). These larger eddies can merge over time and finally form east-west ocean water jets (similar to those observed on the surface of giant gas planets).

We benchmark our codes against previous results from semi-analytical models (Chen et al., 2014) as well as for previous numerical solutions using finite differences (Sears, 1995; Hay and Matsuyama, 2017), finite elements (Rovira-Navarro et al., 2020) and finite volumes (Hay and Matsuyama, 2019). Our Fortran PDE2D codes provided in this research paper are reusable and also extensible to a wide range of potential problems related with tidal heating. We make fully available our codes for the scientific community, and papers treating similar topics can benefit from the PDE2D flexible design.

Section snippets

Laplace Tidal Equations

Assuming that vertical ocean flow is negligible when compared to lateral flow, the LTEs corresponding for a two-dimensional problem for a synchronously rotating satellite consist of the conservation of mass and momentum, as given in (Hay and Matsuyama, 2017, equations (3), (4))):ηt+h0u=0ut+2Ω×u+αu+CDh0|u|u+gη=(1+k2h2)Uwhere h0 + η(φ,θ, t) is the ocean surface height, u(φ,θ, t) is the horizontal ocean velocity vector (depth averaged), φ is longitude and θ is latitude. Equation (1)

PDE2D

For the solution of the LTEs, we employ PDE2D (www.pde2d.com), a finite element program which solves very general partial differential equations (PDEs). PDE2D has been developed by the first author over a period of 45 years (Sewell, 1976). Despite its original name, PDE2D solves nonlinear systems of time-dependent and steady-state PDEs, and linear eigenvalue systems, in 1D intervals, completely general 2D regions with curved boundaries, and a wide range of 3D regions that can be represented by

Bottom drag dissipation results

In shallow water oceans, bottom drag might arise due to turbulent flow interacting with a bottom surface. If an ice cap is modelled, there is also top drag. However, there are no constraints for the drag coefficient in icy satellites, it might depend of the surface roughness and Reynold's number (Massey, 1979; Turcotte and Schubert, 1982; Sohl et al., 1995). On Earth, the bottom drag coefficient for model tidal heating is commonly assumed to range from CD = 0.002 to CD = 0.003 (Lambeck, 1980;

Rayleigh dissipation results

Fig. 4 shows the ocean height variation η, at 0, 25, 50 and 75 percent of the way through the last of 10 periods (sufficient for convergence to a periodic solution), for eccentricity forcing on Enceladus, with CD = 0, α = 10−5/sec, h0 = 500 m. Notice the problem is now linear (CD = 0), and so we are considering Rayleigh dissipation, not bottom drag dissipation. Our results are in good agreement with Hay and Matsuyama (2017), who used the same parameters.

For the linear problem we still have to

The effect of an ice cap

Beuthe (2016), Hay and Matsuyama (2019), and Rovira-Navarro et al. (2020) model the effect of an ice cap by adding a term ∇q/ρocean to the left side of (2), and a new equation, Ehi(R22η + 2η) = R2(R22q+2q – (1+ν)q), which relates the ocean height variation η and the pressure q on the ocean surface from the ice cap. This equation (B.2 in Beuthe (2016) and equation (9) in Rovira-Navarro et al. (2020)) was derived in Beuthe (2008). Here E, ν and hi are the Young's modulus, Poisson ratio and

Accuracy

To check the accuracy of our calculations, we tried various spatial grids and time steps, to see how much our estimates for the total dissipated energy changed, for eccentricity tidal forcing on Titan and parameters CD = 0.1, h0 = 100 m. This choice of parameters ensured that convergence occurred each test within 10 periods, and the results, shown in Table 5, suggest that  =  = 9°, used throughout this paper (so that programs will run with the free version of PDE2D), are sufficiently small

Conclusions

In this paper we use a finite element approach using a general-purpose software (www.pde2d.com) to calculate tidal heating in the subsurface oceans. We solved the Laplace tidal equations for three icy moons, Titan, Enceladus and Europa, and we compared our results with previous efforts, using finite differences (Hay and Matsuyama, 2017), finite volumes (Hay and Matsuyama, 2019), finite elements (Rovira-Navarro et al., 2020) and also analytical solutions (Chen et al., 2014). The numerical

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

We would like to thank Gabriel Fabien-Ouellet and four anonymous reviewers for their constructive suggestions and remarks. Especially, we like to thank anonymous Reviewer #2, for many helpful suggestions, including a way to approximate the effect of self-attraction which we used to obtain the results in Fig. 8b. All numerical computations were performed at the National Laboratory for Advanced Scientific Visualization (LAVIS) at UNAM (www.lavis.unam.mx) computing facility, and the UTEP HPC (//researchcomputing.utep.edu

References (56)

  • M.H. Moore et al.

    Mid- and far-infrared spectroscopic studies of the influence of temperature ultraviolet photolysis and ion irradiation on cosmic-type ices

    Spectrochim. Acta Mol. Biomol. Spectrosc.

    (2001)
  • J.H. Roberts et al.

    Tidal heating and the long-term stability of a subsurface ocean on Enceladus

    Icarus

    (2008)
  • M. Rovira-Navarro et al.

    Do tidally-generated inertial waves heat the subsurface oceans of Europa and Enceladus?

    Icarus

    (2019)
  • W.D. Sears

    Tidal dissipation in oceans on Titan

    Icarus

    (1995)
  • F. Sohl et al.

    Tidal dissipation on titan

    Icarus

    (1995)
  • P.C. Thomas et al.

    Enceladus's measured physical libration requires a global subsurface ocean

    Icarus

    (2016)
  • R.H. Tyler

    Tidal dynamical considerations constrain the state of an ocean on Enceladus

    Icarus

    (2011)
  • R.H. Tyler

    Comparative estimates of the heat generated by ocean tides on icy satellites in the outer Solar System

    Icarus

    (2014)
  • C. Zimmer et al.

    subsurface oceans on Europa and callisto: constraints from Galileo magnetometer observations

    Icarus

    (2000)
  • M. Beuthe

    M. Beuthe. Thin elastic shells with variable thickness for lithospheric flexure of one-plate planets

    Geophys. J. Int.

    (2008)
  • M. Beuthe et al.

    Enceladus' and Dione's floating ice shells supported by minimum stress isostasy

    Geophys. Res. Lett.

    (2016)
  • F. Birch

    Heat from Radioactivity

    (1954)
  • G. Choblet et al.

    Powering prolonged hydrothermal activity inside Enceladus

    Nat. Astronomy

    (2017)
  • Gary D. Egbert et al.

    Estimates of M2 tidal energy dissipation from TOPEX/Poseidon altimeter data

    J. Geophys. Res.: Oceans

    (2001)
  • K. Heier

    Radioactive elements in the continental crust

    Nature

    (1965)
  • L. Iess et al.

    The tides of Titan

    Science

    (2012)
  • L. Iess et al.

    The gravity field and interior structure of Enceladus

    Science

    (2014)
  • M.G. Kivelson et al.

    Galileo magnetometer measurements: a stronger case for a subsurface ocean at Europa

    Science

    (25 Aug 2000)
  • View full text