Numerical simulations of marine hydrothermal plumes for Europa and other icy worlds

Icarus 221 (2012) 970–983

Contents lists available at SciVerse ScienceDirect

Icarus journal homepage: www.elsevier.com/locate/icarus

Numerical simulations of marine hydrothermal plumes for Europa and other icy worlds Jason C. Goodman ⇑, Erik Lenferink 1 Wheaton College, Department of Physics and Astronomy, 26 E. Main Street, Norton, MA 02766, USA

a r t i c l e

i n f o

Article history: Received 17 June 2012 Revised 19 August 2012 Accepted 22 August 2012 Available online 12 September 2012 Keywords: Europa Astrobiology Satellites, Dynamics

a b s t r a c t The liquid water interiors of Europa and other icy moons of the outer Solar System are likely to be driven by geothermal heating from the sea floor, leading to the development of buoyant hydrothermal plumes. These plumes potentially control icy surface geomorphology, and are of interest to astrobiologists. We have performed a series of simulations of these plumes using the MIT GCM ocean circulation model. We assume here that Europa’s ocean is deep (of order 100 km) and unstratified, and that plume buoyancy is controlled by temperature, not composition. Our experiments explore a limited region of parameter space, with ocean depth H ranging from 50 to 100 km deep, source heat flux Q between 0.1 and 10 GW, and Coriolis parameter f corresponding to Europa latitudes between 9° and 47°. As predicted by earlier work, the plumes in our simulations form narrow cylindrical chimneys (a few km across) under the influence of the Coriolis effect. These plumes broaden over time until they become baroclinically unstable, breaking up into cone-shaped eddies when they become 10–35 km in diameter; the shed eddies are of a similar size. Large-scale currents in the region of the plume range between 1 and 5 cm/s; temperature anomalies in the plume far from the seafloor are tiny, varying between 10 and 180 lK. Variations in plume size, shape, speed, and temperature are in excellent agreement with previous laboratory tank experiments, and in rough agreement with theoretical predictions. Plume dynamics and geometry are controlled by a ‘‘natural Rossby number’’ which depends strongly on depth H and Coriolis parameter f, but only weakly on source heat flux Q. However, some specific theoretical predictions are not borne out by these simulations: this may occur because the plumes are ‘‘reingesting’’ their own emissions, a process not considered in our earlier theory. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction One of the major discoveries of the Galileo mission to the Jovian system was the discovery of a thick liquid water layer beneath the surface of Europa and other icy moons. The most clear-cut evidence for this comes from measurements of Europa’s induced magnetic field (Kivelson et al., 2000), which requires the existence of a highly conductive material within a few km of the surface—salty water being the only likely cause. Corroborating evidence comes by combining inferences from gravity data (Anderson et al., 1998), which require an ice + water layer 80–170 km thick, with estimates of ice shell thickness from thermal modeling (Hussmann et al., 2002; Spohn and Schubert, 2003), crater shape (Schenk, 2002; Turtle and Pierazo, 2001), and other elements of surface morphology (Schenk and Pappalardo, 2004; Greenberg et al., 2002), which suggest an ice layer thickness between a few km ⇑ Corresponding author. E-mail address: [email protected] (J.C. Goodman). Present address: Northwestern University, Department of Physics and Astronomy, 2145 Sheridan Road, Evanston, IL 60208, USA. 1

0019-1035/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.icarus.2012.08.027

and a few tens of km. This suggests the liquid water ocean may be 70–160 km thick. Oceans like Europa’s may be common throughout the outer Solar System. Comparable magnetic field data have been observed at Ganymede (Kivelson et al., 2002) and Callisto (Kivelson et al., 2000). A liquid layer sandwiched between ice phases has been proposed for Titan based on thermal modeling (Tobie et al., 2005) and measurements of asynchronous rotation (Stiles et al., 2008, 2010); hot spots (Spencer et al., 2006) and an unusual shape (Collins and Goodman, 2007) point to a south polar ‘‘sea’’ within Enceladus Spencer et al. (2009). Using these observations as a starting point, exciting questions are now being asked in the new field of ‘‘planetary oceanography’’ (Vance and Goodman, 2009). What is the composition of the oceans of icy worlds? (Kargel et al., 2000) How do they interact with the seafloor? (Vance et al., 2007) How do they interact with the ice ceiling? (Collins and Nimmo, 2009; O’Brien et al., 2002) What sort of currents are driven by heat sources? (Goodman et al., 2004) What sort of currents are driven by tidal action? (Tyler, 2008) Does the global ocean have a time-mean general circulation? (Vance and Goodman, 2009).

J.C. Goodman, E. Lenferink / Icarus 221 (2012) 970–983

In this paper, our primary focus is on small-scale thermally-driven circulation. Europa and its siblings are heated by a mixture of radiogenic and tidal heating (Moore and Hussmann, 2009). In many thermal models, a large amount of heat is generated directly within the ice shell (Hussmann et al., 2002), but there is always a significant contribution from the silicate interior. Travis et al. (2012) suggest that internal radiogenic heating active over Europa’s entire history is enough to maintain a liquid layer on its own, but has been supplemented by tidal heating in the past 0.5 billion years. To allow thermal equilibrium, this heat must be transferred through the ocean, into the ice shell, and out to space. Vance et al. (2007) describes how hydrothermal flows within the silicate interior will organize and collect the seafloor heat flux into localized vent sources. These same hydrothermal vents have been proposed as drivers of surface geology (Greenberg et al., 1999; O’Brien et al., 2002), providing localized heat sources which may create the small ‘‘lenticulae’’ or ‘‘microchaos’’ seen at the surface. Hydrothermal vents are also of singular interest to astrobiologists (Schulze-Makuch and Irwin, 2002; Chyba and Phillips, 2002; Hand et al., 2007) as sources of reductants for the metabolic processes of hypothetical Europan organisms. However, any interaction between the ice crust and the sea floor must take place through the liquid water layer, and so the liquid’s flow dynamics are crucial to discussion of geology or astrobiology. The first attempt to describe the buoyant rise of hydrothermal plumes in Europa’s liquid interior were made by Collins et al. (2000) and Thomson and Delaney (2001). Collins et al. considered a plume rising into an unstratified, nonrotating ocean; Thomson and Delaney observed that because Europa rotates rapidly, Coriolis forces will confine the plume and control its width, flow, and temperature. Their predictions assumed Europa’s ocean was stratified and rapidly rotating, like Earth’s. These were followed by Goodman et al. (2004), a fairly detailed study based on dimensional analysis and small-scale laboratory experiment. Unlike Thomson and Delaney, Goodman et al. assumed Europa’s ocean is unstratified (see below). This work forms the foundation of our present study. 1.1. Previous work Following techniques commonly used in engineering fluid mechanics and terrestrial oceanography, Goodman et al. (2004) used dimensional analysis to predict the size, velocity, and temperature scales of buoyant hydrothermal plumes arising from a point source of heat at Europa’s seafloor. Their key assumption, repeated here, was that Europa’s ocean was typically unstratified, that is, that the ocean’s buoyancy and (potential) density are uniform with height. This assumption was justified by arguing that if Europa’s ocean is heated at the bottom, buoyancy is generated at the seafloor. In the steady state, the heat and buoyancy must be carried to the ice crust somehow: if this happens convectively, as is typical in geophysical fluids of all types, then convective mixing will act to homogenize the water column and create an unstratified environment. In the work described here, we continue to assume Europa’s ocean is unstratified; in an ongoing parallel study we are testing the limits of this assumption. However, for the moment its worth noting that several possible alternative models produce stratification only at the very top of the ocean. For example, if the ocean is fresh, the ‘‘Melosh et al. stratosphere’’ (Melosh et al., 2004) model may apply, in which the negative thermal expansion coefficient of fresh water at low temperatures and salinities leads to a shallow stratified layer. If the ocean is salty, melting the overlying ice

971

would produce buoyant freshwater, again near the top. But in both cases, the bulk of the water column would remain unstratified, and our analysis would still be relevant there. Goodman et al. (2004) found that all important plume parameters scaled with a single nondimensional number, the ‘‘Natural Rossby Number’’:

Ro ¼

ðBf

3 1=4

Þ

H

In this equation, H is the ocean depth; f is the Coriolis parameter caused by planetary rotation (f = 2Xsin(h), where X is the planetary rotation rate in rad/s, and h is the latitude); and B represents the buoyancy flux emitted by the seafloor source.

B¼

ga Q qw C pw

where Q is the heat flux emitted by the source, g is the local gravity, a is the thermal expansion coefficient, qw is the density of water and Cpw its heat capacity. Like the traditional Rossby number Ro = U/(fL), the natural Rossby number compares the importance of fluid momentum transport to Coriolis forces, and is small (implying geostrophic balance) when the flow is weak and planetary rotation is strong. However, it depends only on fixed external parameters, whereas the traditional Rossby number involves the flow velocity, which is unknown ab initio. In this respect, it is analogous to the ‘‘flux Rayleigh number’’ (Carrigan, 1987), which differs from the traditional Rayleigh number in that it depends on externally-imposed heat fluxes rather than unknown internal temperature gradients. Goodman et al. (2004) found that the plume diameter, temperature, current speed, and development timescale can all be written in terms of Ro⁄. One crucial result is that since Ro⁄ is only weakly dependent on B (to the one-quarter power), the shapes and sizes of hydrothermal plumes depend mainly on water depth and planetary rotation rate, which are fairly well known, and depend very little on the seafloor vent heat flux Q, which is poorly constrained. For instance, pffiffiffiffiffiffiffiffithe theory predicts that plume diameter is proportional to Ro , so that increasing the heat flux by a factor of 10 will increase the diameter of the plume by only 33%. Goodman et al. (2004) used a simple laboratory experiment to validate these theoretical predictions, and to find unknown scaling constants. A trickle of dyed dense, salty fluid was injected into a large tank of water mounted on a rotating turntable. The dense plume sank, mixing turbulently with its surroundings, just as a warm buoyant plume would rise: the physics is the same, just flipped upside down. The observations loosely bore out the theoretical scaling laws. However, the Goodman et al. (2004) study had a number of shortcomings:  Can the laboratory results from cm-scale plumes be directly applied to 100-km structures on a different world? Rotating tank experiments of this type have a long and distinguished pedigree (Fultz, 1952), but those unfamiliar with the technique may worry about its applicability, and those familiar with it may worry about scale-dependent behaviors invalidating the result.  Goodman et al. (2004) was primarily a theoretical study, and performed only a bare minimum of experiments to ‘‘sanity check’’ the theory, 8 in all. These varied only water depth and rotation rate: the plume buoyancy source B was not varied.  It is rather difficult to measure fluid properties and velocity in situ in a tank experiment without disrupting the flow: Goodman et al. (2004) only measured plume size and formation time, plus a rough estimate of fluid velocity.

972

J.C. Goodman, E. Lenferink / Icarus 221 (2012) 970–983

 Tank experiments are limited in the physical situations they can model. We cannot easily mimic the complicated equation of state of a 100-km-deep high-pressure sulfate brine ocean in the lab. Also, in a real planet or Moon the planetary rotation vector is at an angle to the gravity vector, while in a tank experiment they must be parallel: see Section 5.4 for more information on this. To address these issues, and to enable future exploration of this unique geophysical system, we have conducted a series of numerical simulations of hydrothermal plumes using the MIT GCM ocean circulation model.

2. The numerical model The MIT GCM (Adcroft et al., 2004) is a general circulation model designed to efficiently simulate geophysical fluids using a highperformance computing cluster. Originally designed for terrestrial oceanography, the model has evolved into a general-purpose atmosphere/ocean simulator. The MIT GCM is unique in its ability to handle non-hydrostatic 3-dimensional convection in an efficient way. Most ocean models ignore vertical acceleration, and assume water pressure depends only on the weight of the overlying fluid column. This is appropriate for large ocean basins, but not for the highly turbulent, rapid vertical motions we expect in our hydrothermal plume problem. Unlike traditional ocean models, the MIT GCM is able to handle these geophysical turbulence problems (Jones and Marshall, 1997), and it is able to do so faster and more efficiently than general-purpose fluid mechanics solvers. The MIT GCM does have one limitation for our present use: because of its roots in terrestrial oceanography, it expects to be heated and cooled at the top rather than the bottom. It also assumes a free-slip upper boundary condition: for an ocean overlain by an ice sheet, a no-slip upper boundary condition would be more appropriate. Rather than modifying the model to be heated at the bottom, we have simply repeated the tactic used in our earlier experimental work: instead of heating from below, we cool from the top. Because we are using a linear equation of state in which buoyancy depends linearly on temperature and is independent of pressure, the fluid buoyancy is the same with the sign reversed, so the fluid flow is identical, but flipped upside down. In Appendix A, we demonstrate that this inversion has no significant side effects for the present study: with the exception of that section, from this point forward our graphics and word choices are inverted to describe a buoyant, rising plume as would occur in the real ocean. For our hydrothermal plume simulations, we set up a square model domain Lx = Ly = 180 km across, with the plume heat source at the center of the domain. Periodic boundary conditions are used: since in all cases the flow is essentially zero outside the central 1/3 of the domain, the impact of the boundary conditions is negligible. The horizontal resolution of our model is not constant: it varies from a minimum of 500 m  500 m at the center of the domain to a maximum of 1500 m  1500 m at the outer corners, for a total of Nx = Ny = 180 gridpoints. This allows us to adequately resolve turbulent structures in the plume core while minimizing computing effort spent on insignificant flows far from the plume source. The vertical resolution also varies. For our experiments with a water depth H = 100 km, the vertical resolution is 500 m at the sea floor near the vent source, up to 1000 m at the ice/water interface, for a total of Nz = 133 gridpoints. For our experiments with a water depth H = 50 km, we only use the lower Nz = 77 of these. Our model uses a simple linear equation of state: salinity is taken to be a constant, and temperature variations cause a linear var-

iation of density with a thermal expansion coefficient of 3  104 K1. This is a fairly typical value for terrestrial seawater: if the Ro⁄ scaling offered by Goodman et al. (2004) is correct, model behavior should be insensitive to the choice of this value, so long as it is positive and of the right order of magnitude. We are currently developing code to describe more realistic equations of state for sulfate brines (Goodman, 2010); these will be included in future work. Water is slightly compressible: in a real ocean, adiabatic heating upon compression warms or cools fluid parcels moving vertically, just like the much larger adiabatic gradients seen in planetary tropospheres. The MIT GCM uses ‘‘potential temperature’’ (Pedlosky, 1987) (denoted h) as the fundamental thermodynamic variable rather than temperature. Potential temperature corrects for adiabatic compression, so that a parcel moving adiabatically from one depth to another maintains a constant h. For the rest of this article, when we refer to ‘‘temperature’’ we mean potential temperature. Other key model parameters are set as follows: Following timehonored oceanographic modeling practice (Boussinesq, 1877), the tracer diffusion and viscosity are not meant to represent molecular diffusion or viscosity, but rather the mixing effect of eddies smaller than the model grid scale. These values are empirically set to the minimums needed to ensure numerical stability in the model. Many processes in terrestrial ocean modeling are sensitive to the model resolution chosen, and the value of eddy viscosity and diffusivity required. To ensure that our results are not resolutiondependent, we performed simulations with identical physical parameters, but with twice the resolution in each dimension (and half the timestep and half the viscosity). These simulations take 16 times as long to complete, however, so only a few were performed: we cross-compare them with the standard-resolution runs in Section 5.1.1. 2.1. Model forcing Flow in the domain is driven by a constant seafloor heat flow Q. 50% of this heat is applied to a single model gridpoint at the horizontal center of the models sea floor; the remainder is applied to the 4 neighboring gridpoints to provide a less abrupt transition from unheated to heated zones. A tiny spatially-random heating pattern, 1/1000th of the total, is applied near the source to break the model’s artificial perfect symmetry, and allow turbulent motion to begin. Note that our horizontal resolution (500 m) is comparable to a terrestrial hydrothermal vent system: we cannot resolve the decimeter scales of individual source vents. We choose values for Q appropriate for entire hydrothermal systems rather than individual vents—on the order of 1–10 GW (Thomson et al., 1992; Lowell and DuBose, 2005). Vance and Goodman (2009) suggests that the percolating flow mechanics in a hydrothermal system on Europa may be somewhat lower (up to 0.1 GW), so we consider a few values in this range as well. We consider only the thermal forcing provided by the hydrothermal vent system, and neglect the mass flux, momentum, and kinetic energy of the water exiting the vent. This is entirely justified: if such a jet is 200 °C warmer than its surroundings and exits at a speed of 1 m/s (typical for terrestrial systems), the momentum gained through the buoyancy of the hot water will begin to exceed the fluid’s initial momentum just 20 s after it enters the ocean. On the timescales of interest to us (hours to days), thermal buoyancy is totally dominant. An implicit assumption in this work is that while mixing and entrainment is certain to occur on centimeter to meter scales near the vent source, this mixing does not control the shape and size of the final plume. We assume that the hot water exiting the vent is

J.C. Goodman, E. Lenferink / Icarus 221 (2012) 970–983

973

uniformly mixed with ambient fluid within a single 500  500  500 m gridbox in the GCM, and that mixing and entrainment on larger scales control the large-scale plume behavior. This appears to be the case: resolved large-scale mixing dilutes the fluid in the source volume by several orders of magnitude as it ascends to the ice/water interface. 2.2. Time domain Initially, the water in our model has uniform potential temperature (h = 0 °C everywhere). At t = 0, we switch on the heat source, and study the time evolution of the buoyant plume thereafter. We end the experiment once the plume has developed into a well-defined stable structure, and has begun to shed large eddies through the process of baroclinic instability. The 2004 paper expected that these eddies would be ejected from the plume source region into the ocean at large: in our model the separation between ejected eddies and the plume core is a little less clear. However, we can clearly identify the point at which the plume first splits into two separate vortices. This typically occurs between 2 and 10 million seconds into the simulation (20–100 Earth days; 6–30 Europa days). We end the simulation once this transition to baroclinic turbulence has become well developed. The initial spin-up of the plume system is not of great realworld significance: we are more interested in the steady-state delivery of heat to the ice and surrounding ocean. However, in a domain of finite size, we cannot run the simulation forever: the plume fluid will expand to fill the domain. We end the experiment when the plume occupies no more than 1/2 of the domain, to minimize the plume’s interaction with itself across the periodic boundary of the domain.

Fig. 1. Ocean depth, Coriolis parameter, and heat flux for simulations carried out in this work.



2.3. Parameter space In this study, we consider two different ocean depths, h = 50 km and 100 km, in rough agreement with the ranges of H2O layer thicknesses calculated from gravity data by Anderson et al. (1998), and the upper range of ice thicknesses calculated by heat transfer calculations (Hussmann et al., 2002). We consider a range of heat outputs from individual hydrothermal plumes in the range Q = 0.1–10 GW—this covers most of the range considered by Goodman et al. (2004), but with 100 GW plumes ruled out in light of the hydrothermal system modeling results of Lowell and DuBose (2005) and Vance and Goodman (2009). These authors suggest that even lower heat outputs of 1–10 MW might be possible, but the present model lacks the resolution to investigate such small plume structures, so we shall rely on scaling to investigate these microplume systems. The Coriolis parameter f depends on the planetary rotation rate and on latitude. We consider a range of values corresponding to latitudes from 9° to 47° on Europa’s surface. Fig. 1 shows the values chosen in each model run considered in this paper.







2.4. Processes not considered There are a wide range of possibly important processes which are beyond the scope of the present simulation, and it is important to recognize them at the outset. Many of these issues are discussed in more detail in Sections 5.3 and 5.4.  Plume fluid does not cool in contact with the overlying ice layer: the boundary condition there is adiabatic. This can be justified by assuming either that the rate of heat transfer through the ice is slow compared to the motion and evolution of the



plume, or by recognizing that as soon as warm water melts a small amount of ice, highly buoyant freshwater will be formed: this fresh water forms a boundary layer which impedes ocean ? ice heat transfer. We recognize that this process may be important, but a detailed treatment must be postponed for future investigation. We assume a linearized equation of state, in which fluid density is proportional to potential temperature, with a constant of proportionality appropriate for terrestrial seawater. While writing this article, we have been working to compute a polynomial approximation to the full equation of state for a magnesium sulfate brine, which may be more similar to Europa’s true composition. (Vance, S., personal communication, 2010) This detailed equation of state will be included in future simulations. In traditional oceanography and meteorology, only the compo! nent of the planetary rotation vector X which is aligned with the gravity vector ~ g is retained; a scaling analysis is used to drop ! the horizontal component of X (Pedlosky, 1987). This scaling analysis assumes that the aspect ratio of the fluid flow (depth to length) is small, which is true for terrestrial flows, but is definitely not true for the narrow plumes in deep water we are interested in here. The MIT GCM does allow one to include ! the horizontal component of X , and we will do so in future work. See Section 5.4 for details. We assume the plumes are rising into a stagnant ocean, with no large-scale currents. In reality, Europa’s ocean will certainly have tides Tyler (2008), and may have time-average flows as well. We discuss the likely impact of large-scale motions in Section 5.3. As mentioned above, the range of heat fluxes we consider is narrower than one would like, and the model resolution is coarser. These constraints are due to our limited computing resources. However, our resolution is sufficient to permit intense turbulent motion across a spectrum of length scales, and the range of heat fluxes is broad enough to quantitatively describe the effect this parameter has on the flow. These issues are also discussed further in Section 5.3 Limited computing resources require another sacrifice: we do not attempt to simulate the long-term behavior of a hydrothermal plume over years or centuries. Thus, since the plume will continue to grow over the long term, our estimate of plume diameter should be considered a minimum. As for other measured quantities, we discuss their likely long-term trends in Section 5.3, and argue that our results are likely to be robust.

974

J.C. Goodman, E. Lenferink / Icarus 221 (2012) 970–983

3. Data analysis From the model data, we extract several dimensional quantities to compare against theoretical predictions. Several of these quantities depend on the shape of the outer boundary of the plume: we take this boundary to be the 1.0 lK isosurface (the surface on which ambient water has been warmed by 1 millionth of a degree). The boundary is very sharp: choosing 0.1 or 10 lK as the boundary criterion makes little difference to the results.  Time to Baroclinic Instability. This measure of the timescale of plume evolution is estimated by eye using the horizontal temperature and velocity fields. We do not rigorously demonstrate that the various criteria for baroclinic instability are met: what we are actually measuring is the time at which large eddies begin to form. We report our best guess estimate of this time, along with error bars indicating an upper and a lower limit. The lower limit on this time is taken when the plume boundary has clearly departed from a subcircular shape; the upper limit occurs when the plume core has divided itself into multiple independent vortices. We have explored various quantitative ways of estimating the non-circularity of the plume cross-section, but a visual estimate seems to work best.  Baroclinic Cone Width. This is the largest stable size a plume can reach. The 1.0 lK surface defines the border of the plume. Since this surface is not quite a circular cylinder, we compute an ‘‘effective diameter’’ by finding the average cross-sectional area of the plume over the upper 1/4 of its total height; we then find the diameter of a circle which would have the same area.  Flow velocity. We find the peak horizontal velocity at each depth in the middle 1/3 of the plume; we then average the peak speeds at each of these depths.  Temperature anomaly. We find the maximum temperature anomaly (difference between plume temperature and ambient temperature) at mid-depth within the plume. Our measurement schemes generally avoid the base of the plume, due to the rapidly fluctuating high temperatures and velocities found in the near-field turbulence close to the heat source.

4. Results 4.1. Typical plume behavior We begin by discussing the typical behavior of a Europa hydrothermal plume. This typically proceeds in four stages, as discussed by Goodman et al. (2004). In the first stage, buoyancy forces drive plume ascent, and Coriolis forces are not yet significant. But very quickly (within a fraction of a Europa day), Coriolis forces begin to constrain the horizontal expansion of the plume: confined within a ‘‘Taylor column’’ (Pedlosky, 1987), the plume fluid rises vertically. In this second stage, the thermal wind balance sets up counter-rotating horizontal flow, in the direction of planetary rotation (‘‘cyclonic’’) near the seafloor, opposite planetary rotation (‘‘anticyclonic’’) at the top of the rising plume. In the third stage, the top of the rising cylindrical plume encounters the ice/water interface at the top of the ocean. Since the volume of warm plume fluid is increasing, the plume is forced to expand and increase in diameter. However, there is a maximum possible diameter: once the plume grows beyond a limiting diameter given by the Rossby radius of deformation (Pedlosky, 1987), its perimeter becomes increasingly tortuous, until independent eddies begin to split off from the main plume. The final steady state exhibits a balance between warm plume fluid created through turbulent mixing between hydrothermal fluids and surrounding seawater, and export

of plume fluid via eddy formation and migration. Figs. 3 and 2 show two example model runs. We strongly encourage the reader to view the Supplementary online material, which shows video animations of these simulations. 4.2. Nondimensionalization We report all our results using nondimensional quantities, using the fundamental length and inverse time scales H (ocean depth) and f (Coriolis parameter). So, for instance, plume diameters are reported as D/H (the plume diameter as a fraction of ocean depth), the onset time for baroclinic instability is reported as tbcf, etc. Goodman et al. (2004) predicted that all important nondimensional plume parameters should depend only on Ro⁄, the natural Rossby number. Ro⁄ itself can be viewed as a nondimensionalized heat flux. Stronger seafloor heating means a larger Ro⁄. Using nondimensional quantities allows us to plot the multivariate set of model runs (see Fig. 1) on a 2-dimensional graph, more clearly illustrate parameter dependence, and compare results from a wide range of space- and timescales on a single graph. As a demonstration of this, where possible we plot both measurements taken from the numerical hydrothermal plume studies in the present work and also measurements from laboratory tank experiments performed in 2004. Despite the vast difference in spatial scale (cm vs. km) and methodology (water tanks vs. numerical simulation), all data can be compared on the same nondimensional graphs. To help the reader understand the rationale for these nondimensional quantities, we include Table 2. The verbal descriptions are intended to be illustrative rather than precise: in particular, 1/f is proportional to the rotation period, but not equal to it. 4.3. Parameter dependence of plume properties 4.3.1. Diameter As predicted by the 2004 theory, the maximum diameter of the plume at the time it begins to break up due to baroclinic instability increases with Ro⁄ (Fig. 4). That is to say, plumes driven by hotter seafloor sources, or at lower latitude are larger in diameter. The 2004 theory pffiffiffiffiffiffiffiffi predicted that plume diameter should be proportional to Ro (the slope of the solid green line in Fig. 4). This relationship is consistent with our experiments both in the present work and in 2004—the best-fit power law for the 2011 experimental data has an exponent between 0.4 and 0.8. The large scatter in our diameter measurements suggests we may be underestimating their uncertainty. However, we find no apparent relationship between experimental parameters and plume diameter which is independent of Ro⁄. In particular, note that model runs with an ocean depth H = 50 km fall along the same curve as runs with H = 100 km. In dimensional units, the plume diameter ranges from 5 km (for a 0.1 GW plume in a 50 km deep ocean with f = 3  105 s1) to 35 km (for a 10 GW plume in a 100 km deep ocean with f = 1.3  105 s1). Remember that this is the diameter of the plume when it begins to break up; once a steady-state balance is achieved, the complex of plume and surrounding eddies may be somewhat bigger. 4.3.2. Temperature As predicted by the 2004 theory (and by common sense), plumes driven by a stronger seafloor heat source (large Ro⁄) have higher temperatures at mid-depth (Fig. 5). But note that even for plumes driven by the same heat flux, varying Ro⁄ by changing ocean depth or rotation rate affects plume temperature in a consistent way. The fact that plume temperature obeys a smooth power law dependence on Ro⁄ regardless of which external parameter is

J.C. Goodman, E. Lenferink / Icarus 221 (2012) 970–983

975

Fig. 2. Evolution of a typical hydrothermal plume in our simulations. The red surface marks the outer boundary of the plume, where the water temperature is 106 K above ambient. Initially, Coriolis forces constrain the warm buoyant plume fluid to a narrow cylindrical structure. Once this structure reaches the ice/water interface, it begins to spread out laterally, forming a broader cylinder or cone. When the structure broadens beyond the critical Rossby radius of deformation, parts of it begin to break off, forming baroclinic eddies—this process has just begun in the third panel, and is well underway in the fourth. In this simulation, H = 50 km, Q = 1 GW, f = 1.3  105 s1. See supplementary online material for a video animation of this experiment. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

changed demonstrates its status as the core nondimensional variable in the system. However, our 2004 theory predicted that nondimensionalized plume temperature should depend linearly on Ro⁄, but our experiments suggest a much stronger power-law dependence with an exponent of about 2.3. We will explore this discrepancy in Section 5.1. In dimensional units, temperatures range from 105 K above ambient (for a 0.1 GW plume in a 50 km deep ocean with f = 3  105 s1) to 2  104 K above ambient (for a 4 GW plume in 50 km deep ocean with f = 1.95  105 s1).

4.3.3. Horizontal velocity Since horizontal velocities in a geostrophic flow are determined by temperature gradients and plume diameter using the thermal wind equation, our 2004 theory that the horizontal flow ffi pffiffiffiffiffiffiffipredicts velocities should depend on Ro . Fig. 6 demonstrates that while plume swirl velocity does indeed depend strictly on Ro⁄, it does so with a power law exponent between 1.1 and 1.5 – again, significantly greater than predicted by our theory.

Plume velocities vary in our experiments between 0.9 cm/s (for a 0.1 GW plume in a 50 km deep ocean with f = 3  105 s1) to 5 cm/s (for a 10 GW plume in a 50 km deep ocean with f = 1.9  105 s1). While these current speeds are quite modest, they are an order of magnitude larger than predicted by our 2004 theory. This may have implications for their detectability (see Section 6 below).

4.3.4. Time to instability The final parameter predicted by our theory is the time until baroclinic instability begins (Fig. 7). We see once again that this parameter has a strong dependence on Ro⁄, but in this case the dependence is weaker in our experiments than predicted by theory (1.2 to 0.9, as compared to the 2 predicted by theory). In dimensional terms, the time until the onset of baroclinic instability ranges from 2 to 8 megaseconds (20–80 terrestrial days). In all parameters measured, we find a good agreement between our 2004 laboratory tank experiments and our numerical experiments. This agreement helps to validate the numerical experiments.

976

J.C. Goodman, E. Lenferink / Icarus 221 (2012) 970–983

1.4 Europa days

5 Europa days

17 Europa days

Flow speeds: 5-15 cm/s

Temperature: 0.14 mK

21 Europa days

60

40-

km

Fig. 3. Another plume simulation: in this one, f = 1.3  105 s1 as in Fig. 2, but the ocean is deeper (H = 100 km) and the heat source is more intense (Q = 10 GW). As a result, we see more turbulence, and the transition to baroclinic instability is more complicated. Blue cones in the second and third panels show the direction of flow: cyclonic (counterclockwise) near the seafloor, and anticyclonic (clockwise) at the top of the ocean. See supplementary online material for a video animation of this experiment. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

5. Discussion In this section, we look more closely at the following issues:  Why do our simulations disagree with theoretical predictions by Goodman et al. (2004)?  What do our results imply for the study of Europa geophysics and astrobiology?  What limitations and shortcomings remain in our simulations?  Where will future work on this subject take us?

5.1. Possible explanations for simulation – theory mismatch In our 2004 paper, we were able to find a rough agreement between theory and our tank simulations, but were unable to reach a definite conclusion due to the limitations of the tank data. Here we can draw definite conclusions:  Plume parameters do, in fact, depend strictly on the natural Rossby number, Ro⁄ = (Bf3)0.25H1. We find no obvious dependence on B, f, or H in isolation, and experiments with quite

J.C. Goodman, E. Lenferink / Icarus 221 (2012) 970–983

977

Fig. 4. Variation of nondimensionalized plume diameters at time of baroclinic instability with natural Rossby number. Solid green line: 2004 theory predictions (only slope is predicted). Red circles: laboratory tank experimental data from Goodman et al. (2004). Other symbols: simulation results (diamonds: 100 km ocean depth; squares: 50 km ocean depth; purple triangles: high-resolution runs). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6. Variation of nondimensionalized plume velocities with natural Rossby number. Solid green line: 2004 theory predictions (only slope is predicted). Symbols: simulation results (diamonds: 100 km ocean depth; squares: 50 km ocean depth; purple triangles: high-resolution runs). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. Variation of nondimensionalized plume temperatures with natural Rossby number. Solid green line: 2004 theory predictions (only slope is predicted). Symbols: simulation results (diamonds: 100 km ocean depth; squares: 50 km ocean depth; purple triangles: high-resolution runs). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7. Variation of nondimensionalized time to baroclinic instability (eddy shedding) with natural Rossby number. Solid green line: 2004 theory predictions (only slope is predicted). Red circles: laboratory tank experimental data from Goodman et al. (2004). Other symbols: simulation results (diamonds: 100 km ocean depth; squares: 50 km ocean depth; purple triangles: high-resolution runs). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

different values of these parameters but the same value of Ro⁄ behave almost identically.  However, the specific power-law dependence on Ro⁄ is often at odds with the theoretical predictions from our 2004 paper. Plume diameter and instability time depend more weakly on Ro⁄ than theory predicts; temperature and horizontal velocity depend more strongly. In this section we look for an explanation for this disagreement. We are able to rule out the effects of model resolution or observational bias, but find an unexpected behavior (plume reingestion) which might explain the results. 5.1.1. Model resolution The most obvious problem with any numerical simulation is a lack of spatiotemporal resolution. This is especially true in studies of turbulence, where the real-world flow has structures from centimeter to kilometer scale. In practice, numerical simulations operate by resolving the smallest practical scales, and using a

‘‘turbulent closure scheme’’ (Kundu, 1990) to express the energy and momentum transport of sub-gridscale turbulence. In a successful numerical simulation, the dynamics should be wellresolved enough that how we handle the smallest scales should not control the behavior of the large-scale flow. Because we are trying to simulate long-term behavior of 3dimensional turbulence on a fairly small computing cluster, this study uses a much coarser grid resolution than we would prefer. Is a 500-m grid scale small enough to adequately resolve the turbulent transport processes which set the properties of our hydrothermal plumes? If not, perhaps the disagreement between theory and simulation is simply because the simulation is inadequate. An easy way to check is to do a limited number of simulations at higher resolution. If these behave significantly differently than the standard-resolution runs, we should conclude that our model resolution is inadequate, and our results are called into question. We have performed four simulations with twice the model resolution: a minimum of 250 m at the plume source, up to 750 m at the far corners. Following standard numerical scaling procedures, we also

978

J.C. Goodman, E. Lenferink / Icarus 221 (2012) 970–983

Table 1 Key parameters for MIT GCM hydrothermal plume simulations. See Adcroft et al. (2004) for details. Timestep

153–409 s

Advection scheme Diffusion constant

33 0.06 m2/s

Viscosity 4th-order ‘‘hyperviscosity’’ Thermal expansion coefficient

0.09 m2/s 3  104

Shorter timesteps used for runs with stronger heating 3rd-order DST with flux limiting Required for numerical stability: see comment below Required for numerical stability Required for numerical stability

3  104 K1

Linear equation of state

Table 2 Nondimensionalizations used to compare experiments. Quantity

Nondimensionalization

Rationale

Diameter

lcone H

Time to instability Velocity

Tbcf

Plume diameter compared to ocean depth Time compared to one rotation period Velocity compared to velocity needed to rise through ocean depth in one rotation period Plume buoyancy compared to buoyancy needed to rise ballistically through ocean depth in one rotation period

Temperature

V Hf

Tg a Hf 2

halve the timestep, halve the eddy viscosity and diffusivity, and reduce the hyperviscosity by a factor of 8, as compared with the values given in Table 1. These simulations come at significantly greater computational cost, requiring 24 = 16 times as much computing time – thus, we have kept these simulations to the minimum needed to test the resolution-dependence hypothesis. These high-resolution runs are shown as purple triangles in Figs. 4–7. We see no significant difference between these simulations and the equivalent low-resolution runs—in fact, often the agreement is so good that the symbols overlap and obscure one another. This alone does not prove that the resolution we are using is adequate, but it does demonstrate that our simulation–theory misfit is not a straightforward consequence of inadequate resolution.

5.1.2. Observer measurement bias There are several qualitative ‘‘judgement calls’’ in measuring the plumes in these experiments. At what point, precisely, does the plume go baroclinically unstable? It is a sea of geostrophic turbulence at all times: when can we truly say that true baroclinic eddies have formed? The plume edge is sharp, but not perfectly so. Do our results change when our measurements take the plume margin at 107 K above ambient, rather than 106? As a consistency check, the measurements reported here were performed by both authors separately. Differences between the data reported by each author were generally smaller than the uncertainties caused by the chaotic fluctuation of the plumes themselves. Spot checks were performed varying the plume-edge criterion by 2 orders of magnitude: this changed the measured plume diameter by <20%, but did not change the slope of the data in Fig. 4; other measured values did not change significantly.

5.1.3. Plume reingestion Our 2004 theory, following earlier work by Fernando et al. (1998), assumes that at the seafloor heat source, the hydrothermal fluid is mixing and entraining with cold ambient seawater. However, this is not quite what we see in either numerical simulations or tank experiments. Turbulent motion is strong enough to mix warm water throughout the region near the plume, such that the buoyant plume is entraining a mix of ambient seawater and already-warmed plume fluid pulled down from above—the plume ‘‘reingests’’ its effluent (see Fig. 8). The consequences of this unexpected behavior are not obvious. Certainly, reingestion will make the plume fluid somewhat warmer than it would otherwise be; this will cause greater horizontal swirl velocities and larger baroclinic eddies (since the Rossby radius of deformation will be larger). However, it is not clear whether reingestion will affect the dynamics of large or small plumes most strongly, and thus it is not clear how it affects the scaling of plume parameters with Ro⁄. Plume reingestion is an important unexplored phenomenon. Additional work is needed to assess whether it can explain the theory-vs.-experiment mismatch observed here.

5.2. Implications for Europa geoscience Our results have important implications for the study of Europa’s surface geology and astrobiology, and may be observable with future spacecraft missions.

Fig. 8. Left: entrainment flow assumed by Goodman et al. (2004). Plume source fluid mixes with ambient cold water and rises away from source region. Right: entrainment flow observed in our model. Plume source mixes with both ambient cold water and reingested mixed fluid.

J.C. Goodman, E. Lenferink / Icarus 221 (2012) 970–983

5.2.1. Surface geology Goodman et al. (2004) concluded that the range of diameters of hydrothermal plumes in Europa’s ocean should be 25–50 km, and argued that as a result, the various small subcircular features on Europa’s surface (called ‘pits, domes, and spots’, ‘microchaos’, or ‘lenticulae’ by various authors (Spaun et al., 2002; Collins and Nimmo, 2009; Greenberg et al., 1999)), with sizes down to a few km, cannot be created by basal melting of the ice layer (‘‘meltthrough’’) caused by the broad plumes. Our computer simulations find a similar but somewhat smaller range of diameters (10–35 km) under the same range of natural Rossby number. But for the very smallest plumes we considered (108 W) and the shallower water depth (50 km), we did find plume diameters of as little as 5 km, which are comparable to lenticulae on Europa. Keep in mind, however, that as eddies are shed from the columnar plume and migrate outward, they will spread their heat more widely, so our diameters should be considered lower bounds. Also, the very smallest plumes observed here have very small heat fluxes per unit area, which may be insufficient to cause significant amounts of shell melting. However, our largest plumes are quite comparable in size to mid-sized to large chaos regions. This is not strong evidence in favor of the idea that large chaoses are caused by plumes, but a drastic mismatch in size would be evidence against the idea. 5.2.2. Astrobiology Astrobiologists may be interested in the velocity and efficiency of vertical transport from the sea floor to the ice/water interface. This transport could carry living organisms or nutrients, and have implications for both the metabolism of hypothetical native organisms and the survival of terrestrial species introduced by visiting spacecraft—a key question in decisionmaking regarding planetary protection (Sogin et al., 2012). We hope to work more closely with astrobiologists to evaluate and use the transport data generated by our model. We present the following data (Fig. 9) as an illustration of what’s possible. For our

979

simulation with f = 1.95  105 s1, H = 100 km, Q = 1 GW, we find typical average upward velocities within the core of the plume of 104–102 m/s; this net upward flow is returned by a downward flow far from the plume with a mean velocity on the order of 105 m/s in our model—however, in the real world this value would depend on the spacing between plumes. A lagrangian tracer moving at these average flow velocities would travel from seafloor to ice/water interface on a timescale of years; the return trip to the seafloor would take centuries. However, this is misleading, because vertical velocity is dominated by turbulence. The typical standard deviation of vertical flow velocity is 5  103 to 5  102 m/s within the plume, and 5  104 m/s far away from it. The time variation of the flow is large relative to the mean, so vertical transport of nutrients and organisms will occur via a ‘‘random walk’’ process. The MIT GCM allows one to follow tracer particles within the flow, if more detailed trajectories are needed. 5.2.3. Spacecraft detection of hydrothermal plumes Jupiter’s magnetic field permeates Europa, and Europa’s response to variations in that magnetic field are our strongest evidence for the existence of an ocean (Kivelson et al., 2000). However, if Europa’s ocean moves, it will create an induced magnetic field which could potentially be detected by a spacecraft magnetometer. Could hydrothermal plumes be detected in this way? A full answer to this question is well beyond the scope of this paper, but we can ask whether future work on this subject is justified. Tyler (2011) provided the most complete and thorough description of the effect of fluid flow on Europa’s time-varying magnetic field. He focuses on tidal motions, which he predicts will have a velocity amplitude of a few cm/s: this leads to magnetic field signatures which are very small (of order 1 nT), but nevertheless should be detectable due to their periodic nature. Goodman et al. (2004) predicted very slow ocean horizontal flow velocities of a few mm/s. However, the GCM simulations in

Fig. 9. Vertical velocity statistics for a typical hydrothermal plume simulation (f = 1.95  105 s1, H = 100 km, Q = 1 GW). Left: time-average vertical velocity along the plume centerline (dashed line) and outside the plume, 25 km away (solid line). Right: standard deviation of vertical velocities, inside and outside the plume.

980

J.C. Goodman, E. Lenferink / Icarus 221 (2012) 970–983

this paper find typical horizontal flow velocities an order of magnitude larger than the 2004 predictions, comparable to Tyler’s tidal flow velocity estimate. Does this mean that hydrothermal plumes should also be detectable via magnetometer? Not necessarily. Unlike tidal flows, our hydrothermal plumes are very small in horizontal scale, so their magnetic signatures should decrease rapidly with distance. Their flows are time-varying but not periodic, so ‘‘signal stacking’’ techniques for detecting periodic signals will not work. However, it is worth noting that horizontal flow is only one possible source for magnetic field perturbations. Tyler also discusses the tightening of magnetic field lines caused by vertical stretching of fluid columns. In Tyler’s tidal model, this is driven by tidal heaving of the sea surface, which he argues would have vertical amplitudes of 19 m, and thus velocities on the order of 0.1 mm/s. Our hydrothermal plumes have vertical velocities 10–100 times greater than this, and so would cause much larger vertical compression of fluid columns, with a potentially much stronger magnetic field effect. This is by no means a conclusive argument, but it does suggest that further work be done on this subject. However, a full treatment might require integrating magnetohydrodynamic equations into the MIT GCM, which is no small task. 5.3. Limitations Here, we remind the reader of the limitations and shortcomings of our results. Some of these are reiterated from Section 2.4, while others are more general observations which motivate future research. 5.3.1. Long-term behavior Our simulations are not intended to describe the full life history of a hydrothermal plume. We end the experiment once baroclinic eddy shedding has become well developed, typically after a few months or a year, but a real seafloor heat source might be active for orders of magnitude longer. We do not run our simulations to a steady state: in fact, since there is nowhere for the heat to escape, a steady state is impossible. What elements of our results are likely to change over the long term? The volume of water heated by the plume will clearly increase over time, as the central plume continuously sheds baroclinic eddies. Thus, the plume diameters here should be considered minima when considering the long-term area of influence. Exactly how large the plume-heated area will become depends on the processes which remove heat at the ice interface, a topic of future research. One might expect the plumes to gradually become warmer, but this is not the case. We find in our model runs that the median plume temperature is set quite quickly (within a day) by smallscale turbulence near the heat source; as time goes on, continued heating increases the volume of affected fluid, but not its temperature. Velocity is strictly dependent on the temperature anomaly, and also does not change. To summarize, as a plume continues to develop over the long term, we find it gets larger, but not hotter or faster. Thus, our most important conclusions (that plumes are big, cool, and slowmoving) are robust. 5.3.2. Heat output In our experiments, we used fairly high heat fluxes from our plumes: Vance and Goodman (2009) suggests a range of plume outputs from 102 to 108 W (compared to our 108–1010 W). Our choice of heat flux was motivated by a desire to remain consistent with values used by Goodman et al. (2004), and by computational constraints. (Weaker plumes have smaller length scales and longer time scales, requiring significantly more computing time.) The low end of the Vance and Goodman (2009) range is unlikely to play a

significant role in geology or geochemistry on Europa; the upper end should be reachable with our computer models. One would like to simply apply the Goodman et al. (2004) scaling law to predict the behavior of 1 MW plumes, but per Section 5.1, the theory is not up to the task. While we consider heat flows beyond the Vance and Goodman (2009) range, these may be quite appropriate for groups of multiple plume sources located within a vent field hundreds of meters across.

5.3.3. Large-scale currents and tides Our model’s plumes rise within a stagnant ocean. Large-scale barotropic ocean currents would move the ocean with respect to the plume source, complicating our results. For continuous ocean currents, the effect would be like aiming a can of spray paint at a conveyor belt: a long stripe of ocean will be heated by the plume. For tidal flows, it is more like aiming a spray can at a vibrating plate. The key factor is the amplitude of horizontal tidal displacement: if it is small compared to the width of the plume, tidal motion can be ignored. Tyler (2011) suggests tidal flow velocity amplitudes of 8.4 cm/s, though that estimate does not include dissipative effects. Including dissipation, tidal velocities are likely to lie between 1 and 10 cm/s (Tyler, personal communication, 2012) Dividing by the angular frequency of Europa’s tides (X = 2.05  105 s1), this implies a displacement amplitude of 0.5–5 km: fairly small compared to the plumes’ diameters. Timerectified tidal flow, shear flow, or seafloor tidal turbulence might still disrupt the plumes, but the large-scale depth-independent tidal flow will not carry them far from the heat source that created them. A related question: if Europa’s ocean has time-mean global ocean currents, how large would these currents need to be to disrupt the plume behavior discussed here? Since we have previously assumed an unstratified ocean, we must also assume the putative background flow is independent of depth according to the thermal wind equation Pedlosky (1987). In our model, the seafloor heat source continues to warm the same overlying body of water until the plume becomes baroclinically unstable. If external ocean currents carry the plume downstream more than one plume radius during that time, it will no longer be warmed by the heat source. This suggests that external currents are likely to disrupt the plume if they are larger than (plume radius/baroclinic instability time). For our experimental plumes, this works out to a few mm/s. This may seem incredibly small, but keep in mind that mean velocities in Earth’s abyssal oceans are typically [1 mm/s Vallis (2012), and Europa’s global circulation may be more weakly forced than Earth’s.

5.3.4. Heat flow at ice interface Our model ignores heat transfer from the ocean to the overlying ice. We believe this heat transfer is likely to be slow due to the formation of a buoyant boundary layer of meltwater at the ice/water interface which will inhibit convection: however, this is as yet an untested hypothesis.

5.3.5. Equation of state We prescribe a very simple linear equation of state, with density varying linearly with temperature. Salinity effects are ignored entirely, and the fact that the chemical composition of Europa’s ocean is unlikely to resemble Earth’s has not been considered. These choices were made to avoid making specific claims about the detailed composition of Europa’s ocean, which are not well understood. In ongoing work, we consider specific endmember compositions for Europa’s ocean, with appropriate equations of state (see Section 5.4).

J.C. Goodman, E. Lenferink / Icarus 221 (2012) 970–983

5.3.6. Resolution Our model’s fairly coarse resolution is adequate to resolve some of the turbulent mixing within the plume, and the model passes simple resolution sensitivity tests, but it is quite possible that resolving turbulence on much finer scales may affect the results.

981

5.4. Further research 5.4.1. Realistic equation of state For this work, we used a simple linear equation of state, in which density (and thus buoyancy) depended linearly on temperature. In reality, density is a nonlinear function of temperature, salinity, pressure, and the types of salt ions present. Recently, Vance and Brown (2011) performed measurements of the equation of state of MgSO4 brines and NH3 mixtures, which allow us to consider a full realistic equation of state for icy world oceans. Preliminary descriptions of this equation of state (Vance and Goodman, 2009) suggest that qualitatively new behaviors may arise from the nonlinear behavior. However, the amount and type of salt in Europa’s ocean remains a free parameter in our studies. As the present paper was being completed, we were able to successfully add the experimentally-determined equations of state from Vance and Brown (2011) into the MIT GCM, allowing us to investigate these issues. Results from this analysis will be described in a future paper. 5.4.2. Melting at the ice base and double diffusion We intend to build a parameterization for the complicated melting, conduction, and convection processes that occur at the ice/water boundary. Doing so will inevitably entail a parameterization of ‘‘double diffusion’’ (Schmidt, 1994), a process whereby salt and heat can be transported vertically even though the water column is convectively stable. These mixing processes occur on a much smaller scale than our model can resolve; we intend to use the Zhang et al. (1998) parameterization to deal with this issue.

Fig. 10. A hypothetical structure for Europa’s global ocean circulation, driven by turbulent cascade of small-scale turbulence caused by hydrothermal plumes. Future work will test the validity of this hypothesis.

5.4.3. Horizontal component of planetary rotation In terrestrial oceanography, one typically calculates the Coriolis force using only the component of the planetary angular momentum vector parallel to gravity (Xz). This is easily justified by a scaling analysis (Pedlosky, 1987) when the aspect ratio of the flow

Fig. A.1. Snapshot comparison of plumes heated from below (left) and cooled from above and shown inverted (right). For both plumes, Q = 1 GW, f = 3  105 s1, H = 100 km; simulation time = 4.09  106 s. While the two plumes show different instances of chaotic turbulence, their size, shape, velocity, and temperature structure are the same.

982

J.C. Goodman, E. Lenferink / Icarus 221 (2012) 970–983

(depth to width) is small—typically the case for large-scale terrestrial oceanography. However, this is not the case for our Europan hydrothermal plumes. The importance of the component of planetary rotation perpendicular to gravity (Xy) has not been extensively studied for oceanic convection, with the exception of a remarkable (and remarkably dangerous) experiment by Sheremet (2004), which mounted a large water tank on the end of a rapidly rotating 2.5-m support arm, so that the centripetal acceleration of the tank created a ‘‘gravity’’ component perpendicular to the axis of rotation. The MIT GCM allows one to include the Xy term in the simulation, which should allow us to reproduce and expand on Sheremet’s work, in contexts appropriate to both icy worlds and Earth’s oceans, in comparative safety. 5.4.4. Global ocean circulation Finally, we have also begun to investigate the role of hydrothermal plumes in shaping the global ocean circulation of icy world oceans. Tidal flows have been studied extensively (Tyler, 2008, 2011), but other aspects of the global circulation remain a mystery. In previous work (Vance and Goodman, 2009), we hypothesized that small-scale turbulence created by hydrothermal plumes should coalesce and cascade to larger scales, more strongly in the horizontal than the vertical, leading to the creation of a horizontally banded structure of alternating east–west jets (Fig. 10) quite similar in pattern and dynamics to the atmosphere of Jupiter. We have begun to build a global-scale ocean model for Europa using the MIT GCM to investigate this hypothesis. 6. Conclusion Our experiments broadly confirm the theoretical description of hydrothermal plumes given by Goodman et al. (2004). Plume shapes, sizes, temperatures, velocities, etc. are roughly in accord. However, plume velocities are larger, and velocities and temperatures are more sensitive to the intensity of the hydrothermal source (as measured by Ro⁄) than predicted by theory. Care should be taken in extrapolating these data to hydrothermal plumes much weaker or stronger than what we have modeled. We have ruled out a variety of hypotheses to explain the model/ theory mismatch, but have no satisfactory explanation as yet. Figuring this out is a primary focus for our future research. We also intend to continue investigation into the nonlinear equation of state of icy world ocean brines, the effect of horizontal components of planetary rotation, the impact of hydrothermal plumes on the global ocean circulation, and the detectability of these plumes using magnetometers. 6.1. Planetary oceanography It should be clear that the work done for this paper only begins to scratch the surface of the exciting new field of ‘‘planetary oceanography’’, the study of ocean dynamics throughout our Solar System and beyond. There are enough new challenges and exciting opportunities to keep us busy for years to come. Acknowledgements The authors gratefully acknowledge support from the NASA Astrobiology Institute’s ‘‘Astrobiology of Icy Worlds’’ program at the Jet Propulsion Laboratory. Appendix A. Effect of inverting the plume dynamics This paper attempts to model the behavior of an ocean heated from below. But since the MIT GCM is designed as a terrestrial

ocean model, it was easier to perform our simulations in an ocean that was cooled from above. We argued in Goodman et al. (2004) and Section 2 that for a linear equation of state, this should give the same result, and throughout this paper we have flipped the output from our model upside-down. This was easier to do than adding basal heating to the MIT GCM. For future work, though, we will consider nonlinear, pressuredependent equations of state, where the vertical symmetry is broken. After most of the simulations reported in the main text were completed, we modified the MIT GCM to be heated from below, and were able to confirm that it does not affect the parameters measured here. Fig. A.1 shows an example of this: while the two model runs are different in the details of their turbulence—to be expected in a chaotic system—they have no significant difference in shape or dynamics, for any of the parameters measured in Figs. 4–7. Appendix B. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.icarus.2012. 08.027. References Adcroft, A., Hill, C., Campin, J.M., Marshall, J., Heimbach, P., 2004. Overview of the formulation and numerics of the MIT GCM. In: Proceedings of the ECMWF Seminar series on Numerical Methods, Recent Developments in Numerical Methods for Atmosphere and Ocean Modelling. Anderson, J.D., Schubert, G., Jacobson, R.A., Lau, E.L., Moore, W.B., Sjogren, W.L., 1998. Europa’s differentiated internal structure: Inferences from four Galileo encounters. Science 281, 2019–2022. Boussinesq, J., 1877. Essai sur la théorie des eaux courantes. Mémoires présentés par divers savants à l’Académie des Sciences XXIII. Carrigan, C.R., 1987. The magmatic rayleigh number and time dependent convection in cooling lava lakes. Geophys. Res. Lett. 14, 915–918. Chyba, C.F., Phillips, C.B., 2002. Europa as an abode of life. Origins Life Evol. 32, 47– 68. Collins, G.C., Goodman, J.C., 2007. Enceladus’ south polar sea. Icarus 189, 72–82. Collins, G.C., Nimmo, F., 2009. Chaos on Europa. In: Pappalardo, R.T., McKinnon, W.B., Khurana, K.K. (Eds.), Europa. University of Arizona Press, pp. 259–281. Collins, G.C., Head, J.W., Pappalardo, R.T., Spaun, N.A., 2000. Evaluation of models for the formation of chaotic terrain on Europa. J. Geophys. Res. – Planets 105, 1709– 1716. Fernando, H.J.S., Chen, R., Ayotte, B.A., 1998. Development of a point plume in the presence of background rotation. Phys. Fluids 10, 2369–2383. Fultz, D., 1952. On the possibility of experimental models of the polar-front wave. J. Meteorol. 9, 379–384. Goodman, J.C., 2010. Modeling vertical structure and heat transport within the oceans of ice-covered worlds, P24A-03. In: American Geophysical Union, Fall Meeting; 2010. Goodman, J.C., Collins, G.C., Marshall, J., Pierrehumbert, R.T., 2004. Hydrothermal plume dynamics on Europa: Implications for chaos formation. J. Geophys. Res. – Planets 109, E030008. http://dx.doi.org/10.1029/2003JE002073. Greenberg, R., Hoppa, G.V., Tufts, B.R., Geissler, P., Riley, J., Kadel, S., 1999. Chaos on Europa. Icarus 141, 263–286. Greenberg, R., Geissler, P., Hoppa, G.V., Tufts, B.R., 2002. Tidal-tectonic processes and their implications for the character of Europa’s crust. Rev. Geophys. 40, 1004. http://dx.doi.org/10.1029/2000RG000096. Hand, K.P., Carlson, R.W., Chyba, C.F., 2007. Energy, chemical disequilibrium, and geological constraints on Europa. Astrobiology 7, 1006–1022. Hussmann, H., Spohn, T., Wieczerkowski, K., 2002. Thermal equilibrium states of Europa’s ice shell: Implications for internal ocean thickness and surface heat flow. Icarus 156, 143–151. Jones, H., Marshall, J., 1997. Restratification after deep convection. J. Phys. Oceanogr. 27, 2276–2287. Kargel, J.S., Kaye, J.Z., Head, J.W., Marion, G.M., Sassen, R., Crowley, J.K., Ballesteros, O.P., Grant, S.A., Hogenboom, D.L., 2000. Europa’s crust and ocean: Origin, composition, and the prospects for life. Icarus 148, 226–265. Kivelson, M.G., Khurana, K.K., Russell, C.T., Volwerk, M., Walker, R.J., Zimmer, C., 2000. Galileo magnetometer measurements: A stronger case for a subsurface ocean at Europa. Science 289, 1340–1343. Kivelson, M.G., Khurana, K.K., Volwerk, M., 2002. The permanent and inductive magnetic moments of Ganymede. Icarus 157, 507–522. http://dx.doi.org/ 10.1006/icar.2002.6834. Kundu, P., 1990. Fluid Mechanics. Academic Press. Lowell, R.P., DuBose, M., 2005. Hydrothermal systems on Europa. Geophys. Res. Lett. 32. http://dx.doi.org/10.1029/2005GL022375.

J.C. Goodman, E. Lenferink / Icarus 221 (2012) 970–983 Melosh, H., Ekholm, A., Showman, A., Lorenz, R., 2004. The temperature of Europa’s subsurface water ocean. Icarus 168, 498–502. Moore, W.B., Hussmann, H., 2009. Thermal evolution of Europa’s silicate interior. In: Pappalardo, R.T., McKinnon, W.B., Khurana, K.K. (Eds.), Europa. University of Arizona Press. O’Brien, D.P., Geissler, P., Greenberg, R., 2002. A melt-through model for chaos formation on Europa. Icarus 156, 152–161. Pedlosky, J., 1987. Geophysical Fluid Dynamics. Springer-Verlag. Schenk, P.M., 2002. Thickness constraints on the icy shells of the Galilean satellites from a comparison of crater shapes. Nature 417, 419–421. Schenk, P.M., Pappalardo, R.T., 2004. Topographic variations in chaos on Europa: Implications for diapiric formation. Geophys. Res. Lett. 31, L16703. http:// dx.doi.org/10.1029/2004GL019978. Schmidt, R.W., 1994. Double diffusion in oceanography. Annu. Rev. Fluid Mech. 26, 255–285. Schulze-Makuch, D., Irwin, L.N., 2002. Energy cycling and hypothetical organisms in Europa’s ocean. Astrobiology 2, 105–121. Sheremet, V.A., 2004. Laboratory experiments with tilted convective plumes on a centrifuge: A finite angle between the buoyancy force and the axis of rotation. J. Fluid Mech. 506, 217–244. Sogin, M.L., Collins, G.C., et al., 2012. Assessment of Planetary Protection Requirements for Spacecraft Missions to Icy Solar System Bodies. The National Academies Press. Spaun, N.A., Head, J.W., Pappalardo, R.T., 2002. The spacing distances of chaos and lenticulae on Europa. Lunar Planet. Sci. XXXIII. Abstract No. 1723. Spencer, J.R. et al., 2006. Cassini encounters enceladus: Background and the discovery of a south polar hot spot. Science 311, 1401–1405. Spencer, J. et al., 2009. Enceladus: An active cryovolcanic satellite. In: Dougherty, M., Esposito, L.W., Krimigis, S.M., et al. (Eds.), Saturn from Cassini–Huygens. Springer, pp. 683–724. Spohn, T., Schubert, G., 2003. Oceans in the icy Galilean satellites of Jupiter? Icarus 161, 456–467. http://dx.doi.org/10.1016/S0019-1035(02)00048-9.

983

Stiles, B.W. et al., 2008. Determining Titan’s spin state from Cassini radar images. Astron. J. 135, 1669–1680 (Note erratum (Stiles et al., 2010)). Stiles, B.W. et al., 2010. Erratum: ‘‘Determining Titan’s spin state from Cassini radar images’’. Astron. J. 139, 311. http://dx.doi.org/10.1088/0004-6256/139/1/311. Thomson, R.E., Delaney, J.R., 2001. Evidence for a weakly stratified Europan ocean sustained by seafloor heat flux. J. Geophys. Res. – Planets 106, 12335–12365. Thomson, R.E., Delaney, J.R., McDuff, R.E., Janecky, D.R., McClain, J.S., 1992. Physical characteristics of the Endeavour Ridge hydrothermal plume during July 1988. Earth Planet. Sci. Lett. 111, 141–154. Tobie, G., Grasset, O., Lunine, J.I., Mocquet, A., Sotin, C., 2005. Titan’s internal structure inferred from a coupled thermal–orbital model. Icarus 175, 496–502. Travis, B.J., Palguta, J., Schubert, G., 2012. A whole-Moon thermal history model of Europa: Impact of hydrothermal circulation and salt transport. Icarus 218, 1006–1019. Turtle, E.P., Pierazo, E., 2001. Thickness of a Europan ice shell from impact crater simulations. Science 294, 1326. http://dx.doi.org/10.1126/science.1062492. Tyler, R.H., 2008. Strong ocean tidal flow and heating on moons of the outer planets. Nature 456, 770–772. Tyler,R.H.,2011.MagneticremotesensingofEuropa’soceantides.Icarus211,906–908. Vallis, G.K., 2012. Climate and the Oceans. Princeton University Press. Vance, S., Brown, J., 2011. Exploring deep icy world oceans through new experimental equations of state for aqueous MgSO4 and NH3. In: American Geophysical Union, Fall Meeting; 2011. Abstract #P23D-1735. Vance, S., Goodman, J.C., 2009. Oceanography of an ice-covered Moon. In: Pappalardo, R.T., McKinnon, W.B., Khurana, K.K. (Eds.), Europa. University of Arizona Press. Vance, S., Harnemeijer, J., Kimura, J., Hussmann, H., deMartin, B., Brown, J.M., 2007. Hydrothermal systems in small ocean planets. Astrobiology 7, 987–1005. http://dx.doi.org/10.1089/ast.2007.0075. Zhang, J., Schmidt, R.W., Huang, R.X., 1998. Sensitivity of the GFDL modular ocean model to parameterization of double-diffusive processes. J. Phys. Oceanogr. 28, 589–605.