Seismic wave propagation in icy ocean worlds

TL;DR

This paper develops a seismological framework for icy ocean worlds like Europa and Enceladus, using waveform analysis and simulations to classify seismic waves and explore interior properties relevant to habitability.

Contribution

It introduces a new seismic wave classification scheme and adapts seismological concepts specifically for icy ocean worlds, aiding future planetary seismology efforts.

Findings

Developed a concise naming scheme for seismic waves on icy worlds

Analyzed seismic wavefield features on Europa, Titan, Ganymede, and Enceladus

Simulated seismic measurements to constrain interior geochemical parameters

Abstract

Seismology was developed on Earth and shaped our model of the Earth's interior over the 20th century. With the exception of the Philae lander, all in situ extraterrestrial seismological effort to date was limited to other terrestrial planets. All have in common a rigid crust above a solid mantle. The coming years may see the installation of seismometers on Europa, Titan and Enceladus, so it is necessary to adapt seismological concepts to the setting of worlds with global oceans covered in ice. Here we use waveform analyses to identify and classify wave types, developing a lexicon for icy ocean world seismology intended to be useful to both seismologists and planetary scientists. We use results from spectral-element simulations of broadband seismic wavefields to adapt seismological concepts to icy ocean worlds. We present a concise naming scheme for seismic waves and an overview of the features of the seismic wavefield on Europa, Titan, Ganymede and Enceladus. In close connection with geophysical interior models, we analyze simulated seismic measurements of Europa and Titan that might be used to constrain geochemical parameters governing the habitability of a sub-ice ocean.

Figures

Tables

Table 1. Table 1: Properties of ocean worlds considered here. a Schubert et al. ( 2004 ) b Thomas ( 2010 ); Iess et al. ( 2014 ) c Jacobson et al. ( 2006 ); Iess et al. ( 2012 ) d Vance et al. ( 2014 ) e Vance et al. ( 2017b )

	Radius	Bulk density	Moment of Inertia	ice thickness	ocean depth	high-P ice?
	km	(kg m-3)		km	km
Europaa	1565$\pm 8$	2989$\pm$46	0.346$\pm 0.005$	5–30	¿100	no
Ganymedea	2631$\pm 1.7$	1942.0$\pm 4.8$	0.3115$\pm 0.0028$	50–150d	¿100d	yesd
Enceladusb,c	252.1$\pm$0.2	1609$\pm$5	0.335	10–55e	12–50e	no
Titanc	2574.73$\pm 0.09$	1879.8$\pm 0.004$	0.3438$\pm 0.0005$	45–120e	¿80	maybe

Table 2. Table 2: Examples of phase names in our proposed naming scheme, the scheme proposed by Knapmeyer ( 2003 ) and the one proposed by Lee et al. ( 2003 ) . The latter is only applicable to local phases that did not enter the rocky part below the ocean. Surface-wave-like phases like the Crary phase cannot be named in any of the other schemes, only TauP can describe it by its phase velocity.

Our scheme (short)	our scheme (full)	Knapmeyer	TauP	Lee
Pn	PFPFP	POPoPopOp	P20P136P136P20P	N/A
Pe	P	P	P	P
PoP	PFoFP	POPo+POP	P20Pv136P	PCP
PnPn	PFPFPPFPFP	POPoPopOpPOPoPopOp	P136PP136P	N/A
PFSFP	PFSFP	POPoSopOp	P20P136S136P20P	N/A
Pf	PFoFeFoFP	POPo+pO-Po+pOp	P20Pv136p^20Pv136p20p	N/A
C	C	N/A	2950kms	N/A

Table 3. Table 3: Wave types encountered in ice floating on water. v P , v S subscript 𝑣 P subscript 𝑣 S v_{\mathrm{P}},v_{\mathrm{S}} are the wave velocities in ice, which are assumed to be constant over the depth. ρ ice , ρ water subscript 𝜌 ice subscript 𝜌 water \rho_{\mathrm{ice}},\rho_{\mathrm{water}} are the densities of ice and water The definition of symmetry is analog to Lamb waves, i.e. symmetry with respect the the center axis of the ice layer.

Wave type	frequencies	polarization	symmetry	group vel. $U$	phase vel. $c$
Love wave	high	transversal	antisymmetric	$v_S\sqrt{1-{\left(\frac{n{v}_S}{2fd}\right)}^2}$	$v_S/\sqrt{1-{\left(\frac{n{v}_S}{2fd}\right)}^2}$
Toroidal mode	very low	transversal	symmetric	$\approx v_S$	$\approx v_S$
Crary wave	monochromatic	radial	n/a	$\approx v_P\sqrt{1-v_S^2/v_P^2}$	$v_P$
Longitudinal	low	radial/vertical	symmetric	$\approx 2v_S\sqrt{1-v_S^2/v_P^2}$	$2v_S\sqrt{1-v_S^2/v_P^2}$
Flexural wave	low	vertical	antisymmetric	$c-\lambda \frac{dc}{dk}$	$v_S\sqrt{\frac{8{\rho }_{\mathrm{ice}}{\left(kd\right)}^3}{3{\rho }_{\mathrm{water}}}\frac{1-v_S^2/v_P^2}{1+2kd{\rho }_{\mathrm{ice}}/{\rho }_{\mathrm{water}}}}$
Rayleigh wave	high	radial/vertical	n/a	$0.9194v_S$	$0.9194v_S$

Table 4. Table 4: Distance ranges of seismograms regimes on different icy worlds. Epicentral distances between earthquake and receiver are noted in degrees. The number in parentheses notes the fraction of the total planetary surface that is covered by this distance regime. If we assume a random distribution of events on the surface, this is also the probability that any given event will fall into this distance range from a lander. For comparison, the distance ranges on Earth are also noted, where ”ocean shadow” corresponds to ”far regional earthquakes”.

Moon		regional	ocean shadow	teleseismic
Europa	5 km ice	$0-5^{\circ }$ (0.2%)	$5-{30}^{\circ }$ (6.5%)	$30-{150}^{\circ }$ (86.6%)
Titan	33 km ice	$0-{20}^{\circ }$ (3%)	$20-{55}^{\circ }$ (18.3%)	$55-{160}^{\circ }$ (75.7%)
	114 km ice	$0-{50}^{\circ }$ (17.9%)	- (0%)	$50-{150}^{\circ }$ (75.4%)
Ganymede 104 km		$0-{60}^{\circ }$ (25%)	- (0%)	$60-{120}^{\circ }$ (50%)
Enceladus 15 km			$5-{80}^{\circ }$ (41.1%)	$80-{180}^{\circ }$ (58.7%)
Earth		$1-{10}^{\circ }$ (0.75%)	$10-{30}^{\circ }$ (5.9%)	$30-{105}^{\circ }$ (56.4%)

Table 5. Table 5: Titan model used for detectability test of an IceVI layer as shown in fig. 11

Layer	thickness [km]	vP [km/s]	vS [km/s]	$\rho$ [kg/m3]
Iceh	46	3.9	2	925
ocean	410	1.55 – 2.55	–	1020 – 1170
IceVI	5.5	4.4	2.26	1394
mantle		7.99	4.54	3526

Table 6. Table 6: Potential scientific objectives of a future icy ocean world seismometer and seismic observables that address it.

Scientific objective	Seismic observable	signal strength	distance range	reference
Tectonic activity	Location of seismic events	strong	global	appendix A
Ice thickness	Resonant frequency of Crary phase	strong	global	fig. 8
	Transition frequency between Rayleigh and flexural surface wave	intermediate	global	fig. 6
	Reverberations in body wave coda	strong	teleseismic	fig. 7
	Autocorrelation of ambient noise	strong	ambient noise	(Panning et al., 2017)
Ocean depth	Reverberations in body wave coda	strong	teleseismic	fig. 7
	Autocorrelation of noise	intermediate	ambient noise	(Panning et al., 2017)
Ocean chemistry (from sound velocity in water)	Reverberations in body wave coda	strong	teleseismic	fig. 7
High pressure ice phases	Coda of Sn-waves	strong	teleseismic
fig. 11
	Scholte waves	strong	teleseismic
	P-to S ratio	intermediate	teleseismic, depends on focal mechanism	fig. 11
Core diameter	Autocorrelation of noise	weak	ambient noise	(Panning et al., 2017)
	Autocorrelation of seismogram	intermediate	$>100$ deg	fig. 7

Table 7. Table 7: Properties of the AxiSEM simulations for this paper. Since AxiSEM is a 2D method, the numerical cost scales roughly with the inverse cube of the minimum period. NCPU is the optimal number of CPUs for mesh decomposition; nelement the number of elements of the SEM mesh. nsteps is the number of elements times the number of time steps for a 1 hour simulation. This is directly proportional to the computational cost, which is shown here on Piz Daint, a Cray XC40 system at the Swiss National Supercomputing Centre CSCS. For comparison, the cost of a 1 hour simulation of Mars (model DWAK in Clinton et al. ( 2017 ) ) and Earth (PREM) are shown. These costs are based on using gfortran 6.1 and may be lower for more optimized compilers

Body	Radius	Period	nelement	time step	nsteps	cost
	km	s		ms	${10}^9$	CPUh
Enceladus		0.25	1406304	2.73	1854.5	3999
		0.5	355224	5.6	228.4	492
		1	90630	11.4	28.6	62
Europa	1565	1	2139072	8.88	867.4	1871
(5km ice)		2	546000	15.2	129.1	278
Europa		1	2047690	8.95	824.0	1777
(20km ice)		2	521136	15.4	121.6	262
Titan	2574	1	7663680	12.6	2189.5	4722
(33km ice)		2	1929600	25.2	276.1	595
Titan		1	7416318	9.94	2687.3	5795
(124km ice)		2	1873808	19.9	339.6	732
Ganymede	2631	1	8281306	11.5	2584.6	5574
		2	2090130	23	326.8	705
Mars	3396	1	5098080	7.7	2371.5	5114
Earth (PREM)	6371	1	9520000	7.7	4436.8	9568

Equations

$$\frac{Q_{\mu }}{{\omega }^{\gamma }}$$

$$H(P)$$

$$v_L=2v_S\sqrt{1-{\left(\frac{v_{\mathrm{S}}}{v_{\mathrm{P}}}\right)}^2},$$

$$U(f)=\frac{v_{\mathrm{S},\mathrm{ice}}}{\sqrt{1-{\left(\frac{n{v}_{\mathrm{S},\mathrm{ice}}}{2fd}\right)}^2}},$$

$$f\min =\frac{n{v}_{\mathrm{S},\mathrm{ice}}}{2d}\approx \frac{n}{d}\mathrm{mHz},$$

$$f_{\mathrm{Cr}}=\frac{(n+1)v_{\mathrm{S}}}{2d\sqrt{1-{\left(\frac{v_{\mathrm{S}}}{v_{\mathrm{P}}}\right)}^2}},$$

$$x=\frac{t_S-t_P}{1/v_S-1/v_P}\approx (t_S-t_P)\cdot 4km/s.$$

$$\Delta =180^{\circ }\left(1-\frac{t_2-t_1}{t_3-t_1}\right).$$

Full text

1 Introduction

Most of the larger icy moons of the solar system’s gas planets are considered to harbor liquid oceans below their surface (Nimmo and Pappalardo, 2016). These oceans have been predicted from the estimated internal energy budget by radioactive decay and tides (Lewis, 1971; Cassen et al., 1979) and are very likely to form in a multi-moon situation around a giant planet. In specific settings, even small worlds like Enceladus can have enough energy available to sustain a liquid ocean at least temporarily. The existence of oceans at Europa, Titan and Enceladus has been identified using density structure inferred from gravitational moment-of-inertia measurements (Schubert et al., 2004), tectonic interpretations of the ice thickness, and electromagnetic effects consistent with conductive fluids beneath the ice (Kivelson et al., 2000; Béghin et al., 2010). Due to the low density contrast between ice and water and ambiguity of electrical conductance measurements, these methods alone cannot constrain the thicknesses of ice and ocean layers uniquely, which are crucial to the ocean chemistry and therefore habitability. Especially in the large worlds Ganymede, Callisto and Titan, high-pressure ice phases at the ocean bottom would could limit water-rock interaction which regulate redox and ocean composition(Vance et al., 2014, 2017b), with a potentially negative effect on habitability. Briny fluids (Hogenboom et al., 1995) or ice convection might mitigate this effect (Choblet et al., 2017), nevertheless the presence, structure and dynamics of these ices is therefore of great scientific interest.

The thickness of a surface ice layer less than a few kilometers can be estimated from radar measurements, as planned for the upcoming ESA JUICE and NASA Europa Clipper missions (Bruzzone et al., 2013; Phillips and Pappalardo, 2014). Attenuation in the water will make it difficult to measure the ocean depth with radar; seismological measurements would complement the existing methods. Where radar is sensitive to electrical conductivity, and crystal orientation fabric (COF), seismic waves are influenced by elastic parameters and indirectly by material chemistry, COF, and temperature (Diez, 2013). A detailed discussion of geophysical measurements linked to habitability can be found in the companion paper by Vance et al. (2017b).

The seismology of an icy ocean world differs from that of Earth. Many of the classical concepts and terminology of seismology become meaningless with the addition of a thick global ocean and ice covering. The global ocean, with its low seismic velocity and absence of S-waves, partially decouples the seismic wavefield in the ice and the underlying solid interior, creating seismic wave types that are not observed at a global scale on Earth. The coming decades will likely bring the first seismological station on an icy moon [e.g.][](Hand et al., 2017), and mission design begins now. This paper gives an introduction to the general features of icy ocean world seismology, which should help terrestrial seismologists adapt their vocabulary and thinking to this environment. With the possibility of landers on Europa and Titan in the coming decades, these two worlds are discussed in more detail. Enceladus is also discussed, even though there is currently no strong proposal for a lander and its small, irregular size distinguish it from the other ocean worlds. Ganymede is also discussed as a target for a possible future lander (Vijendran et al., 2010), and as end-member in size and ocean thickness. To keep the paper concise, we will not discuss Callisto and potential ocean worlds in the Uranus and Neptune systems in detail, even though the general concepts presented here apply to them as well. For the four worlds, potential seismic tests for habitability are presented and tested using synthetic seismograms.

Previous studies of seismology on Europa, mostly older than a decade, were therefore limited by the computational methods of the time. The detailed study of (Kovach and Chyba, 2001) was based on ray-theoretical methods and analytic derivations and did not calculate any full seismograms. Lee et al. (2003) used ray theory to discuss body waves and surface waves in the ice and ocean layer of Europa, Panning et al. (2006) used the normal mode method (Woodhouse, 1988) to calculate seismograms, but were limited to long periods (¿10s), where some of the ice layer effects are not evident. Computational seismology has progressed rapidly in the last decade (for an overview, see Igel, 2017), enabling full modeling of the broadband global wavefield for the first time now. We used the spectral-element solver AxiSEM (Nissen-Meyer et al., 2014), which is extremely efficient for layered models, since it reduces the 3D wave propagation to a 2D problem. A second companion paper (Panning et al., 2017), couples these results with a seismicity model to estimate the seismic background noise caused by tidal cracking on Europa.

Any realistic lander mission on an icy world will face severe constraints in terms of seismometer mass: Generally, the noise level of a seismometer is controlled by the size of its proof mass: Smaller masses mean higher eigenfrequencies (above 1 Hz) and a higher noise level at long periods. Classical broadband seismometers, like the InSight very Broadband (VBB) main instrument (Lognonne et al., 2012; Mimoun et al., 2012; Lognonné et al., 2016), have a very good signal-to-noise ratio over frequency range from below 0.01 Hz to 10 Hz, but have total masses of several kg. Also, their sensitivity to shock and tilted installation impacts mission planning, which is prohibitive as a secondary payload on a lander. Geophones, as proposed for the Mars96 mission (Lognonné et al., 1998) are very lightweight, passive instruments, but have a low sensitivity below one Hz and will only record local events and structure. The SESAME instrument on the Philae lander (Seidensticker et al., 2007) used piezoelectric sensors, sensitive in the mHz range to see structure on a meter scale around the landing site on 67P/Churyumov–Gerasimenko (Knapmeyer et al., 2016). The second instrument on InSight, SP for short period, is a MEMS-based seismometer with a flat response to 120 s period and a noise level slightly higher than that of a Trillium compact borehole instrument on Earth (Pike et al., 2016). Therefore, the report of the Europa Lander Science Definition Team (Hand et al., 2017) defined the baseline instrument as equivalent to the Short Period (SP) instrument of the InSight Mars lander.

Figure 1 shows a seismogram from Earth, compared with synthetic seismograms for an identical source-receiver configuration on Europa, Enceladus, Titan and Ganymede. This comparison highlights the diverse waveforms that are created by the icy-ocean environment. The measured seismogram on Earth is different from the PREM reference synthetic (Dziewoński and Anderson, 1981), since the surface waves traveled along a path through thick continental crust, while PREM contains an oceanic, thin-crust model. No simple global reference model can explain the complete shape of the seismic waveform. Also, ocean-generated microseismic noise appears in the data, with periods of 7 seconds, and another noise signal at around 100s period. This may act as a rough guide for the differences one could expect between simulations and real data from icy worlds.

Seismograms on icy worlds will contain exotic phases, especially the prominent long-period flexural waves in the ice and a strong, non-dispersive longitudinal phase on the R-component (both of which were identified and explained in floating terrestrial ice sheets by Press and Ewing (1951)). On Enceladus, the Rayleigh (on Z and R) and Love (on T) surface waves orbit the moon in a relatively brief 800 seconds, and they show little energy decrease due to the small size of the object. The seismogram is therefore dominated by its repeated occurrences, which permit determination of the distance between an ice-quake and the receiver (for an application to Mars, see Panning et al., 2015). On Europa, the ice layer is so thin relative to the size of the moon that amplitude and shape of the seismogram depends very strongly on the thickness of the ice layer. On Ganymede, with a 116km ice layer, the overall seismogram shares many features with the Earth: Body waves with low amplitudes arrive first, followed by a complex surface wave train, but very weak exotic phases.

The paper starts with a description of the models used. In section 3, the seismic wavefield is described, a concise naming scheme is introduced and the various surface wave types in the ice layer are discussed. Section 4 applies the methodology to Europa, Titan, Ganymede and Enceladus, and identifies potential seismic measurements on this representative set of icy ocean worlds.

2 Models and Modeling

2.1 Planetary models

The structural models used here were adopted from the companion paper Vance et al. (2017b). As described in detail there, the models include up-to-date properties for ices, saline oceans, rocky interiors and iron cores, and are thermodynamically self-consistent within the limits of available equations of state. Radial structures are computed as per Vance et al. (2014), with self-consistent ice and ocean thermodynamics, using boundary conditions of surface and ice-ocean interface temperature. Thermodynamics for rocks have been added as per Cammarano et al. (2006). All modeling tools are freely available via GitHub (http://github.com/vancesteven/PlanetProfile).

Temperature profiles can be tuned to produce different ice shell thicknesses. It is also possible to consider a range of internal compositions and temperature profiles for the interior below the ocean. We consider the range of plausible ice thicknesses inferred from isostatic surface features, gravitational moment of inertia constraints, and possible radiogenic and tidal heat inputs.

For basic estimates of attenuation, we followed the approach of Cammarano et al. (2006) to obtain temperature and frequency dependent estimates of shear quality factor, $Q_{\mu}$ with the expression

[TABLE]

in which $B_{a}=0.56$ is a normalization factor, $\omega$ is the seismic frequency, exponent $\gamma=0.2$ is the frequency dependence of attenuation, and $R$ is the ideal gas constant. $H$, the activation enthalpy, scales with the melting temperature $T_{m}$ and with the anisotropy coefficient $g_{a}$, and the values of $g_{a}$ chosen for various ices are described in Vance et al. (2017b). The bulk quality factor, $Q_{\kappa}$, is neglected. Values of $Q_{\mu}/(\omega^{\gamma})$ in the ice range from $10^{7}$ in the colder upper regions of the ice, to $10^{2}$ approaching the freezing point, consistent with measurements in glaciers on Earth (Peters et al., 2012).

It is important to note that the mineral composition and structure of the rocky interior may be very different from Earth’s. For simplicity, we will use the term ”mantle” for the rocky layer between water and a potential iron core, which may contain parts compositionally similar to terrestrial crust or even contain unknown minerals created by exotic pressure-temperature regimes. The reader can find relevant literature in the companion paper (Vance et al., 2017b).

2.2 Wavefield modeling

This article deals with modeling the full wavefield of a layered, i.e. spherically symmetric icy world. To this end, the spectral-element solver AxiSEM was used (Nissen-Meyer et al., 2014). It separates the problem of wave propagation in a cylindrically or spherically symmetric object into an analytical solution of the problem in the azimuthal ($\varphi$) direction perpendicular to the source-receiver plane, and a numerical spectral-element discretisation within the in-plane $r,\theta$, which reduces the numerical cost to that of a 2D method (Nissen-Meyer et al., 2008, 2007) and includes attenuation (van Driel and Nissen-Meyer, 2014a) and anisotropy (van Driel and Nissen-Meyer, 2014b). The software was modified to handle general 1D velocity profiles. The changes are included in the recent release 1.4. The properties of the SEM meshes and the resulting numerical cost are shown in appendix B Compared to the typical meshes for global wave simulation on terrestrial planets, smaller elements are needed due to the low wave velocities in the ocean. As these are fluid elements, the maximum allowed time step is larger than on earth, so that the computational cost is about a factor 2 larger than for an terrestrial planet of similar radius.

The reciprocity of the Green’s function permits switching of the locations of source and receiver of a seismic wavefield. This has been used in the Python package Instaseis (van Driel et al., 2015), which allows to reconstruct seismograms for arbitrary source and receiver locations from a precalculated wavefield database. Instaseis uses the stored displacement wavefield to calculate strain, which allows to simulate arbitrary moment tensor and single force point sources with arbitrary source time functions. Theoretical arrival times were calculated with TauP software package (Crotwell et al., 1999) in the implementation of ObsPy 1.0 (Krischer et al., 2015).

To estimate the effects of three-dimensional heterogeneities and scattering, we introduced lateral heterogeneities into some runs. Since AxiSEM assumes an axially-symmetric model, these heterogeneities have to be symmetric around the source axis. This creates a limited scattering effect, especially since no off-path scattering is simulated, but it serves as a good first-order approximation to estimate the total effect of heterogeneities. For the scattering, we implement a von Kármán random medium with a correlation length of 5 km and velocity variations of 10%. This means that the correlation length is similar to the shortest wave length of P-waves. In models with lateral heterogeneities to simulate scattering, these structures become symmetric around the receiver axis.

Input models for all planets discussed in this article are available as an electronic supplement. The seismic wavefields for the models discussed here are available as Instaseis databases on http://instaseis.ethz.ch/icy_ocean_worlds/ and can be used freely for further analysis. Waveform databases for Earth models can be accessed from the IRIS Syngine service (Krischer et al., 2017).

3 The seismic wavefield

3.1 Phase naming

[TABLE]

To facilitate seismic analysis, it is necessary to be able to name phases (i.e. distinct signals in the seismogram) unambiguously. The terrestrial naming scheme adopted by IASPEI (Storchak et al., 2003) needs to be extended in the context of icy moons to account for the various ice phases. A completely non-ambiguous scheme, which marks every interface crossing, was proposed by Knapmeyer (2003), but since it produces rather cumbersome phase names, we propose a simplified version here. Lee et al. (2003) proposed a scheme specifically for local waves, but since it cannot be generalized to global wave propagation, we decided not to use it here. In our scheme, all names from the IASPEI scheme keep their meaning, where applicable, but we add a few new letters for ice-specific phases. Because “i” and “I” have already been used for inner core phases, we propose “e” as the letter generally related to ice phases, using the German word Eis (analogous to core phases being called “K” from Kern).

The main issue at hand is that there is only the Moho as a strong, global near-surface interface on Earth, but there are at least 3 in ocean worlds (a Moho, plus the ice-bottom and the ocean floor and possibly an interface between high-pressure ice and the crust). Standard IASPEI nomenclature would denote these interfaces by their depths. Since this depth is poorly constrained on ocean worlds, we propose to give these interfaces descriptive letters: “e” for the bottom of a surface ice layer, “o” for the bottom of an ocean, “g” for the top of the rocky layers and “m” for the top of the mantle. To use this scheme for glaciers on Earth or worlds like Ceres, the bottom of the ice layer should also be called “e” if it is in contact with a solid rocky interior, see 3 for an overview.

The IASPEI standard for marking reflections on these interfaces (appended “+” or “-”) is not widely used in the community, so we stick to the nomenclature of Aki and Richards (2002), where topside reflections are noted by the interface descriptor (or depth), as in PeP, and underside reflections are noted by the descriptor as a superscript (e.g. PmP). Where the usage of a typographic feature is problematic, the mnemonic TauP convention of P^mP should be used instead.

Terrestrial seismology calls compressional waves in the ocean either ”H”, for a hydroacoustic wave from a source in the water, or ”T”, for a hydroacoustic wave from a source in the ground. Both wave types are linked to the SOFAR channel of low sound speed between 500 and 1500 meter depth in Earth’s ocean (SOund Fixing And Ranging; Tolstoy and Ewing, 1950; Ewing et al., 1951). This channel is an effect of the specific contributions of depth-dependent salinity and pressure to sound speed in earth’s ocean and most likely does not exist in ocean worlds. We therefore propose to mark a leg of the path that traverses the ocean generally with “F” (fluid). The terrestrial “H”-wave would then be “FP” and the “T”-wave “pFP”.

Waves in the uppermost ice layer and in the mantle are both called P or S, depending on their type. A PFPFP wave is therefore one that started within the ice layer, crossed the ocean, traveled through crust and potentially mantle, turned there and went back through the ocean to the surface. To retain a certain elegance in naming phases, we propose to use “Pn” for this wave type, analogous to the mantle phases on Earth as proposed by (Mohorovičić, 1909). To specifically highlight local phases that never left the ice, we propose to add an “e”: “Pe”, “Se”. Note the distinction from “PeP”, which was reflected at the ice bottom.

To keep phase names concise, we introduce a few more abbreviations: “Po” is short for a topside reflection from the ocean floor (full: PFoFP), “Pm” is short for a topside reflection from the mantle top (full: PFPmPFP). “Pf” is a single reverberation within the water (PFoFeFoFP). These abbreviations can also be used as prefixes to mark multiples on other phases, as in “PfPn” = “PFoFeFoFPFP”. See fig. 4 for an overview of ice- and ocean-interacting body wave paths. As on Earth, a number following one descriptor means a multiple reflection. “Pf2Pn” was reflected twice at the ocean bottom and twice at the bottom of the ice layer before traveling as a P-wave through the mantle.

Up-going seismic waves from an event in or below the water layer should be distinguished. This differs from the terrestrial usage, where source depth is explicitly not coded in the phase name. We therefore propose to write the first letter in lowercase for events below the ice, e.g. pFP for an event originating below the ocean layer or fP for an acoustic source in the ocean. This is a slight extension of the terrestrial usage, where only rays with an upward leg are noted as such. This naming scheme combines the brevity of the IASPEI scheme with the flexibility of Knapmeyer’s scheme, while at the same time not being limited to local waves, as in the scheme by Lee et al. (2003). See table 2 for examples of phase names in our scheme, Knapmeyer’s scheme and the input to the popular TauP software.

We show a visualization of the global wavefield in fig. 5 that follows the example of the IRIS global stack (Astiz et al., 1996). Seismograms for different distances are color coded, where vertical motion is color coded blue, motion in R direction (away from the event) in green and transverse motion in T. Each trace is normalized separately. We see that in the coda of Pn-waves (see also fig. 7), every body wave phase is followed by multiples from the mantle-ocean and ocean-ice boundaries.

3.2 General features of the icy-moon wavefield

3.2.1 Seismic phases in the ice

Seismic events within the ice layer will create a range of seismic wave types that are not observed in global seismograms on Earth (see fig. 6 for a spectrogram, which separates the spectral characteristic better; see Table 3 for a summary of the different wave types, including their group and phase velocities.). This is especially the case for ice layers thinner than 40 km, as expected for Europa. The three most prominent new phases are flexural waves, longitudinal waves, and the Crary phase. As shown by Press and Ewing (1951), the response of a floating ice layer of thickness $d$ is similar to that of an infinite slab. Kovach and Chyba (2001) applied this to the situation on Europa, so that we can summarize the most important results here, assuming a compressional wave velocity $v_{\mathrm{P,ice}}=4~{}\mathrm{km/s}$ and a shear velocity $v_{\mathrm{S,ice}}=2~{}\mathrm{km/s}$:

Rayleigh waves:

Short period surface waves with vertical/radial polarization do not see the bottom of the ice layer and can therefore be described as a Rayleigh wave with the usual properties, if the inverse wavenumber $k^{-1}$ is smaller than the ice thickness $d$: ($kd=2\pi\frac{d}{\lambda}>1$).

The velocity gradient with depth is very small in the ice, so Rayleigh waves have little dispersion. For small wavelengths, the usual result of a phase velocity $c_{R}=0.92v_{\mathrm{S,ice}}\approx 1.8~{}\mathrm{km/s}$ is valid.

Rayleigh wave (ice bottom):

At short periods, an interface wave travels at the ice-water interface, called Stoneley wave by Press and Ewing (1951) and Kovach and Chyba (2001) (even though that term usually describes solid-solid interface waves). An interesting property is that its maximum phase velocity is controlled by the sound velocity of the ocean: $c_{\mathrm{St}}=0.87v_{\mathrm{P,water}}\approx 1.3\mathrm{km/s}$. Its signal is only detectable for ice shells of a few kilometers, but there, the arrival time difference from the surface Rayleigh wave constrains the sound velocity of the ocean.

Flexural waves:

Long-period vertically polarized surface waves, with $kd<1$ are perturbed by the bottom of the ice layer and travel as full-layer flexural waves. The most interesting result is an intrinsic dispersion $c_{\mathrm{Fl}}\propto\sqrt{f}$ that is independent of the velocity gradient in the ice. Due to this frequency dependence, the short periods arrive first, and the group velocity of long-periods goes to zero. The cutoff frequency is where the Rayleigh regime (see above) starts. Panning et al. (2006) showed the effect of different ice thicknesses on this wave train.

Longitudinal wave:

Another long-period solution is a wave with motion in the radial direction, in phase over the whole ice layer, described analytically in Press and Ewing (1951). This wave has a phase and group velocity defined by the expression

[TABLE]

with $v_{L}\approx 3.46$ km/s for the example P and S velocities. It therefore travels as a fast, non-dispersive pulse along the surface. This means that it arrives much before all other ice phases at long distances, and even before mantle body waves at intermediate distances (see tab. 4). It is the most prominent feature in the broadband seismogram. Both flexural waves and longitudinal waves are solutions for wavelengths long relative to the ice shell thickness, and so only show up at relatively low frequencies ($kd<1$).

Love waves:

No shear waves exist in the fluid ocean, so the ice layer is a wave guide for SH waves, which are completely reflected at surface and ice-ocean interfaces. This results in a dispersive group velocity of

[TABLE]

where $n=1,2,3,...$ is the mode of the Love wave, $v_{S}$, the shear wave speed and $d$ the ice thickness. Similar to the flexural wave, we have an inverse dispersion, where infinite frequencies travel with $U=v_{\mathrm{S,ice}}$, while lower frequencies are slower. The minimum frequency of the $n$-th mode Love wave is

[TABLE]

for ice with $v_{\mathrm{S,ice}}=2~{}\mathrm{km/s}$. This means that $f_{\mathrm{min}}=0.2$ Hz for the fundamental mode in a 5 km ice shell and 0.05 Hz for a 20 km ice shell. It thereby constrains the minimum ice thickness along the path.

Toroidal wave:

For long periods below the minimum frequency of the guided Love wave, the whole ice layer moves in phase. Analogous to the longitudinal wave, this results in a non-dispersive, low-frequency wave with high amplitude and a phase and group velocity equal to the SH velocity of the ice shell.

Crary phase:

As a result of SV waves being completely reflected at critical angle for SV-P conversion on the boundaries of the ice, another phase arises with interesting properties (Kovach and Chyba, 2001, p. 282): As first described by Crary (1954), with application to terrestrial floating ice in the Arctic, it has mainly radial displacement, a phase velocity of $v_{\mathrm{P}}$ (for homogeneous ice layers) and a harmonic frequency spectrum characterized by

[TABLE]

where $n=1$ defines the characteristic frequency, and larger integer values of $n$ representing overtones. Its group velocity is $U\approx 0.75v_{\mathrm{P}}$, which means it arrives slightly after the longitudinal wave. Measurement of $f_{\mathrm{Cr}}$ constrains the ice thickness $d$, if $v_{\mathrm{S}}$ and $v_{\mathrm{P}}$ are known. Since $f_{\mathrm{Cr}}(d=5\mathrm{km})=0.46~{}\mathrm{Hz}$ and $f_{\mathrm{Cr}}(d=20\mathrm{km})=0.11~{}\mathrm{Hz}$, the thickness of Europa’s ice layer can be measured with a seismometer that is sensitive at these frequencies. For ice thicknesses above 40 km, the maximum frequency becomes too low to be observable for realistic events with a space-ready seismometer.

Especially for thin ice layers, the wavefield has two frequency regimes: Most of the phases specifically related to the ice layer (Crary phase, longitudinal and flexural waves) have a maximum frequency between 0.1 and 0.5 Hz (for thicknesses of 20 and 5 km resp.), depending on ice thickness. Body waves from the mantle are not limited in maximum frequency and should thus be clearly distinguishable (see the spectrogram in fig. 6). Scattering within the ice layer enhances this effect by damping high frequency ice body waves, with long paths in the ice layer, generally ice-body waves fade away after a short distance.

3.2.2 Water phases

Due to the large velocity contrast between compressional wave velocity in the ice and in the water, no phases like PFP exist in shallow ocean worlds like Europa, while they would be a dominant signal in the late seismogram on Ganymede, Titan. On Enceladus, they will be present, if its ice is thin enough. Attenuation in water is very small, so phases that are multiply reflected at the top and bottom of the water layer (Pf, Pf2…) are observable almost over the whole distance range. Their polarization is mainly on the radial component.

3.2.3 Crustal and Mantle phases

Body waves in the rocky interior cover the whole frequency range. On a thin-ice world like Europa, they dominate the seismic signal above 0.1 Hz. Because of this, a broadband seismometer, sensitive enough to measure in the whole frequency range, can estimate the wave type from its frequency content alone. It also means that, while displacement amplitudes are dominated by ice phases after their arrival, body waves might still be recovered from the high-frequency content of the seismogram.

The Scholte wave is one particular crustal phase occurring at the ocean bottom. This wave shares most characteristics of a Rayleigh wave, but is defined on a solid-fluid interface (Scholte, 1947). Its phase velocity is frequency-dependent, but below the shear velocity in the solid and the sound velocity in the fluid, whatever is lower (Rauch, 1980), i.e. below 1.8 km/s. It has been used in marine exploration on Earth to estimate the shear wave velocity of the uppermost layers below the sea floor (Nolet and Dorman, 1996; Kugler et al., 2005). If it can be detected, a detailed dispersion analysis might be used to constrain the existence of hydrated low-velocity layers on the sea floor. As it can be seen in Figure 5, it shows strong water multiples.

3.2.4 Distance regimes

In terrestrial seismology, several distance regimes are defined, because seismograms change their first order shape considerably from one to the next (Kennett, 2002):

Local (0 to 100 km)

For shallow events, the ground displacement is controlled by the near field of the source and significant permanent deformation occurs at the surface. The local regime can be considerably larger for mega-thrust events (Grapenthin and Freymueller, 2011).

Regional (100–1000 km)

Crustal body wave phases arrive first. Crustal S-waves and surface waves form one wave train (Lg).

Far regional (1000–3000 km)

Mantle body wave phases arrive first and are separated from the surface wave train, but due to the upper mantle discontinuities, multiple (triplicated) phases arrive within a short time window.

Teleseismic (3000–9000 km)

Mantle body wave phases are clearly separated from each other and are visible in the seismogram as well-defined sharp pulses. Direct body waves do not sense the core.

Core shadow (beyond 9000 km)

Direct P-waves are shadowed by the core. Multiple complex core phases arrive first, followed by surface-reflections (PP).

To summarize, the distance regimes define whether the first arriving phases have mainly travelled through the crust, the complex upper mantle, the relatively homogeneous lower mantle or the core. The situation on icy ocean worlds is similar. We propose to separate the following regimes:

Local

Near field effects are dominant, including permanent displacement.

Regional

Body waves that traveled through the ice layer only arrive first.

Ocean shadow

Direct body waves are masked by the curvature of the planet’s surface. The first- arriving phases are longitudinal and Crary waves and Love waves.

Teleseismic range

Mantle body waves arrive first

Core shadow

The planetary body’s core masks direct body waves.

The distance range of these regimes depends mostly on the thickness of the ice layer. An exception is the local regime, whose size depends mainly on the magnitude of the event, but also on the local velocity structure. On thin-ice worlds like Europa, considerable coseismic displacement will occur over hundreds of kilometers, but its long periods will not be detectable with a space-ready seismometer. Therefore, we do not discuss this region in detail here.

The distance regimes for a set of icy ocean worlds can be found in Table 4. For the thin-ice worlds like Europa, it is unlikely that an event will randomly fall into the regional area around a lander, where direct ice phases are well-recordable (unless a specifically active region has been chosen for the lander). Roughly 5% of randomly distributed events will be in the ocean shadow range, where special ice layer phases will be the first arrivals. The remaining events will be beyond that in the range, where mantle and core body waves arrive first.

3.3 Single-station seismology on icy moons

The first efforts at seismology on icy moons will probably be conducted from a single station. Single-station seismology was common on Earth for decades, and many fundamental discoveries, like the existence of the inner core, were derived from observations of earthquakes on single observatories (Lehmann, 1936). Other fundamental work studied one earthquake at multiple receivers, which is conceptually similar (e.g. the discoveries of the core and the Mohorovicic discontinuity by Oldham, 1906; Mohorovičić, 1909). It has long been known that the parameters of earthquakes can be determined from one station alone (Ekström et al., 1986; Wu et al., 2006). A large body of work on modern single-station seismology has been done in the context of the Mars seismometer on the InSight lander (Panning et al., 2015; Khan et al., 2016; Böse et al., 2016). The location of a marsquake is determined by the polarization of P and surface waves and the time difference between specific phases.

Because the longitudinal and the Crary wave are prominent features in the seismogram of Europa, they can be used for fast determination of the backazimuth of the event (the direction as seen from the seismometer; see appendix A for a worked example). Their energy loss due to attenuation is relatively low, so they propagate around the planet multiple times. From the time difference between these orbits, the distance can be determined. On thick-ice worlds, the Rayleigh wave can take this role. This constraint is best achieved using distant events, using the following work flow:

Estimate the ice thickness from converted waves in the coda of body waves of distant events. This analysis needs no event localization and not even a clear identification of phases. Figure 7 shows the coda of Pn waves in different Europa models. Its coda contains clearly recognizable multiples, whose time difference is the ice thickness divided by the wave velocity. The effect is clearest for very large distances (135 degree), where the incident angle is almost vertical. The horizontal component seismograms have their maximum energy several seconds after the vertical seismogram. The two lower models in fig. 7 have been modified to simulate the effect of a highly attenuative ice layer (lowQ) or of a low-attenuation, high-scattering model, similar to the lunar crust (scat). Here, the reverberations become much less pronounced, but are still visible, especially on the long-distance horizontal components. This method will work best for large distances, where Pn is a clear arrival. The 45 degree seismograms contain triplicated phases, which creates the very complex waveform. 2. 2.

With the constraint from step 1, refine the thickness estimate from the characteristic harmonic frequency of the Crary wave (for ice layers ¡ 50 km) and estimate P- and S-wave velocities in the ice. Figure 8 shows spectra of the Crary wave from simulated waveforms for two Europa models, Titan and Ganymede. A time window after the high amplitude longitudinal wave LL was selected, where the Crary wave can easily be found as a harmonic signal, even without identifying any other phases. For all realistic ice thicknesses of Europa, the Crary peaks are very prominent in the spectra of the R-component and can be used to determine the ice thickness reliably. For ice layers of more than 100 km, the peaks are less prominent and the implicit assumption of a thin ice layer with constant wave speeds is not valid anymore, as seen for the Ganymede model. 3. 3.

With an ice layer model, locate local events, using the time difference between Pe and Se phases, as well as reflections from the ice-ocean boundary (PeS or PeP), once the ice thickness has been constrained. The distance of a local event with P and S arrival times $t_{P},t_{S}$ is

[TABLE] 4. 4.

From multiple-orbit Crary, longitudinal or Rayleigh waves, locate distant events. Appendix A shows an example of determining the direction of an Europaquake using the polarization of body waves and the energy content of the horizontal channels. If minor arc ($t_{1}$), major arc ($t_{2}$), and single orbit ($t_{3}$) arrival times of any surface wave are known, the distance of the event is (Panning et al., 2015)

[TABLE]

Note that the group velocity of the wave is not needed. If it is known, only two of the arrival times are needed. 5. 5.

Identify ocean multiples in the coda of body waves of distant events and constrain the ocean depth by estimating the two-way travel time of reflected phases. 6. 6.

Identify Scholte waves at the ocean bottom to constrain the uppermost crustal structure below the sea floor from spectral peaks or dispersion curves (see Sect. 4.2 for an example on Titan). 7. 7.

Identify as many body wave arrivals as possible to constrain overall mantle velocity structure with the event locations from step 4, analogous to the method proposed by Khan et al. (2016) for Mars. 8. 8.

If direct ocean phases (PFP, SFS) have been identified together with other body waves, the speed of sound in the ocean can be estimated, constraining the ocean chemistry.

The process sketched here will be an iterative one in practice. Steps 1 and 2 can be done even in noisy seismograms without clearly identified phases or located events. Steps 3 and 4 need some identified phases and benefit from first estimates on ice thickness and ice velocities. The steps 5-8 need first estimates of event distances to work, but their results can be used to better constrain the velocity models in turn. While the process is described here only for a single event, observations of multiple events will be necessary to obtain reliable results. Recent estimates of the seismicity induced by tidal cracking on Europa (Panning et al., 2017) suggest at least one globally observable event per week. The seismicity on other ocean worlds is not yet estimated.

4 Icy ocean worlds in the solar system

4.1 Europa

Europa has been the subject of the most detailed studies to date, and a lander mission could take place in the coming decades (Hand et al., 2017). Therefore, we discuss wave propagation and measurability of seismic waves in detail here.

4.1.1 Overview

Compared to its radius, Europa has a thinner ice layer than all other icy ocean worlds. This creates pronounced ice-layer-specific seismic phases. The transition from a ice-Rayleigh wave to flexural waves occurs between 0.1 Hz (for 20 km ice thickness) and 0.5 Hz (for 5 km ice) (Panning et al., 2006), which would be well observable with a wide-band seismometer, as the SP instrument used for InSight (Pike et al., 2016) and proposed as the baseline instrument for the NASA Europa lander (Hand et al., 2017).

Apart from the thickness of the ice layer, which can be measured as described above, or from the Crary phase, key science questions on Europa are the depth and composition of the ocean, the presence of liquid or mushy pockets in the ice and the potential of material transfer between the rocky interior, the ocean and the ice layer. The depth of the ocean can be determined well from the time delay between ocean multiples (PfPn, Pf2Pn, …) in the P-wave coda of teleseismic Europaquakes. The ocean chemistry can be estimated from the sound speed in the ocean, which is increased by higher salt contents. This increases the curvature of the SFS path, resulting in higher amplitudes for SFS at short periods; a phase that would not be observable in a pure water ocean. This could be used to constrain the ocean chemistry, once the ocean depth and thickness have been determined. Liquid pockets close to the lander would result in scattering of seismic energy from certain incidence angles. Motion of liquids in the ice shell would create seismic signals similar to geysers on Earth (Kedar et al., 1998). The sharpness of the ocean-rock interface, including potential hydrated layers could be estimated from Scholte waves on the seafloor.

Due to Europa’s modest size, multiple planetary orbits of seismic ice-phases, especially the longitudinal waves, would be observable with an SP-like instrument, allowing for a straightforward inversion for the distance of a seismic event. Due to the thin ice layer, no direct SH waves are possible beyond 100 km, so that the transverse component will be very distinct from the radial component, allowing an estimation of the event backazimuth.

In a companion paper (Panning et al., 2017), the background seismicity from tidal cracking is estimated, showing that it is detectable with a seismometer. This seismicity will create a significant seismic hum, which exceeds noise floor of industry standard Trillium compact instruments (Ringler and Hutt, 2010) at periods between 1 and 10 seconds. It could be be used for ambient noise analysis techniques, to constrain the thickness of the ice layer or mantle discontinuities.

Another important seismic source will be surface tectonics. The young surface of Europa shows clear signals of subduction Kattenhorn and Prockter (2014), collapses of cavities (Walker and Schmidt, 2015) and potentially plume activity, which would all have specific seismic signatures (Vance et al., 2017a). While no atmosphere exists to create microseismic noise, the vicious circulation in the ocean (Soderlund et al., 2013) might just be strong enough to create a detectable long-period signal (Panning et al., 2017).

4.1.2 Measurability

Since a seismometer on Europa will probably have characteristics not exceeding those of a Trillium compact instrument, we tested which seismic phases could be detected for a magnitude 3.3 event at a distance of 90 degree (fig. 9). Body waves exceed the instrument self-noise by 20 dB between 0.1 and 2 Hz, and will therefore be observable. The longitudinal phase, which is crucial for event location, exceeds the self-noise by 30 dB. The late seismogram on R and Z, including Rayleigh and flexural waves, can be measured for periods up to 50s. The Love wave on T, which is confined to the ice layer and therefore strongly affected by attenuation and scattering in the ice, is not observable for most focal mechanisms. Based on the seismicity models in Panning et al. (2017), a magnitude 3.3 event can be expected to occur somewhere between once per week and once per month. Compared to Mars or Earth, an instrument on Europa will not be affected by atmospheric noise, so instrument self-noise is the reference to measure signal strength against.

4.2 Titan

Titan is predicted to have water and ice layers with total thickness around 480 km above a low-density silicate mantle. The ice covering is interpreted to be 55-80 km based on the observed Schumann resonance (Béghin et al., 2012). Following Vance et al. (2017b), we consider thicknesses in a slightly wider range between 46 and 118 km. Depending on the ocean’s temperature profile with depth, high-pressure ice phases might be present below the ocean, with a total thickness of up to 230 km or not present at all. Depending on the ocean temperature, the ice thickness may be so high that the longitudinal, Crary, and flexural waves are restricted to very long periods, unobservable with a realistic instrument. Because the rocky mantle makes up more than 80% of the diameter in all models, the teleseismic range dominates on Titan (see fig. 10).

An important question for habitability in Titan is the existence of high-pressure ice phases at the bottom of the ocean. If they are present, they may form a barrier to water-rock interactions that regulate redox and ocean composition, and which are thereby critical for supporting life. While intra-ice convection might be able to transfer ions (Choblet et al., 2017), the absence of a high-pressure ice layer could benefit habitability by providing a direct water-rock interface, potentially with life-supporting hydrothermal systems analogous to those on Earth’s seafloor. Because such a high-pressure ice layer is not detectable with radar or gravitational moment-of-inertia measurements, seismological measurements provide an unique probe. As was argued in Vance et al. (2017b), the thickness of surface ice above the ocean and high-pressure ice below are strongly correlated and both mainly controlled by the ocean ion content and temperature profile. According to these simulations, every model with a surface ice thickness of more than 50 km must also contain high-pressure ice layers. Seismic methods could constrain surface ice and high-pressure ice thicknesses independently, and thereby constrain ocean chemistry and temperature.

As one example, Figure 11 shows a comparison of waveforms for two thin surface ice models (46 km) on top of a 410 km deep ocean with 3.3 weight percent NH3 (similar to the 3 wt percent NH3 model with 264 K in table 8 of (Vance et al., 2017b)). For this surface ice thickness, high pressure ice can neither be confirmed nor excluded. The figure compares seismograms of this model with one where the iceVI layer has been replaced by mantle (see table 5). Because there is no distance where high-pressure ice phases arrive first, their presence will have to be inferred from the coda of mantle body waves. The strong velocity and density increase at the ocean-mantle interface (assuming that no iceVI layer is present) causes a strong impedance contrast, which is weakened if there is an intermediating layer. The lower two plots in fig. 11 show acceleration spectra of a window around the P-wave arrival (left gray patch in upper subfigure) and the S-wave arrival (right gray patch in upper subfigure). The amplitude of the Pn wave train in a model with an ice layer is two times higher than in the model without the layer. Since it will be difficult to reliably determine the magnitude and focal mechanism of a seismic event with just one station, this effect may go unnoticed. The spectrum of the S-waves (lower right) however, shows clear effects of the iceVI layer: While the spectrum of the simple model is almost flat between 1 and 10 seconds, the iceVI layer produces two clear peaks at 2 and 5.1 seconds, corresponding to reverberations of P and SV waves. This signal should be easily distinguishable, even without a detailed source model.

Titan is the only icy ocean world with an atmosphere. Due to the atmosphere’s low density in comparison with the solid and liquid portions, the direct effect on the seismic wave propagation in the planet will be minimal; however, it could serve as a source for seismic signals. Wind waves on Titan’s methane lakes have been predicted to reach significant wave heights of 0.2 m at periods of 4 seconds (Lorenz and Hayes, 2012) for winds of 1m/s. The handful of Cassini measurements available through spring did not show evidence of waves above a few millimeters in height e.g. (Hofgartner et al., 2014; Zebker et al., 2014) but these may simply have been on calm days. If or when 0.2m waves form, they create pressure on the sea floor of the order of 100 Pa, which is comparable to pressure variations on the abyssal sea floor on Earth due to ocean wave interference (Cox et al., 1984; Davy et al., 2014; Stähler et al., 2016). Because the latter creates a globally observable seismic signal (Kedar et al., 2008), it must be assumed that measurable ocean-generated microseisms are possible on Titan. Surface waves from these signals could be used to constrain the uppermost layers below the surface. Note, that while Titan’s seas occupy 1% of Titan’s surface area total, they cover some 12% of the terrain northward of 55 degree N (Hayes, 2016), so microseisms might be expected to be a much more prominent feature of the seismic environment in Titan’s high northern latitudes than elsewhere, and more so during the summer when winds are expected to be strongest.

Pressure changes in the atmosphere might also create seismic signals directly, as observed on long periods on Earth (Peterson, 1993; Beauduin et al., 1996). Local phenomena, such as ”dust devils” on Earth and potentially Mars (Lorenz et al., 2015) might also be observable. Seismic waves may also create atmospheric signals that might be detectable from orbit, as observed on Earth and proposed for Venus (Lognonné et al., 2016).

4.3 Ganymede

The density structure inferred from moment of inertia measurements suggests that Ganymede has water layers in different forms, with total thickness of about 800 km above a rocky mantle and an iron core (Vance et al., 2014). The uppermost ice layer is assumed to have a thickness of more than 50 km, floating on a liquid ocean of more than 100 km depth. Due to the high gravity and abundance of water, pressure at the ocean floor creates layers of ice V and VI 10s to 100s of km thick, between the liquid ocean and the rocky mantle. Similar to Titan, the existence and thickness of these layers can only be measured by seismology. Since both are controlled by ocean chemistry and temperature, this creates a direct geophysical observable for habitability conditions of the ocean (Vance et al., 2017b).

The thick ice layer means that longitudinal and Crary phases are relatively weak and constrained to long periods, which makes them unobservable without a true high-sensitivity broadband seismometer. The ice thickness therefore would have to be constrained from multiples in the coda of body waves. Compared to Titan, the thickness of the combined water layers is 30% greater, which clears the coda of body waves. As fig. 12 shows, the wavefield is nevertheless dominated by water multiples. The analysis presented in the previous section to detect high-pressure ice phases could be done in a similar fashion for Ganymede, though the high-pressure ice layers are much thicker and more heterogeneous than on Titan which will shift the resonances to longer periods.

Ganymede’s interior is likely to be at least as seismically active as that of Earth’s moon, where continuous seismicity was observed by Apollo instruments (Vance et al., 2017b). A major question at Ganymede involves the presence of a liquid iron core, which seems to be required by Galileo observations of an intrinsic dipolar magnetic field (Kivelson et al., 2002). Melt could be present in the overlying rock as well. This study necessarily focuses on the initial problem of identifying radial boundaries closer to the surface, but a study focused on the goal of long-term exploration of Ganymede would need to evaluate the ability to probe the deeper interior to understand Ganymede’s composition and current thermal state.

4.4 Enceladus

Enceladus’ global subsurface ocean has been confirmed from its measured shape (McKinnon, 2015) and libration (Thomas et al., 2016), but the ocean thickness probably varies globally, exceeding 15 km at the south pole due to a thinner ice shell, and less than 10 km everywhere else. For the simulation, we assumed a constant global thickness of 10 km.

Due to its small size, geometrical spreading and attenuation of seismic energy plays a small role on Enceladus. As fig. 1 and the global wavefield stack in fig. 13 shows, even for magnitude 3 events, more than 5 orbits of longitudinal phases are observable, which would allow for easy determination of ice thickness and event locations. Tidal heating from the resonance with Dione suggests a high energy budget for seismicity, consistent with the active tectonic deformation of the South Pole region (Porco et al., 2006).

The plumes at Enceladus’ south pole (Porco et al., 2006; Hansen et al., 2006) are specifically caused by geysers or cryovolcanism. Since geysers are sources of seismic signals on Earth (e.g. Kedar et al., 1998), the south pole region of Enceladus can be expected to have a very high seismic background noise. At the same time, the subterranean volcanic system means that the seismic velocity structure will be very heterogeneous, which may mean that the concepts of global seismology derived from spherically symmetric models presented here are not applicable.

In addition to these sources, Enceladus should also regularly generate seismic signals in the ice as a result of tidal flexing (Vance et al., 2017a). The combination of these sources merits further investigation building on the study of Panning et al. (2017).

5 Discussion

This study gives an overview of general characteristics of icy moon seismology on a global scale, but is in no way exhaustive. An obvious limitation is that only spherical models were studied so far. The above-mentioned methods to determine the ice thickness are likely sensitive to ice undulations (e.g., Nimmo et al., 2007; Lefèvre et al., 2014). The characteristic frequency of the Crary wave is sensitive to the average thickness. The Love wave cutoff is sensitive to the minimum thickness, while body wave coda measures the thickness at the location of the source and the receiver. Enceladus for example has an ellipticity of 1/50, and probably a variation in ice and ocean thickness of more than 50% between the south pole and the equator (Thomas et al., 2016), which will strongly affect propagation, especially of multiply orbiting surface waves. From simulations using cylindrically symmetric models, we found that the longitudinal wave is a very robust feature and multiple orbits are detectable, even with strong heterogeneities in the ice. Undulations of the ice bottom may, however, have a very strong effect on the amplitude and characteristic frequency of the Crary wave and therefore on its effectiveness in constraining ice thickness.

On Europa and Enceladus, the existence of plumes suggests that liquid channels exist in the ice that will strongly scatter seismic waves. Therefore, full-3D simulations will be necessary as a next step.

Seismic wave propagation on icy ocean worlds is a numerically challenging problem, since it requires high frequencies of up to one Hz to resolve seismic features like the Crary wave for thin ice layers. Since the numerical cost in full 3D methods typically scales with frequency $f$ as $O(f^{4})$, this prohibits their application on planet scales. This is problematic, given that the ice waves may be very sensitive to three-dimensional structure, like varying ice thickness. This is rather different from the situation on Earth, were surface waves are strongest at longer periods (typically above 15 seconds), which are possible to simulate on large HPC systems. Seismic wavefield solvers suitable for smooth three-dimensional structure at planetary scales at frequencies of up to one Hz are just becoming available (Leng et al., 2016) and will offer more insight into the question which phases are usable. Local structure around a lander may strongly affect the measured seismic waveforms, which may be simulated by including local 3D wavefield simulations into global 1D simulations.

This article is intended to set the stage of global seismic wave propagation in icy moons on which more detailed studies of effects of three-dimensional structure can build.

Since the presented analysis of the seismic wavefields is in all likelihood incomplete, we make the Instaseis databases available at http://instaseis.ethz.ch/icy_ocean_worlds/. The scripts to reproduce the figures are available on the seismo-live website https://www.seismo-live.org.

Future research should also explore the parameter space of possible mantle compositions and thermal states. As was shown here, seismology in principle provides tools to constrain both, but detailed studies are necessary for each planetary body, in close collaboration with planetary geophysics (see table 6 for a brief summary).

6 Conclusion

In this paper, we discussed the general characteristics of global seismology on icy ocean worlds. The existence of the liquid-filled gap between the mantle and the icy crust creates a different wavefield than the one that is known on earth. In general, this means that a large part of seismic energy remains close to the surface, where it can be measured, which makes even relatively small events observable globally. Also, several seismic measurables are specifically sensitive to parameters related to habitability. This is a promising perspective for the short term installations of seismometers on icy moons.

7 Acknowledgements

The authors acknowledge computational support in the project pr63qo ”3D wave propagation and rupture: forward and inverse problem” at Leibniz-Rechenzentrum Garching. SCS was supported by grant SI1538/4-1 of Deutsche Forschungsgemeinschaft DFG. MvD was supported by grants from the Swiss National Science Foundation (SNF-ANR project 157133 ”Seismology on Mars”) and the Swiss National Supercomputing Center (CSCS) under project ID sm682.

This work was partially supported by strategic research and technology funds from the Jet Propulsion Laboratory, Caltech, and by the Icy Worlds node of NASA’s Astrobiology Institute (13-13NAI7_2-0024). A part of the research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.

The Instaseis wavefield databases available at http://instaseis.ethz.ch. The scripts to reproduce the figures are available on the seismo-live website https://www.seismo-live.org.

All rights reserved prior to publication.

Appendix A Determine the backazimuth of an event from a Europan seismogram

As described in section 3.3, the radial and transverse component of ocean world seismograms differ strongly, which can be used to estimate the backazimuth of an earthquake (i.e. the direction of the earthquake as seen from the receiver). We tested a simple automated implementation of two methods.

Maximize the energy of the Longitudinal wave. The Longitudinal wave is purely longitudinally polarized. Therefore, the horizontal components of the seismogram can be rotated such that energy in the time window around the longitudinal wave is maximal one one horizontal component (R) and minimal on the one perpendicular to it (T). Note that in the presence of noise, the energy in the T-component will not be zero. Also, this method has a 180 degree ambiguity, since the polarity of the longitudinal wave can be positive or negative. 2. 2.

Maximize coherence between vertical and radial component for body waves. The first arriving body waves are compressional waves with a high incidence angle, which have motion purely along the direction of propagation (which is a linear combination of the Z-direction and the R-direction pointing away from the event). The direction away from the event is the one, which has maximal correlation to the vertical component. This method has a much lower resolution than the first one, but has no 180 degree ambiguity and so can used to remove that ambiguity from the higher reolution method 1.

Figure 14 shows the application of both methods to a M3 Europaquake in 83 degree great circle arc distance with a true backazimuth of 76 degree for a Europa model with 5 km ice thickness. It shows that both methods together correctly identify the backazimuth within a few degree.

Appendix B Properties of the numerical simulations

The numerical simulations presented in this manuscript were done using the open source spectral-element wavefield solver AxiSEM. Some modifications were done on the mesher, to make it more stable around very strong velocity contrasts close to the surface. These modifications were added to the main development branch of the software and are included in the recent version 1.4. The simulations were run on the Linux HPC systems SuperMUC of the Leibniz-Rechenzentrum Garching and Piz Daint of the Swiss National Supercomputing Centre CSCS.

Due to the low velocity in the ocean, more (small) elements are needed and simulation is about 50% more expensive compared to a terrestrial planet of the same radius. A 1 Hz simulation for Ganymede or Titan is therefore slightly more expensive than one of the Mars runs done for the InSight blindtest (Ceylan et al., 2017; Clinton et al., 2017).

Appendix C Retrieving seismic waveforms for ocean worlds

While the forward simulations of the seismic wavefield require thousands of CPU-hours; the Instaseis-database-approach allows to reuse the stored wavefield to calculate seismograms for arbitrary source-receiver locations on the fly (van Driel et al., 2015). The disk space requirement of the databases depends on the maximum depth of an event, but is generally in the range of several ten Gigabyte. Databases for the models shown in table 7 are stored on a server at ETH Zürich. Similarly to the terrestrial databases stored at the IRIS service Syngine (Krischer et al., 2017), they can be openly accessed via the Python package Instaseis (http://www.instaseis.net). The following example shows the retrieval of the waveform for a 5.2, pure $M_{rr}$ event in 40 degree distance on Ganymede, including storage into a MiniSEED file.

import instaseis db_path = ’http://instaseis.ethz.ch/icy_ocean_worlds/Gan126km-00pMS-hQ_hyd30km_2s’ db = instaseis.open_db(db_path) src = instaseis.Source(latitude=0.0, longitude=0.0, m_rr=1e17) rec = instaseis.Receiver(latitude=40, longitude=00.0) st = db.get_seismograms(src, rec) st.write(’Ganymede_event.mseed’, format=’MSEED’)