Bulletin of Glaciological Research ,0 (,**2) ,-ῌ-,ῌJapanese Society of Snow and Ice

Simulation of cryovolcanism on Saturn’s moon Enceladus with the Green-Naghdi
theory of thermoelasticity
Swantje BARGMANN+*, Ralf GREVE, and Paul STEINMANN+**+ Chair of Applied Mechanics, University of Kaiserslautern, P.O. Box -*.3, 010/- Kaiserslautern, Germany, Institute of Low Temperature Science, Hokkaido University, Kita-+3, Nishi-2, Kita-ku, Sapporo *0*ῌ*2+3, Japan(Received April 1, ,**2; Revised manuscript accepted August +*, ,**2)

Abstract
In ,**/, the Cassini spacecraft proved the existence of cryovolcanism, i.e., the icy counterpart of
volcanism on Earth, on Saturn’s moon Enceladus during its close ﬂy-bys. In particular, water-richplume venting was discovered in the south polar region. Thus, Enceladus was found to be one out of
three outer solar bodies to be geologically active. This contribution is concerned with the modellingand computation of this phenomenon. For the underlying thermoelastic description of ice at cryogenic temperatures, we resort to the Green-Naghdi approach. The Green-Naghdi theory includes theclassical Fourier approach, but, in addition to that, it is a lot more general as it also allows for othertypes of heat propagation. The numerical implementation is carried out with the help of the ﬁnite
element method. Results show that lateral spreading of internal and surface warming away from anactive volcanic vent increases strongly with increasing contribution of the non-classical heat ﬂux.
Agreement with available high-resolution surface temperature data based on infrared spectrometry
seems to be best if the non-classical heat ﬂux contributes signiﬁcantly to the total heat transport.
Complementary laboratory studies would be required in order to strengthen this speculative, yetpromising idea.

+.

Introduction

Enceladus is one of Saturn’s inner satellites. It
was discovered in +123 and orbits Saturn in an almost
circular orbit with a semi-major axis of ,-1,3.2 km anda period of +.-1 days. Enceladus has the shape of a
ﬂattened ellipsoid with a mean diameter of /*. km,
which makes it the sixth-largest Saturnian moon. It
consists of a rocky core and an icy mantle and surface,
the mean surface temperature being at a cryogenic
value of 11 K. In many ways, the general properties
of Enceladus are very di#erent from those of othersatellites, and thus, it is a very active research topic.
With an albedo of *.33, it has the most reﬂective
surface of any body in the solar system. Moreover, it
shows recent geological activity, which is very unusual for such a small body.
Of particular interest is Enceladus’ south polar
region. Active cryovolcanism was found there in ,**/
during a ﬂy-by of the Cassini spacecraft (Porco et al.
,**0), which makes Enceladus the fourth body in thesolar system (along with the Earth, Jupiter’s moon Io
and Neptune’s moon Triton) where volcanic eruptions
have been observed. This ﬁnding goes along withthe spectrometric detection of a distinctive warm spot
centered on the south pole with a temperature of
approximately 2/ K, which is +/ K more than expectedfrom a simple radiation balance (Fig. +; Spencer et al.
,**0).
Cryovolcanoes are icy equivalents of the wellknown terrestrial volcanoes. Their main features areillustrated in Fig. ,. Instead of magma, Enceladus’
cryovolcanoes erupt water, which has its source inpressurized sub-surface water chambers (Porco et al.
,**0). No ammonia, which would lower the melting
point, was found (Buratti et al. +33*), suggesting thatthe cryomagma is pure H, O. Cryovolcanoes on othericy moons may also erupt a mixture of water and, forexample, ammonia or methane.
The exact source which produces the cryovolcanic water plumes is still subject to research. In theliterature, di#erent mechanisms have been proposed
to explain the plume origin. Collins and Goodman* Now at Department of Applied Mechanics, Chalmers University of Technology, .+,30 Gothenburg, Sweden** Now at Chair of Applied Mechanics, University of Erlangen-Nuremberg, Egerlandstraße /, 3+*/2 Erlangen, Germany

24

Bulletin of Glaciological Research(,**1) suggest that localized subsurface melting on
Enceladus has produced an internal south polar sea,
whereas Nimmo et al. (,**1) ascribe the heat and vapour production to shear heating by tidally driven
lateral fault motion. Another possible explanation is
given by Kie#er et al. (,**0), stating that due to tectonic processes, a clathrate reservoir might be theplume’s source. In any case, the existence of thecryovolcanic eruptions is undoubted, and the details
on their source are not relevant for this study.
Most problems of heat transport in ice, as well as
in other common materials, are accurately describedby the classical theory developed by Fourier. However, at cryogenic temperatures material behavior candi#er greatly from that at room temperature. One ofthe properties that might change is the way heatpropagates. In the course of time, considerations led

Fig. +. Predicted (radiation balance) vs. observed (composite infrared spectrometer CIRS onboard theCassini spacecraft) surface temperatures of Enceladus. The south polar region shows a warm spot
due to active cryovolcanism. Credit: NASA/JPL/
GSFC.

to hyperbolic heat transfer equations that are able tomodel thermal wave propagation and to overcome theparadox of inﬁnite wave speed inherent in Fourier’s
theory. Among the di#erent approaches which havebeen suggested, the one of Green and Naghdi (+33+)
stands out by being very consistent and, most important for our case, very general. They formulated atheory which comprises Fourier’s classical law, a nonclassical law leading to a hyperbolic heat transfer
equation as well as a third type of heat conduction
which combines both limits. This theory has provento be well suited for modelling thermoelastic problems at cryogenic temperatures; see Bargmann andSteinmann (,**/, ,**0). Consequently, in this study
we apply the Green-Naghdi thermoelastic theory to
ice in order to simulate the phenomenon of cryovolcanism on Enceladus. We state explicitly that thishas so far not been directly supported by laboratory
studies or ﬁeld measurements on ice, so that the study
is clearly of a speculative nature and aims at exploring a promising possibility.
In +33+, Green and Naghdi published their thermoelastic approach which is capable of modelling awider range of phenomena than Fourier’s theory.
The latter is fully included in the Green-Naghdi approach ῌ however, more possibilities are admitted.
Thus, classical problems can be simulated as well as
thermoelastic phenomena at low temperatures; see
the works of Bargmann and Steinmann (,**/, ,**0).
The theory of Green-Naghdi is introduced in Sect. , inorder to familiarize the reader with this very generalapproach. Subsequently, a space-time ﬁnite element
discretization of the resulting initial boundary value
problem is sketched in Sect. -. This is followed by
setting up a simpliﬁed model for the material properties, geometry, initial and boundary values of Enceladus’ cryovolcanism (Sect. .). Finally, in Sect. /, computational results for the temperature evolution in theice surrounding the cryovolcano are presented anddiscussed for a number of di#erent scenarios.

,.

Fig. ,. “Cold geyser model” for cryovolcanism on
Enceladus: Pressurized sub-surface water erupts
through a volcanic pipe (vent) into the atmosphere.
Credit: NASA/JPL/SpaceScience Institute/CassiniHuygens.

Green-Naghdi thermoelasticity

The very general theory developed by Green andNaghdi was only published recently, but a number ofstudies devoted to its investigation have already beenreleased. Since in this contribution we restrict ourselves to the so-called Green-Naghdi theory of type III(see below), we list only some papers concerning thistype in the following. For example, Puri and Jordan(,**.) examine the propagation of plane waves in a
type-III-thermoelastic continuum. Surface waves inthermoelastic solids are studied by Kumar and Deswal(,**+). They apply the theory to an academic example of a micropolar half-space in order to build a basis
for the extension to study the role of Green-Naghditype surface waves in earth quakes and other geo-Bargmann et al.

physical applications. Bargmann and Steinmann(,**/, ,**0, ,**2) and Bargmann et al. (,**2) investigate
the general numerical behavior as well as the behavior of temperature proﬁles of the Green-Naghdi theory.
First, we reiterate the basic ideas of the GreenNaghdi theory in order to introduce its basic properties to the reader. The theory developed by Greenand Naghdi (+33+) is subdivided into three di#erent
types, labeled type I, II and III. The linearized theoryof type I is equivalent to Fourier’s law, resulting in aparabolic heat equation. In case of type II, no energy
dissipation is involved and the resulting heat transfer
equation is hyperbolic. Type III is the most general
type. It contains both type I and II as limiting cases,
thus it is capable of modelling classical Fourier heat
conduction as well as undamped thermal wave propagation ῍ and, in addition to that, many more phenomena. The heat ﬂux of type III is an extension of thoseof type I and II. Both types II and III can overcome
the unnatural property of Fourier’s law of inﬁnite
propagation speed and imply a ﬁnite wave propagation speed instead. For modelling cryovolcanism on
Enceladus, we apply the universal type III and couple
the heat transfer equation with the balance of linear
momentum. Type III seems to be the best choice as itis the most general case and includes the other two
types.
Thermoelasticity in a ﬁnite, isotropic and homogeneous body (an elastic deformable solid) B is considered. First, we summarize the basic concept brieﬂy,since we have already discussed the derivation of thethermal problem in detail in Bargmann and Steinmann(,**/), see also Green and Naghdi (+33+). Subsequently, the balance of linear momentum is introduced in order to describe the thermomechanical stateof body properly.
Green and Naghdi’s theory of type III is a very
general approach based on the introduction of a scalarquantity a deﬁned by῍a῎Xῌt῏῔῎ T῎Xῌt῏dtῒa*ῌ
῏῍ῌ

(+)

where X denotes the position vector in the reference
conﬁguration, t the time, T the absolute temperature,
and a* is the initial value of a at the reference time t*.
a is called thermal displacement, and
ῌ

a῔T

(,)

holds, where the dot represents the material time
derivative.
At this point, we only state the two equations
governing Green-Naghdi thermoelasticity of type III.
For a derivation of these equations, the reader isreferred to the Appendix. The thermal problem isgoverned by the heat transfer equation25

ῌ

rcT῔div῎k+῕a῏ῒdiv῎k,῕T῏ῒrrΐ-T* wKI: eῌῌ (-)
where r denotes the constant mass density, c is the
non-negative speciﬁc heat capacity, k+, , are the nonclassical and classical thermal conductivities, respectively, T* is the reference temperature, w the thermal
expansion coe$cient, K the bulk modulus, I the unity
tensor, r the external heat source and e the mechanicalstrain tensor. Furthermore, the colon operator (:) denotes double contraction. The last term in Eq. (-)
represents a coupling between the thermal and themechanical problem, i.e., mechanical strains lead totemperature changes.
The mechanical (deformation) problem is governed by the balance of linear momentum. For geometrically linear deformations, it reads
rü῔div῎E : eΐ-ῐTΐT*ῑwKI῏ῒbῌ(.)

where u is the mechanical displacement vector, E theelasticity tensor and b the volume force. Note thatEq. (.), like Eq. (-), includes a coupling between thethermal and the mechanical problem, so that a fully
thermo-mechanically coupled problem is at hand.

-.

Finite element discretization

In this section, we shortly describe the chosen
discretization method for the initial boundary value
problem arising from the Green-Naghdi thermoelasticity described above. The solution method is based
on ﬁnite elements for the spatial as well as for the
temporal problem. Bargmann and Steinmann (,**/,,**0) have already shown that this discretization approach suits well for Green-Naghdi thermoelasticity
of type III. Those readers who are interested in amathematical description, we refer to Bargmann andSteinmann (,**/, ,**0) where we explain the procedure in detail and list all mathematical equations. In
the following, only the basic steps are reiterated inorder not to repeat everything.
We apply a semi-discretization technique. Thismeans that the spatial and the temporal domain are
discretized one after the other. We start with thespatial one and use a standard Bubnov-Galerkin ﬁnite
element method. The so-called weak form of the equation of interest is one prerequisite. Therefore, theheat transfer equation (-) and the balance of momentum (.) are weighted by test functions and, subsequently, all equations are integrated over the reference conﬁguration B. The unknowns, i.e., the thermal
displacement a, the temperature T and the mechanicaldisplacement u, are interpolated element-wise by theshape functions Ni. Moreover, the test functions arealso discretized with the shape functions Ni, according
to the Bubnov-Galerkin method. Inserting the approximations results in semi-discretized equations
meaning that the nodal values of the unknowns are26

Bulletin of Glaciological Research

now only time-dependent.
After the spatial discretization, we apply a ﬁnite
element method for the temporal one as well. We
resort to a continuous Galerkin time ﬁnite element
method because it has proved to suit well for GreenNaghdi type III thermoelasticity in the past, see
Bargmann and Steinmann (,**/, ,**0). The reader isalso referred to the book of Eriksson et al. (+330). The
continuous Galerkin method approximates the trial
function piecewise and continuously in time by polynomials of degree k. The test functions are approximated piecewise by polynomials of degree k+ which
are discontinuous across the element boundaries.
The considered time interval I[t*, tf] (t*: initial time,
tf: ﬁnal time) is divided into a ﬁnite number of subin
tervals. In analogy to the spatial ﬁnite element method, the weak form is a prerequisite for the discretization process. Thus, we multiply by the temporal test
functions and integrate over the corresponding time
subinterval. Afterwards, the ﬁnite element approximations of the unknowns, i.e., the thermal displacement a, the temperature T and the mechanical displacement u, and those of the test functions are inserted into the system of equations. As a result, we
obtain a fully discrete algebraic set of equations.
We seek for the temporal weak form of the semidiscrete equations on one time subinterval after theother. The solution of the preceding subinterval isused as an initial condition for the current one. In
each time step, the coupled system is solved monolithically.

..

Set-up for modelling cryovolcanism
on Enceladus

First, the material parameters for cryogenic ice on
Enceladus are introduced  followed by the description of the computational domain, the initial and theboundary conditions.
Model parameters for cryogenic ice
According to Powell (+3/2), Slack (+32*) and Greve(,**0), the thermal conductivity, the thermal expansion coe$cient and the speciﬁc heat capacity of ice are

the ice density, we employ the standard value r3+*
kg m-. While Buratti et al. (+33*) state that Enceladus’ ice is indeed pure H,O, newer investigations tend
to assume the existence of very small portions ofother substances, see Matson et al. (,**1). However,even in the latter case, the above value for r seems
reasonable. Furthermore, the acceleration due togravity on Enceladus is set to g*.*13 m s,.
The Green-Naghdi theory of type III includes thetwo thermal conductivities k+ and k,, to which weshall refer as the non-classical and classical thermal
conductivities, respectively. Green and Naghdi (+33+)
do not narrow the ratio of k+ and k,. Since theydescribe thermal conductivities (or at least somethingsimilar), they should be non-negative, i.e., k+, k,ῌ*.
Having another look at the heat transfer equation oftype III, Eq. (-), we see that if k+*, the remaining
equation is equivalent to the classical heat transfer
equation resulting from Fourier’s law. Thus, the classical theory is fully incorporated in the Green-Naghdi
theory of type III. The other extreme case, i.e., k,*,leads to an undamped wave equation, and that is what
Green and Naghdi refer to as type II. One of its
properties is that the constitutive equations are usually chosen in such a way that they do not include anyenergy dissipation. Such a property might be reasonable for modelling the pure second sound phenomenon, but even in this special case, the debate
whether or not energy is dissipated is still going on(Green and Naghdi +33-, Quintanilla and Straughan,***).
Although no experimental data exist on energy
dissipation during the eruption of a cryovolcano on
Enceladus, one can certainly conclude that the model
should include energy dissipation due to the experiences on Earth, and due to the fact that so far a purethermal wave propagation has only been proven below ,* K. We therefore identify the classical thermal
conductivity k, with the measured thermal conductivity given in Eq.(/)+,

..+

0+/
Wm+K+
T+*0 K+
wT /0./ *.,/ T,1+ +cT +.0.- 1.,/-T Jkg K 

k,T

where the temperature T is taken in units of [K].
Natural ice is composed of a large number of individual crystals, and we assume that the orientation distribution is random. Thus, its macroscopic mechanicalbehavior of the ice is isotropic, and the elasticity
tensor E is fully determined by Young’s modulus, E
+* GPa, and Poisson’s ratio, n*.--. Both quantities
are assumed to be independent of temperature. For(0)

For the non-classical thermal conductivity k+, we assume that it is proportional to k,,

kT

(/)

0+/
Wm+K+T

k+T

k,T
tproc

f

(1)

where tproc is a time-scale for the modelled process, andf is an adjustable, dimensionless factor. It is clear
that, for small values of f, the Fourier-type (di#usive)
heat ﬂux will be dominant, while for large values of f,heat transport will be governed mainly by the nonclassical a-type ﬂux.
In order to constrain possible values for the nonclassical thermal conductivity k+, we employ the following line of argumentation. One of the beneﬁts ofBargmann et al.

27

the thermoelastic theory of Green and Naghdi is that
it is capable of modelling the physical low temperature phenomenon of second sound. For materials
which show this phenomenon, there exists a slower
thermal wave called “second sound” in addition to thewell-known mechanical (elastic) wave termed “ﬁrstsound”. The velocity of the thermal wave was foundto be typically about ῒ῍ times smaller than the velocity of the elastic wave (Bargmann and Steinmann,**/, ,**0). If we apply this approximation for ice, we
obtain an independent estimate for the value of k+,

Fig. -. Set-up of the model cryovolcano. The adjacent ice layer is /* km long and /** m thick. The
dashed square marks the /** m/** m area next tothe volcanic vent for which the evolution of thetemperature ﬁeld will be plotted in Figs. .ῌ1.

uthermalῌuelastic ῌῒ῍
ῐῌ
῎
(2)
ῌ
ῌῒ῍
ΐ ῑ῏ ΐ ῑEck+ῌ
For +** K, this yields k+,.3+*+, Wm+ K+, whilek,0.+/ Wm+ K+. If we choose the process time-scale
as tproc+ d20.** s (see below), we obtain a factor f aslarge as f.+*+0. However, since it is not likelythat ice shows fully developed second-sound phenomena (otherwise it would most likely have been observed already), we argue that reasonable values for fmust be at least several orders of magnitude smaller.
..,Computational domain, initial and boundary
conditions
As a simpliﬁed representation of a cryovolcano
and its environment on Enceladus, we consider a twodimensional domain of /* km (approx. +* degrees oflatitude) length and /** m thickness, with the volcanic
vent at the right edge. The problem is supposed tobe symmetric around the vent, and the curvature of
the surface is neglected. We assume that the waterchamber next to the volcanic pipe is + km long andover-pressurized by a factor of +* compared to thehydrostatic pressure of the ice. The set-up is illustrated in Fig. -.
The initial conditions at the time t** are a surface
temperature of 2/ K (see Fig. + and its discussion in the
introduction) and a bottom temperature of ,1- K. The
reference temperature T* required for the solution of
the heat transfer equation (-) is set to the initial surface temperature, that is, T*2/ K. The simulated
eruption starts abruptly at this time and lasts untilt+ d. During the eruption period, the temperature
over the entire depth of the volcanic vent is set to ,1- K,
and the pressure in the vent falls linearly from theover-pressurized water chamber at the bottom to thezero value at the surface. At t+ d, the eruption is“switched o#”, and the computation is continued untiltf+* d. For the process time-scale tproc in Eq. (1), wechoose the duration of the eruption, that is, tproc+ d.
There are no data available on the duration of areal eruption, so that the choice of one day is arbitrary. However, the general development of the heatpropagation will not be a#ected by the eruption time.
Here, we want to discuss the e#ect of the GreenNaghdi model, and therefore, we do not focus on
details of the temporal evolution.
The surface temperature Ts is computed from aradiation balance, which takes into account the incoming solar radiation, the black-body radiation fromthe ice surface and the geothermal heat ﬂux approaching the ice surface from below. It reads
S*+A.sSBT .s .qy*(3)

where S* is the solar constant for Enceladus, A thealbedo (reﬂexivity of the ice surface), sSB the StefanBoltzmann constant (/.01+*2 Wm,K.) and qy the
upward geothermal heat ﬂux. As already mentioned
in the introduction, the albedo is as high as *.33, whichis the highest value of any known body in our solar
system (Buratti et al. +33*). The solar constant can becomputed from its terrestrial counterpart, S*, Earth+-02 Wm,, via
S*S*, Earthr,SE+..3Wm,
r,SS

(+*)

where rSE and rSS are the distances Sun-Earth (+ AU) and
Sun-Saturn (3./2 AU), respectively. Note also that theassumed thickness of the ice layer of /** m can bejustiﬁed by inserting a surface temperature of 2/ K anda bottom temperature of ,1- K into the radiation balance (3) and expressing the heat ﬂux with the classical(Fourier-type) part of Eq. (+.), as qyk, T/ z.

/.

Results and discussion

The dimensionless factor f in Eq. (1) still needs tobe speciﬁed. Since the available data on Enceladus
do not give a hint on how to choose f in the ﬁrst place,
we have decided to vary this parameter starting withf+ and going up to f+*2. As mentioned above, forsmall values of f, our model will be close to the classical, Fourier-type heat conduction. By increasing f,
the model will become more and more non-classicalbecause of the increasing impact of the a-type ﬂux.
However, note that the maximum value f+*2 is still
eight orders of magnitude smaller than the estimate

28

Bulletin of Glaciological Research

given in the discussion of Eq. (2), so that in neithercase we assume fully developed second-sound properties for ice.
In the following, the solutions for f῏+ (“case +”),
f῏+*0 (“case ,”), f῏+*1 (“case -”) and f῏+*2 (“case .”) will
be discussed. All simulations have been carried out
with a numerical time-step of ῑt῏+*῍- d῏+... min.
Figures .ῌ1 illustrate the temperature development of the /** m῎/** m square next to the volcanic
pipe. The solutions are plotted at four di#erent
times of the observation for each of the four cases.
First, the temperature is depicted right after the beginning of the volcanic eruption (i.e., at time t῏+...
min). The next plot shows the temperature distribution + hour and +, minutes after the beginning of the
volcanic eruption ῌ at that time the volcano is still
erupting. Immediately after the end of the eruption,i.e., after + day, we look at the temperature solution
for the third time. Fourth, the temperature is plotted+* days after the beginning of the eruption, that is, 3days after its end.
In case +, which is very close to the classicalFourier approach, no change in temperature can beseen, see Fig. .. Neither the +-day-lasting volcanic
eruption nor the cooling down period of 3 days haveany noticeable impact on the temperature distribution within the ice.
In case ,, where f῏+*0, hardly any changes intemperature can be seen during the volcanic eruption,
see Fig. /. However, 3 days after the end of the eruption, it can be clearly seen that the temperature in the
volcanic pipe is decreasing, starting at the top, i.e., at
Enceladus’ surface.
Figures 0 and 1 show that, with an increasing
impact of the non-classical ῐa-type ﬂux (cases - and.), the volcanic eruption starts to heat the interior of
the ice next to the pipe, and the surface temperaturerises. Also, immediately after the end of the eruption, the ice starts to cool down, which is best visible
in the volcanic pipe next to the surface. Dependingon the dimensionless factor f, the ice temperatureequilibrates slower or faster: The higher the value off, the faster the system equilibrates. Furthermore,
due to the increasing impact of the ῐa-type heat ﬂux,

Fig. .. Case + (f῏+): Evolution of the temperature ﬁeld (in K) in the /** m῎/** m area next to thevolcanic vent shown in Fig. -.

Bargmann et al.

29

Fig. /. Case , (f῍+*0): Evolution of the temperature ﬁeld (in K) in the /** mῌ/** m area next to thevolcanic vent shown in Fig. -.

the temperature rises within the ice during the eruption. Owing to the radiation balance (3), the temperature change is most pronounced in the interior of theice layer and smaller at the surface; however, a surface warming of the order of several degrees occurs in
the vicinity of the vent during the eruption and fades
away afterwards.
Spencer et al. (,**0) discuss the temperature on
Enceladus’ surface observed during the Cassini ﬂybys in ,**/ based on infrared spectrometry. Highresolution data reveal that the occurrence of hightemperatures is spatially correlated with the geological features termed “tiger stripes” (linear troughs),
which are most likely the source of the cryomagmadelivered during the volcanic eruptions. With regard to the results of our simulations, it is particularlyinteresting that the data also indicate some increased
thermal emission from regions adjacent to the warm
troughs (several kilometers away). This spreading of
surface warming agrees roughly with the results ofour cases - and ., where the non-classical, ῎a-type

heat ﬂux contributes signiﬁcantly to the heat transfer
in the ice. By contrast, in case +, which is very close
to the classical Fourier approach, the warming extends to only a few meters away from the volcanic
pipe during the simulation time (and is therefore invisible in Fig. .). Of course, this depends on the assumed duration of the eruption, but the time-scale [t]for purely di#usive heat transport over a length-scaleof [L]῍+ km is as large as [t]῍rc [L],/k,ῌ.*** a (computed for +** K), and it appears unlikely that a stationary cryovolcanic eruption can be sustained over sucha long time. Therefore, with all due caution arisingfrom our simpliﬁed model and the data uncertainties,
it seems that non-classical heat transport in ice at
cryogenic temperatures may play a role in explaining
the observed temperature distribution in the vicinityof the volcanically active troughs in Enceladus’ southpolar region.

30

Bulletin of Glaciological Research

Fig. 0. Case - (f῍+*1): Evolution of the temperature ﬁeld (in K) in the /** mῌ/** m area next to thevolcanic vent shown in Fig. -.

0.

Conclusions

Motivated by the fully consistent and very general thermoelastic theory of Green and Naghdi (+33+),
the aim of this contribution was the modelling ofcryovolcanism on Enceladus according to their approach. Since their theory of type III includes theclassical thermoelastic theory as a limiting case (k+῍*)and allows for various non-classical models, we stickto this Green-Naghdi sub-theory. Type III allows thermal waves to propagate without energy dissipation
at ﬁnite speed, another limiting case. We modeled
cryovolcanism on Enceladus applying di#erent scenarios lying in between the two extreme cases aforementioned.
First, we reiterated the basic equations of thermoelasticity based on the approach of Green and Naghdi.
Second, we introduced a ﬁnite element method for thediscretization of the governing system of equations.
This is followed by a detailed description of the ap-

plied ice volcano model. Finally, we present numerical results for the Green-Naghdi theory type III for
four di#erent cases with increasing inﬂuence of the
non-classical heat ﬂux. While the approach closest to
the classical theory (our case +) does not allow for any
signiﬁcant spreading of surface warming away fromthe active volcanic vents, the scenarios where thenon-classical heat ﬂux plays a role (cases - and .) show
surface warming of several degrees in the vicinity ofthe vents. This seems to be in better agreement withtemperature data based on infrared spectrometry(Spencer et al. ,**0). Therefore, we tentatively andcarefully conclude that a heat transport modelled by amore general approach than Fourier’s may be signiﬁcant in ice at the very low temperatures which
occur on Enceladus and other icy moons in the outersolar system. Of course, complementary laboratory
studies would be required in order to further support
this idea.

Bargmann et al.

31

Fig. 1. Case . (f+*2): Evolution of the temperature ﬁeld (in K) in the /** m/** m area next to thevolcanic vent shown in Fig. -.

Appendix
Derivation of the governing equations for thermoelasticity of Green-Naghdi type III
The heat transfer equation (Eq. (-)) follows fromthe balance of entropyῌrhdivhrsx(++)

where h, h, s and x denote the speciﬁc entropy (entropy
per mass unit), the entropy ﬂux vector, the external
entropy source and the internal rate of entropy production, respectively. Bargmann and Steinmann(,**1) have shown that the classical thermodynamic
relation between the entropy ﬂux vector h and theheat ﬂux vector q,
Thq

(+,)

can be transferred to the non-classical case.
As mentioned above, one of the main features ofGreen and Naghdi’s theory of type III is the introduction of a new thermal variable, namely the thermal

displacement a. In case of geometrically linear thermoelasticity of type III, the speciﬁc free energy (freeenergy per mass unit) y and the heat ﬂux vector q
depend on the thermal displacement a, the temperature T, the thermal displacement gradient a, thetemperature gradient T and the mechanical strain
tensor e,yy a T a T eqq a T a T e(+-)

The relations for y and q, based on those concluded by
Green and Naghdi (+33+), read+ry e : E : e-TT*wKI : e,
rc TT*,
TT* rh*
,
T*

(+.)

qk+ak,T
where h* is the reference value of the speciﬁc entropy
for TT* and e*. The nonnegative speciﬁc heat
capacity c is equal to cT h / T. Moreover, theentropy reads

32

Bulletin of Glaciological Researchrhry rc TT* -wKI : e rh* (+/)
TT*

and the internal rate of entropy production x equals
rTxhῌT;(+0)

see Bargmann and Steinmann (,**1). Thus, energy is
dissipated during general processes.
The heat transfer equation of type III is derived
by inserting the constitutive assumptions (+.) and therelations (+/) and (+0) into Eq. (++) and, following Greenand Naghdi (+33-), subsequent replacement of T by T*.
After some calculations, we ﬁnally obtain
῍

rcTdivk+a divk,T

rr-T* wKI : e῍ (+1)

where the heat source r is related to the external
entropy source s via rTs.
The balance of linear momentum states that therate of change of the momentum rü is equal to thesum of the divergence of the Cauchy stresses s andthe volume force b,rüdivs b(+2)

The stresses s are thermodynamically conjugated tothe strains e via the speciﬁc free energy y,
s

ry
e

(+3)

Thus, the constitutive equation for the stresses readssE : e-TT*wKI(,*)

including thermo-mechanical coupling due to the lastterm. Inserting Eq. (,*) into the balance of momentum (+2) ﬁnally yields Eq. (.),
rüdivE : e-TT*wKI b(,+)

Acknowledgments
The authors wish to thank N. Azuma and an anonymous referee who provided many helpful comments
and suggestions. This work was supported by theGerman Research Foundation (Deutsche Forschungsgemeinschaft, DFG) under grant STE /../,-. Partsof the work (mainly the numerical simulations) werecarried out during a visit of S.B. at the Institute of
Low Temperature Science in Sapporo, supported by aresearch fellowship from the Japan Society for the
Promotion of Science (JSPS).

References
Bargmann, S. and P. Steinmann. ,**/. Finite element approaches to non-classical heat conduction in solids.
Comp. Meth. Eng. Sci., 3 (,), +--ῌ+/*.

Bargmann, S. and P. Steinmann. ,**0. Theoretical and computational aspects of non-classical thermoelasticity.
Comp. Meth. Appl. Mech. Eng., +30 (+ῌ-), /+0ῌ/,1.
Bargmann, S. and P. Steinmann. ,**1. Classical results for a
non-classical theory: remarks on thermodynamic relations in Green-Naghdi thermo-hyperelasticity. Cont.
Mech. Thermodyn., +3 (+ῌ,), /3ῌ00.
Bargmann, S. and P. Steinmann. ,**2. An incremental variational formulation of dissipative and non-dissipative
coupled thermoelasticity for solids. Heat and Mass Transfer., ./, +*1ῌ++0. doi: +*.+**1/s**,-+ῌ**2ῌ*.*/ῌ/.
Bargmann, S., P. Steinmann and P.M. Jordan. ,**2. On the
propagation of second sound in linear and nonlinear
media: Results from Green-Naghdi theory. Phys. Lett. A,-1, (,.), ..+2ῌ..,..
Buratti, B. J., J.A. Mosher and T. V. Johnson. +33*. Albedo andcolor maps of the Saturnian satellites. Icarus, 21 (,), --3ῌ-/1.
Collins, G.C. and J.C. Goodman. ,**1. Enceladus’ south polarsea. Icarus, +32 (+), 1,ῌ2,.
Eriksson, K., D. Estep, P. Hansbo and C. Johnson. +330. Computational di#erential equations. Studentlitteratur, Lund,
Sweden.
Green, A.E. and P.M. Naghdi. +33+. A reexamination of thebasic postulates of thermomechanics. Proc. R. Soc. Lond.
A, .-, (+22/), +1+ῌ+3..
Green, A.E. and P.M. Naghdi. +33-. Thermoelasticity without
energy dissipation. J. Elasticity, -+ (-), +23ῌ,*2.
Greve, R. ,**0. Fluid dynamics of planetary ices. GAMMMitteilungen, ,3 (+), ,3ῌ/+.
Kie#er, S.W., X. Lu, C.M. Bethke, J.R. Spencer, S. Marshak
and A. Navrotsky. ,**0. A clathrate reservoir hypothesis
for Enceladus’ south polar plume. Science, -+. (/2*0),+10.ῌ+100.
Kumar, R. and S. Deswal. ,**+. Disturbance due to mechanical and thermal sources in a generalized thermomicrostretch elastic half space. Sadhana, ,0 (0), /,3ῌ/.1.
Matson, D.L., J.C. Castillo, J. Lunine and T.V. Johnson. ,**1.
Enceladus’ plume: Compositional evidence for a hot interior. Icarus, +21 (,), /03ῌ/1-.
Nimmo, F., J.R. Spencer, R.T. Pappalardo and M.E. Mullen.
,**1. Shear heating as the origin of the plumes and heatﬂux on Enceladus. Nature, ..1 (1+.,), ,23ῌ,3+.
Porco, C.C., P. Helfenstein, P.C. Thomas, A.P. Ingersoll, J.
Wisdom, R. West, G. Neukum, T. Denk, R. Wagner, T.
Roatsch, S. Kie#er, E. Turtle, A. McEwen, T.V. Johnson,
J. Rathbun, J. Veverka, D. Wilson, J. Perry, J. Spitale, A.
Brahic, J.A. Burns, A.D. DelGenio, L. Dones, C.D. Murray
and S. Squyres. ,**0. Cassini observes the active southpole of Enceladus. Science, -++ (/100), +-3-ῌ+.*+.
Powell, R.W. +3/2. Preliminary measurements of the thermal
conductivity and expansion of ice. Proc. R. Soc. Lond. A,,.1 (+,/+), .0.ῌ.00.
Puri, P. and P.M. Jordan. ,**.. On the propagation of planewaves in type-III thermoelastic media. Proc. R. Soc. Lond.
A, .0* (,*/+), -,*-ῌ-,,+.
Quintanilla, R. and B. Straughan. ,***. Growth and uniqueness in thermoelasticity. Proc.R. Soc. Lond. A, ./0 (+332),+.+3ῌ+.,3.
Slack, G.A. +32*. Thermal conductivity of ice. Phys. Rev. B, ,,(0), -*0/ῌ-*1+.
Spencer, J.R., J.C. Pearl, M. Segura, F.M. Flasar, A. Mamoutkine,
P. Romani, B. J. Buratti, A.R. Hendrix, L. J. Spilker and R.
M.C. Lopes. ,**0. Cassini encounters Enceladus: Background and the discovery of a south polar hot spot.
Science, -++ (/100), +.*+ῌ+.*/.