Analysis of 3D numerical simulations of subsolidus thermal convection: application to Venus and Europa

By Ian W. Bolliger ‘11

Geophysics and Planetary Geosciences, Jet Propulsion Lab
National Aeronautics and Space Administration and California Institute of Technology

Numerical simulations using 3D spherical and Cartesian coordinates are carried out to analyze convective properties of fluids with viscosities of varying dependence on temperature and varying Rayleigh numbers. Several MATLAB codes are developed in order to interpret the results of the numerical models. Scaling laws are predicted and compared to those of studies involving 2D and/or isoviscous simulations. The observed ratio of ΔTTBL= ~1.43ΔTμ, where ΔTTBL is the temperature difference across the hot thermal boundary layer and ΔTμ is the viscous temperature scale, correlates well with 2D predictions. Convective stress for a range of planetary parameters is then compared with shear strength within a planet’s conductive lid to determine whether the vigor

of convection can account for lid faulting. In all but one parameter set, stress data indicates that lid fracturing should occur on Venus, which conflicts with observations that Venus is in a stagnant lid regime. In the one exception, a simulation with a coefficient of friction similar to laboratory predictions for Earth’s mantle produces stresses small enough for the lid to remain intact. High Rayleigh numbers in Europa’s icy shell have led to slow convergence in the simulation, so more work is needed to obtain significant results for Europa.


Several evolving methods of modeling planetary convection processes exist today. Until recently, none of these models has incorporated a temperature dependent viscosity in a three-dimensional spherical space. In the current study, a finite-volume multigrid numerical simulation called ŒDIPUS1 incorporates these qualities. ŒDIPUS is used to model the interiors of Venus and Europa, and the results are analyzed to produce qualitative and quantitative predictions regarding scaling laws as well as to draw conclusions about the behavior of the lithosphere of Venus and the icy shell of Europa.

The outer shell of Europa is an H2O layer with estimates of its thickness ranging up to 120 km2. While it is not known exactly how much of this layer is ice and how much is water, surface fractures are observed in images from Galileo, which suggests that perhaps the ice layer can be broken. If this is the case, and if, as has been postulated, the Europan ocean has the necessities for primitive life3, it is possible that any organic material that exists in the ocean could be transported to the surface through these fractures. The Europa Jupiter System Mission (EJSM), scheduled for launch in 2020, could investigate this possibility; however, a mechanism for the observed fractures, and thus for the possible organic transport, must first be identified.  One possible mechanism is convective stress, and the ŒDIPUS model seeks to provide insight on this possibility. Unfortunately, at the time of this report, fewer results have been obtained than are necessary to draw any significant conclusions on the strength of convective stress within Europa. However, methods for doing so are now in place for when the data does arrive. These will be described later in the report.

Venus, though clearly quite different from Europa, also provides a strong motivation for stress analysis. Though Earth-like in size and composition it does not currently experience plate tectonics. It has been hypothesized that plate tectonics on a given planet exist when convective stresses are greater than the yield strength of the material in a planet’s no-slip lithosphere. This study analyzes the stress and strength values produced for a number of simulations that utilize a range of possible Venusian parameters. By doing so, we more accurately confine the parameter space that should produce the observed “stagnant lid” and provide insight into Venus’ interior dynamics.

In all but one simulation, the study produces Venusian convective stresses that imply a “mobile lid”. However, when the coefficient of friction of the mantle is raised to a value equivalent to that determined from laboratory experiments of Earth-like materials and the rest of the planetary parameters are configured to provide low stress values, a situation is obtained in which the lithosphere does not fracture.

Scaling laws, while not directly dependent on the use of Europan or Venusian parameters, also provide useful information about the behavior of planetary interiors. This study compares laws demonstrated in its own results to those previously determined for models utilizing other geometries and viscosity approximations. It finds that the ratio of temperature difference across the bottom thermal boundary layer to viscous temperature scale (see Section II) observed in a 2D domain applies to 3D space as well, and that the relationship between Nusselt number and Rayleigh number obtained in previous studies does not apply when varying the temperature-dependence of viscosity. Finally, where previous scaling laws do not exist, the study provides qualitative explanations for the observed behavior of convection-related values when varying planetary parameters. Before introducing the methods used in the study, a brief introduction to convection processes and an explanation of terms used throughout the report is now provided.

There are four layers of interest in a convecting space like a planetary interior (Figure 1). The lowermost layer is the hot thermal boundary layer where hot plumes form and convect upwards. The instabilities of this layer, which create the hot plumes, are what drive the convection of the entire layer. On top of this layer is the convective domain, where hot plumes are rising and cold plumes are falling (through both advection and diffusion). Tavg in this layer is adiabatic, though in ŒDIPUS the effects of hydrostatic pressure are not incorporated and thus Tavg is constant in the model. Following this section is the cold thermal boundary layer, where cold plumes form and descend. Finally, in models that incorporate a large enough viscosity gradient, the top layer is the conductive lid. In this area, viscosities are high enough that advection is negligible. The temperature gradient in this layer is linear, with heat transported solely by conduction.


Figure 1. Vertical profile of Tavg for a fluid with a strongly temperature-dependent viscosity in a model that does not incorporate the effects of pressure.



ŒDIPUS, the 3D multigrid flow solver used in this study, is the model described in Ref. 1. For simulations of Venus, in which the depth of the mantle is relatively large compared to the radius of the planet, a spherical model is used in order to realistically represent the convecting region. For simulations of Europa, in which the icy shell represents a small portion of the planet’s radius, a simpler, faster, Cartesian model is used. Both models utilize a finite-volume grid with 128x128x128 elements to numerically integrate the Navier-Stokes equations. The Europan model has a 2x2x1 aspect ratio in the x,y,z directions, respectively, while the width in the xy-direction for the spherical model is related to the curvature of the layer, defined using the ratio rc/rp. One spherical simulation with Europan parameters was attempted in order to compare the results with the Cartesian runs; however, the spherical model encountered convergence issues associated with a large rc/rp value and did not run successfully. We used multiple runs to test free-slip and no-slip boundary conditions for both the spherical and Cartesian models. Ideally, one might want to set a no-slip bottom boundary (representing the mantle/core interface) and a free-slip upper boundary (representing the lithosphere/atmosphere interface). However, because the simulation is currently unable to incorporate heterogeneous boundary conditions, we make inference from the two homogeneous possibilities.

Figure 2. Comparison of the Arrhenius and Avis viscosity laws for a spherical, Venus-based model with avis=20, DT=840 K. The three horizontal lines represent, from top to bottom, the lower boundary of the conductive lid, the lower boundary of the cold thermal boundary layer, and the upper boundary of the hot thermal boundary layer.


To better understand ŒDIPUS, a brief and limited description of the model is provided here. The main input parameter varied between runs was the temperature dependence of the viscosity. In ice and silicate materials, viscosity has a dependence on temperature, shear stress, and grain size. With the assumptions of a Newtonian solid and a constant grain size, viscosity can be approximated with the Arrhenius law:


Unfortunately, this equation leads to large viscosity differences on the scale of 1015 Pa s, which cause divergences in numerical simulations. Therefore, ŒDIPUS utilizes the following approximation:


where h0 is the viscosity at T=T0, and avis is a constant that is related to H through the viscous temperature scale, defined in the following equations:


The non-Newtonian effects on viscosity have previously been investigated, but their effect is much less than that of temperature4,5. Equations (1) and (2) yield similar viscosities in the convective layer but diverge in the thermal boundary layers (Figure 2), a difference that is discussed further in the stresses section.

In order to overcome convergence issues in the model, it was necessary to begin runs with a value of avis lower than that calculated in Equation (3). Specifically, the model was first run with avis=1. Once steady state was reached for a given value of avis, the constant was increased by 1-4. This process was repeated with the goal of reaching the value calculated using Equation (3) for each planet. For Venus, this value is around 20, depending on assumptions made about the enthalpy involved in the deformation of silicates (H), and it was reached in our study. For Europa, the value is slightly less, around 15, but has not yet been reached due to slower steady state convergence caused by high Rayleigh numbers.

Steady state is defined for both models by the difference in Nu values at the top and bottom of the layer. Nu is a non-dimensional value that, in Cartesian coordinates, represents the ratio of convective heat flux (qconv, the sum of advective and diffusive fluxes) to what the conductive heat flux would be if there were no advection {qcond=k(T1-T0)/d}. It is defined as

, where  

        (4, 5),

and is the horizontally averaged temperature at the radius of a given grid element. In spherical coordinates it no longer represents this ratio but rather a non-dimensional value that can be compared to heat flux. In this case it is defined as


For the Cartesian case, we define steady state as a case in which the difference in Nu values was less than 1%. For the spherical case, due to the r2 difference in surface area from the core-mantle boundary (CMB) to the surface, where r=rc/rp, it is necessary to multiply Nutop by r2 before comparing it to Nubot. This difference in surface area causes complications in heat flux convergence; thus, we consider all spherical runs with a difference in Nu values of less than 5% to have reached steady state.

In addition to avis, each model involves one more parameter that is varied from run to run. This value, the “surface” Rayleigh number, is a ratio of the buoyancy forces to the viscous forces in the entire layer and thus corresponds to the vigor of convection. It is calculated using the following equation:


For each value of avis, h0 is not defined directly. Rather h1, the bottom viscosity, is defined for all avis values. For Venus,  h1= 1021 Pa s. This value is based loosely on an Earth analog6. The Europan runs are performed with h1 = 1014, based on findings in Ref. 7. Equation (2) is then inverted to determine h0. Ra varies when avis is changed due to the change in h0.

For each element, ŒDIPUS outputs 8 values. These are comprised of velocities in each of the three dimensions on both the top and bottom edges of the cuboid perpendicular to that dimension in addition to temperature and pressure in the center of the element (Figure 3). All values are non-dimensionalized using the factors in Table 1. To obtain velocities in the center of the elements (so that ux, uy, and uz are all taken from the same point), the mean of the velocities on either edge of the cuboid is taken. Thus, the adjusted output involves 5 values for each element, all calculated at the center of the cuboid. These values are the velocity in each dimension, the pressure, and the temperature. In our study, only the temperature and the velocities are used. Once steady state has been reached for a given avis value, the ŒDIPUS outputs are then run through a variety of MATLAB programs created and/or modified for this study in order to obtain the desired data. These programs are described in Section III.

The interior limits of the thermal boundary layers that develop in the model are defined as the depths where the rates of advection and diffusion are equal. The results of the study show that this definition corresponds well with another definition used frequently in previous studies: the point where a linear approximation of the vertical temperature gradient of the layer (qconv/k) extended toward the interior intersects with Tm8. This correlation is especially accurate for the bottom thermal boundary layer. For the top thermal boundary layer it is accurate for low values of avis and begins to become less so with increasing temperature-dependence of viscosity (Figure 4).

The boundary between the conductive lid and the upper thermal boundary layer is much tougher to define. After deliberating over several possible definitions, it was decided that this boundary should be at the point where the maximum vertical velocity at a given depth was equal to 1/1000th of the maximum vertical velocity for the entire layer. This was based on a qualitative analysis, which suggested that this definition placed the boundary as close as possible to the point where V0 in all dimensions.

Figure 3. Depiction of the ŒDIPUS outputs.


Table 1. Non-dimensionalization factors.
Figure 4. Vertical profiles of spherical, Venus-based simulation at steady state for avis=5 and avis=20, respectively. The black line represents the base of the conductive lid, calculated as the point where Vr,max for a given depth is equal to 1/1000th of Vr,max for the entire space. The green and red lines represent the interior limits of the top and bottom thermal boundary layer, respectively. The purple lines on the Tavg plot represent the interior extension of the vertical temperature gradient at the top and bottom of the layer, while the vertical blue line represents Tm. The correlation between the two boundary layer definitions at the top thermal boundary layer is much greater in the low avis case. All velocities and temperatures are non-dimensional.

MATLAB programs

Several programs are needed in order to interpret the ŒDIPUS data. Some had been created by previous researchers and required only minor modifications, while two were developed specifically for this study. These two deal with viscous stress and scaling law analyses.

Viscous stress analysis

Viscous stress in the ŒDIPUS grid is calculated using the deviatoric stress tensor,


It is true that viscosities in the conductive lid, where faulting is likely to occur, will be underestimated by the Avis law (Figure 2). However, this low viscosity will lead to similarly overestimated velocities, and the values for tvisc should not be significantly affected.

Both tx,z and ty,z are calculated for each gridvolume in the ŒDIPUS model except for the most exterior points. This exception occurs because the velocity derivative calculation involves taking the difference in velocities on either side of a given gridvolume, and this is not possible when considering the outermost gridvolumes.

Four derivatives are calculated for each of these grid elements: ux/¶z, ¶uz/x, ¶uy/¶z, and ¶uz/¶y. These derivatives are calculated as follows:



The ¶uy/¶z, and ¶uz/¶y values are calculated in a similar fashion. With these four derivatives and with the viscosity for each gridvolume, it is possible to calculate both deviatoric stress tensors at each point using Equation (8). We obtain the total strain rate () and viscous shear stress (tvisc) for each point by finding the sum of the two strain and the two stress vectors:


Next, to determine the stress on the conductive lid, the deviatoric horizontal normal stress is calculated. To do this, we follow the approach of Ref. 9, which suggests that this stress multiplied by the depth of the conductive layer must balance the shear stress times the distance over which it is acting (the distance between hot and cold plumes). This correlation gives Equation (12):


D is approximated using the assumption that plumes are symmetrically spaced in the planform. This is typically true for runs at higher avis values (Figure 5), which are the most useful runs to evaluate when analyzing lid faulting. At lower avis values, when the planform is less symmetric, the approximation provides a fairly simple method to obtain the average distance between plumes yet introduces a significant source of error in yield strength calculations.

The method for counting plumesis as follows: first the number of plumes at the base of the lid is counted. The boundary of a plume is defined as the points where the temperature difference from the Tavg is 40% of the maximum temperature difference for that gridlayer. Using this definition for plume boundaries, it is possible to systematically identify hot plumes at the base of the lid. Because cold plumes are in the early stages of formation at this depth, it is often necessary to look at a gridlayer closer to mid-layer depth in order to identify the number of cold plumes. Once the total number of plumes is determined using this method, the following equation is used to find D:


After D and s are determined, s is compared to the yield strength of the ice layer, in Europa’s case, or the lithosphere, in Venus’ case, to determine if the viscous stress is large enough to cause faulting. The yield strength as a function of depth is calculated using the law proposed in Ref. 10,


where z is depth. Values of C0 range from 25 to 100 MPa for silicates. For ice, according to laboratory experiments11, C0=1 MPa and m=0.55 from the surface until pressures reach ~10 MPa, after which the coefficients are equal to 8.5 MPa and 0.2, respectively. For this study, a C0 value for silicates of 80 MPa and values of m from 0.05 and 0.2 are used. The C0 value is in accordance with Ref. 12, while the two m values are analogous to estimations of Earth’s interior. While laboratory experiments predict values near 0.2, the behavior of Earth’s mantle is reproduced in numerical simulations with values around 0.05.

The MATLAB program developed to analyze stress and yield strength plots vertical profiles of seven values related to these two properties. These values are: maximum viscous stress, RMS stress, Maximum strain, deviatoric horizontal normal stress, yield strength, viscosity (using both the Avis and the Arrhenius approximations), and, once again, average temperature. In this case, all values plotted are dimensional. The program plots these values for the entire layer and also magnifies the thermal boundary layers, plotting the values for these areas in a separate figure. Finally, the program creates two more figures in which it overlays a vector field indicating the direction and magnitude of stress on top of a temperature contour at the base of the conductive lid and the base of the cold thermal boundary layer. This program is used to analyze data from both the spherical and Cartesian ŒDIPUS models.

Scaling law analysis

The second program written for this study creates several plots of resultant values for each avis run, aiding in the development of scaling laws. One plot compares Rad,hot and Rad,cold, to Raconv (Figure 9), where

, and 


The DTTBL value is the temperature difference across the appropriate thermal boundary layer and h is the viscosity of the interior (constant-temperature) layer. The second compares DTh to DTTBL,cold and DTTBL,hot. (Figure 6). The third plots Tm and Raconv vs. avis (Figure 7), while the fourth compares Nu to Ra (Figure 8). The fifth and final plot compares observed values of qconv at the top and bottom of the layer with predicted values using Equation (17), which is based on the heat flux across the thermal boundary layers (Figure 11):


This equation is applicable to the current study because it has been determined that, for a given temperature-dependence of viscosity, Rad is constant over a wide range of Ra values13.

Figure 5. Mid-Depth non-dimensional temperature contour for spherical, Venus-based model with avis=20.
Figure 6. Plot of DTTBL vs. DTh for the limited data obtained from the Cartesian, Europa-based ŒDIPUS runs. Both hot and cold DTTBL are plotted. The green line represents DTTBL=1.43DTh, while the black line represents DTTBL = 1/2DT of the entire layer (the case when avis=0).
Figure 7. Non-dimensional Tm and Raconv vs. avis for Europa-based Cartesian simulations.
Figure 8. Logarithmic plot of Nuconv vs. Ra. The points marked with “+” represent data from the current study. The three thin lines represent plots of Equation (20) with different values of B. a is adjusted so that the lines fit the observed data point for avis=5. The thick line is a best-fit power regression, with Nu=1.27Ra0.1623. This line is not matched to fit a specific data point like the others.
Figure 9. Plot of Rad vs. Raconv for the hot and cold TBLs. Raconv decreases with increasing avis. Thus, avis increases from top to bottom of the figure.


Scaling law comparison

Much of the scaling law analysis conducted in this study involves the data from the Venus-based spherical model. As mentioned earlier, there is more data for these runs than for the Europa runs due to convergence issues in the Cartesian model. However, an analysis of spherical scaling laws from the current study appears in Ref. 15 and will not be repeated here. Instead, this section will focus on the brief analysis that was performed on the limited Europa data that was obtained.


As can be seen in Figure 6, as avis grows, the DTTBL,hot/DTh ratio seems to approach 1.43. Specifically, at avis=5, DTTBL,hot = 1.43DTh. This seems to support the previous conclusion; however, the experiments in this study were based on no-slip boundary conditions, which should make the thermal boundary layer more stable. This means it would require a higher DTTBL to form plumes than would be the case with free-slip conditions. While data from higher avis values are needed to confirm this theory, it does appear to be likely that the trendline relating to the bottom thermal boundary layer in Figure 6 will continue across the green DTTBL,hot = 1.43DTh line at values of avis higher than 5. This would correspond to a DTTBL,hot /DTh ratio >1.43, which is predicted for no-slip boundary conditions. Also in agreement with Ref. 8 is the finding that DTTBL,cold does not seem to show a strong proportionality to DTh, though again this finding is based on a very limited dataset.First, we compare DTh to DTTBL for both thermal boundary layers (Figure 6). In Ref. 8, it is observed that for a 2D fluid heated from below with highly temperature-dependent viscosity and free-slip boundary conditions, one obtains the following scaling law:

We next compare the mean convective-layer temperature and the convective Ra value with avis. The results are shown in Figure 7. For both Tm and Raconv, there are too few data points to begin to hypothesize about scaling laws, so this brief analysis simply provides qualitative reasoning for the upward and downward trends of these two values, respectively. The upward trend of Tm can be understood by considering that increasing avis corresponds to increasing the temperature-dependence of the viscosity (Equation 2). Since hbot is kept constant, hsurf and all interior viscosities increase (for a given temperature). Higher viscosities lead to lower velocities and diminished advective heat transfer ability. Thus, it takes longer to remove the heat from the lower boundary and the mean temperature increases. The behavior of Raconv is only slightly more complicated. There are some competing forces acting on the value as avis increases. First, the increased viscosity at the top of the convective layer leads to a thicker conductive lid, which decreases the convective layer thickness, the temperature difference across the convective layer, and consequently, Raconv. h(Tm) may increase or decrease because, although Tm will grow, the increased temperature-dependence will cause h(Tm) to decrease for a given Tm. Regardless, the third degree of dependence on convective layer thickness overcomes any possible decrease in h(Tm) and Raconv decreases with increasing avis.

Next, Nu is compared to Rasurf (Figure 8). Because we wish to evaluate the heat transfer in the convective layer (the entire layer except for the conductive lid), it is necessary to define a convective Nusselt number:


Using Eqs. (4,17), it can be seen that


and it follows that


For isoviscous fluids (in which there is no conductive lid), the amount of heat transferred is controlled by the instability of the thermal boundary layers. Thus Rad is constant for varying Ra, and, as observed in Ref. 13, the exponent b’ in Equation (20) is equal to 1/3. This proportionality has been confirmed with linear stability analysis and in previous numerical simulations of isoviscous fluids8,16,17. For temperature-dependent viscosities, the heat flux is controlled by the instability of the hot thermal boundary layer8. In the current study we are varying the temperature-dependence of viscosity and are consequently decreasing Rad,hot as we increase avis (Figure 9). Thus, because varying Ra implies varying avis, and thus varying Rad,hot, L, and Tc, the value of b obtained with the limited data is different than 1/3. Using a power fitting tool, we obtain coefficients for a and b of 1.270 and 0.1623, respectively (Figure 8), with an R2 value of 0.9956 for the regression. Those two values have 95% confidence ranges of 0.7032-1.837 and 0.1295-0.1951, again respectively. This range of b values falls between the fractions 1/5 and 1/8.

Part of the justification for this low value may have to do with the viscosity used in the Rayleigh number employed in the calculations. This viscosity, calculated in Equation (7), is the surface viscosity, which decreases at a more rapid rate than the mean convective-layer viscosity with increasing avis values. Indeed, when Nuconv is compared with Raconv in the current study, the value of b in Equation (21) surpasses 0.5 (Figure 10). In isoviscous fluids, where much of the work relating to this scaling law has been done, this distinction between surface viscosity and mean convective-layer viscosity does not exist. The no-slip boundaries employed in this study may also contribute to this low value. In previous no-slip-boundary experiments with isoviscous fluids, a value for b slightly less than 1/4 has been obtained8,18, rather than the predicted 1/3 that is confirmed in simulations of isoviscous fluids with free-slip boundary conditions.

Figure 10. Logarithmic plot of Nuconv vs. Raconv. The marked points correspond to simulation data, while the red line is a power regression of this data and the green line corresponds to the relation NuRa1/3.
Figure 11. Plot of observed vs. predicted heat flux for the top and bottom thermal boundary layers.

Finally, the study examines the prediction of qconv using Equation (17) for both the hot and cold thermal boundary layer with the observed heat flux at the top and bottom of the layer. The results are shown in Figure 11. The predicted qconv,botagrees well with the observed data, but the predicted qconv,topis not as well correlated with the observed surface heat flux. This is likely due to an inability to accurately define the boundary between the conductive lid and the upper thermal boundary layer (see Section II). This causes inaccuracies in dcold, DTTBL,cold, and, consequently, the predicted qconv,top.

In order to determine exactly how well the preliminary scaling laws developed using the Cartesian ŒDIPUS data correlate with previous studies, runs with larger values of avis will need to be performed. However, the initial data produced does show an expected ratio of DTTBL to DTh and a low but qualitatively understandable value for b in Equation (21). While these results support findings from earlier work, the results of the spherical ŒDIPUS code used in this study produce scaling laws more relevant to planetary interiors, and it will be useful to compare these results to the current Cartesian laws. For a detailed analysis of these results see Ref. 15.

Viscous stress

This study includes several comparisons between the viscous stress and yield strength properties of different Venus-based simulations. First, it analyzes the difference between these properties in runs with no-slip and with free-slip boundary conditions. Next, it looks at the difference when DT is varied in the Venus runs. Finally, it looks at the differences in shear stress/yield strength ratio when m is altered.

Five separate Venusian simulations are compared. Two runs incorporate a large DT (1140 K) and two a small DT (840 K). For each of these sets, one run involves no-slip boundaries and the other free-slip conditions. The fifth run uses the aforementioned parameters that give the lowest stress values (free-slip, small DT) and incorporates a value of 0.2 for the coefficient of friction in the mantle (Table 2).

Simulation DT, K Boundaries m τvisc,max, Mpa Radius, km Max Strain Rate, s-1 Radius, km smax, Mpa Radius, km ôyield at this depth, MPa
Rbig 1140 No-slip 0.05 299.2 6018.0 3.77E-015 3181.0 1033.3 6018.0 132.8
Rsmall 840 No-slip 0.05 213.7 6018.0 3.14E-015 3181.0 650.6 6018.0 132.8
Fbig 1140 Free-slip 0.05 274.1 5791.0 4.37E-015 3793.0 1064.5 5813.7 449.9
Fsmall 840 Free-slip 0.05 199.3 5768.3 3.50E-015 3589.6 678.0 5791.0 485.1
m0.2 840 Free-slip 0.2 199.3 5768.3 3.50E-015 3589.6 678.0 5791.0 1700.5

Table 2. Maximum viscous stress, strain rate, horizontal normal stress values, with corresponding radii and yield strength at the radius of maximum normal stress. In all runs, avis=20.

Effects of no-slip vs. free-slip boundary conditions (Figure 12)

Free-slip simulations lead to similar values for maximum viscous stress (τvisc,max) and horizontal normal stress (smax) as those obtained in no-slip-boundary simulations. More precisely, τvisc,max values differ by a factor <0.085 between the two boundary conditions cases, and smax values differ by a factor of <0.045. The main difference is in the location of the maximum stresses. No-slip boundaries lead to stresses that continue to increase towards the surface. The radius of the maximum stress in these cases should be equal to Venus’ radius but our model could not determine stresses in the outermost grid elements due to the manner in which we calculated the derivatives of the velocities. This is why values of 6018 km appear in Table 1, rather than the planetary radius of 6052 km that is used in our model. Still, it is clear that while no-slip boundaries cause the radius of the maximum stress values to be equal to the planetary radius, free-slip conditions cause this maximum value to occur >200km into the conductive lid. This is because stress values must approach zero at the surface when material is allowed to move freely along this surface. When the boundaries are no-slip, on the other hand, velocity must approach zero. This leads to larger ¶uxyz values and, consequently, greater stress near the surface. Strain rate values do not differ in the same fashion because velocities are low enough in the conductive lid to make their vertical gradient negligible without the viscosity multiplier (Equation 11).

Below the conductive lid, stresses in the free-slip and the no-slip simulations appear to be virtually equal. This is unsurprising, as the changing behavior of the surface should have less of an effect on material in the interior than it does on material close to the top. Strain rate, on the other hand, while similar above the cold thermal boundary layer, displays a much sharper increase in the hot thermal boundary layer in the no-slip boundary regime than it does in the free-slip regime. Again, this is due to the fact that velocity must reach zero at the boundaries in the no-slip case, thus leading to a larger vertical velocity gradient. The same difference in the stress values is not seen in Figure 12 because viscosity in the lower regions of the layer is low enough that τvisc and s are negligible.

Free-slip conditions also lead to a thinner conductive lid, by approximately 120-140 km. This is again due to the fact that velocities do not approach zero at the surface in the free-slip regime. Because of this, the lid’s boundary, which is defined as the point where horizontally averaged maximum uz is equal to 1/1000th of the global maximum, occurs at a higher point in the free-slip case.

Considering all of the above-mentioned factors, it is clear that at most points in the conductive lid, especially those near the surface where tyield is small, sfree£ sno-slip. Still, for m=0.05 and for either value of DT, s > tyield at some point in the lid for both the free-slip and no-slip regimes (Table 2). This implies lid fracture, which conflicts with Venus’ current stagnant lid regime.


Figure 12. Stress value profiles for spherical, Venus-based model with avis=20, DT=1140 K, no-slip and free-slip boundary conditions, respectively.

Effects of large versus small DT (Figure 13)

The two values tested in our study are 840 and 1140 K. The larger value led to horizontal normal stresses that were approximately 1.6 times greater than those generated with the smaller value. It also caused thinner conductive lids. Both of these phenomena are explained by the fact that a larger temperature difference leads to a greater vigor of convection. More vigorous convection is accompanied by both greater velocities, which lead to greater stresses, and a larger convective layer, which corresponds to a thinner conductive lid. Similar to the no-slip/free-slip comparison, both values chosen for DT cause sto exceed τyield at a point in the lid.


Figure 13. Stress value profiles for spherical, Venus-based model with avis=20, free-slip boundary conditions, and DT = 840 K, 1140 K, respectively.

Effects of large m (Figure 14)

The last case studied is that which involves a larger coefficient of friction. For this case, m = 0.2 was taken instead of 0.05. The larger value is equivalent to that obtained in laboratory experiments of Earth’s interior and is certainly within the (relatively unconstrained) realm of possibilities for Venus’ interior. In this simulation, the temperature and boundary parameters that lead to the lowest stress values in the lid (DT = 840 K, free-slip boundaries) were used. The only value that changes with a higher m is the rate of increase of tyield with depth. As can be seen in Figure 14, the yield strength increases rapidly enough in this case so that τyield ³ s at all points, indicating a stagnant lid. This particular run (DT = 840 K, free-slip boundaries, m = 0.2) was the only simulation in this study that produced conditions indicative of an unbroken lithosphere, which suggests that perhaps the friction coefficient of Venus’ interior is significantly higher than the 0.05 value given by numerical simulations for Earth’s interior. However, a wider parameter space must be examined in order to more accurately determine if this is the case.

Figure 14. Stress value profiles for spherical, Venus-based model with avis=20, low-stress parameters, and m=0.2.



This study presents an initial evaluation of data produced by 3D spherical and Cartesian numerical simulations of a fluid with temperature-dependent viscosity. Several MATLAB codes were created to aid in the analysis of this data. The parameters utilized in the simulations are derived from, and therefore allow us to apply the results to, the interiors of Europa (Cartesian model) and Venus (spherical model). From the preliminary Cartesian, Europa-based results, several conclusions are drawn, most of them qualitative, about scaling laws in this type of fluid and this type of domain. In addition, several case studies involving the spherical, Venus-based model are performed with varying parameters in order to determine the parameter space that produces results indicative of the observed stagnant lid.

It is noted that the previously determined scaling law relating temperature across the hot thermal boundary layer to viscous temperature scale (Equation 18) for 2D fluids applies in 3D Cartesian space (indeed it also applies in the spherical domain15). We also determine that mean convective-layer temperature increases while the convective Rayleigh number decreases with increased temperature-dependence, and these phenomena are qualitatively explained. We explain that the relationship between Nusselt number and Rayleigh number is not NuµRa1/3, as observed in isoviscous models, because a changing Rayleigh number in this study implies a changing temperature-dependence, which leads to a changing DTTBL (see Eqs. 20, 21). Next, it is explained why increasing avis values lead to smaller DTTBL values and smaller hot thermal boundary layer Rayleigh numbers. Finally, predicted and observed heat fluxes are compared, and it is determined that the predicted values for the bottom layer are much more accurate than those for the top. This is likely due to difficulties in determining the location of the boundary between the conductive lid and the upper thermal boundary layer.

With regard to the stresses, it is determined that lower DT values lead to significantly weaker viscous stresses and thicker conductive lids. Free-slip boundary conditions reduce the lid thickness and lead to lower depths for maximum stresses. This is important because yield strength increases linearly with depth. With a value for the friction coefficient of the Venusian mantle equal to that obtained in numerical simulations of Earth’s interior, the entire parameter space examined produces lid fractures. However, with low-stress parameters (low DT, free-slip boundaries) and a friction coefficient equal to that obtained in laboratory experiments with Earth-like materials, a situation is produced in which normal stress does not significantly exceed yield strength at any point. This friction value falls well within the range of possible values for Venus but suggests that the friction in the Venusian mantle may be relatively high compared with values predicted for Earth’s interior by numerical simulations.

The next step in stress analysis of the ŒDIPUS simulations should be to examine runs in which radiogenic heating is incorporated. These runs have been performed15, but the stress analysis is still pending. Additionally, a case in which internal heating is involved and the bottom viscosity is an order of magnitude smaller appears to produce a Venusian lid thickness within the estimated range of 100-200 km. A simulation with these parameters has yet to be run but should provide insight on the behavior of stresses in a domain that may be more closely related to the true interior of Venus.

Simulations of Europa have not as of yet been performed for avis values higher than 5, which makes it difficult to draw any significant results from the current data aside from preliminary scaling analysis. Simulations with higher temperature dependence of viscosity are currently being run, however, and more significant data should be available in the coming months. These results will provide insight on the observed fracturing of the icy Europan shell and on the transport of possible organic material from the satellite’s interior ocean to its surface.


I would like to thank Christophe Sotin and Sue Smrekar for their excellent mentorship, their willingness to always take the time to answer questions, and their ability to give me just the right amounts of direction and freedom that led to a rewarding internship experience. I would also like to thank my colleague Holly Taylor for sharing the burden of this project and for correcting my many coding errors along the way.

This research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, and was sponsored by The Undergraduate Student Research Program and the National Aeronautics and Space Administration.






1014 Pa s

1021 Pa s


5×10-5 K-1

3×10-5 K-1


1.2 m s-2

8.87 m s-2


920 kg m-3

3500 kg m-3


173 K

840, 1140 K


100 K

735 K


273 K

1175, 1575 K


1570 km19

6052 km


1520 km

3120 km


50 km

2932 km


1.3×10-6 m2 s-1

10-6 m2 s-1

Appendix A. Input parameters used in all Europa- and Venus-based simulations





Boundary Conditions

Tm, K

δhot, km

ΔTTBL,hot, K

δcold, km

ΔTTBL,cold, K

L, km

TcT0, K

DTh, K

















































































































Boundary Conditions





Predicted qconv,hot, W/m2

Predicted qconv,cold, W/m2

Observed qconv,hot, W/m2

Observed qconv,cold, W/m2



























































































Appendix B. Table of relevant results from Europa-based Cartesian simulations and Venus-based spherical simulations.


1Choblet, G., Cadek, O., Couturier, F., & Dumoulin, C., “ŒDIPUS: a new tool to study the dynamics of planetary interiors,” Geophysical Journal International, 170(1), 2007, pp. 9-30.

2Billings, S. E., & Kattenhorn, S. A. “The great thickness debate: Ice shell thickness models for Europa and comparisons with estimates based on flexure at ridges,” 2005.

3Hand, K.P., “Unmasking Europa: The Search for Life on Jupiter’s Ocean Moon,” Nature, 457, 2009, pp. 384-385.

4Dumoulin, C., Doin, M. P. & Fleitout, L., “Heat transport in stagnant lid convection with temperature- and pressure-dependent Newtonian or non-Newtonian Rheology,” Journal of Geophysical Research-Solid Earth, 104, 1999, pp. 12759-12777.

5Solomatov, V.S., and Moresi, L.N., “Scaling of time-dependent stagnant lid convection: Application to small-scale convection on Earth and other terrestrial planets,” Journal of Geophysical. Research, 105, 21, 2000, pp. 795-21,817.

6Steinberger, B. & Calderwood, A. R., “Models of large-scale viscous flow in the Earth’s mantle with constraints from mineral physics and surface observations,” Geophysical Journal International , 167, 2006, pp. 1461-1481.

7Sotin, C., Choblet, G., Head, J.W., Mocquet, A., & Tobie, G., “Thermal Evolution of Europa’s Icy Crust,” Universite de Nantes, Nantes, FR & Department of Geological Sciences, Brown University, Providence, RI, 2004.

8Deschamps, F. & Sotin, C., “Inversion of two-dimensional numerical convection experiments for a fluid with a strongly temperature-dependent viscosity,” Geophysical Journal International, 143, 2000, pp. 204-218.

9Valencia, D., O’Connell, R. J. & Sasselov, D. D., “Inevitability of plate tectonics on super-earths,” Astrophysical Journal, 670, 2007, pp. L45-L48.

10Jaeger, J.C., and N.G.W. Cook. Fundamentals  of Rock Mechanics, Chapman and Hall, New York, 1979, pp. 593.

11Beeman, M., Durham, W. B. & Kirby, S. H., “Friction of Ice,” Journal of Geophysical Research-Solid Earth and Planets, 93, 1988, pp. 7625-7633.

12Foster, A. & Nimmo, F., “Comparisons between the rift systems of East Africa, Earth and Beta Regio, Venus,” Earth and Planetary Science Letters, 143, 1996, pp. 183-195.

13Sotin, C. & Labrosse, S., “Three-dimensional thermal convection in an iso-viscous, infinite Prandtl number fluid heated from within and from below: applications to the transfer of heat through planetary mantles,” Physics of the Earth and Planetary Interiors, 112, 1999, pp. 171-190.

14Parmentier, E. M. & Sotin, C., “Three-dimensional numerical experiments on thermal convection in a very viscous fluid: Implications for the dynamics of a thermal boundary layer at high Rayleigh number,” Physics of Fluids, 12, 2000, pp. 609-617.

15Taylor, H., “Applications of 3D Numerical Simulations for Mantle Convection Within Venus,” Jet Propulsion Laboratory, National Aeronautics and Space Administration, Los Angeles, CA, 2009 (unpublished).

16Hansen, U. & Yuen, D. A., “High Rayleigh Number Regime of Temperature-Dependent Viscosity Convection and the Earth’s Early Thermal History,” Geophysical Research Letters, 20, 1993, pp. 2191-2194.

17Schubert, G. & Anderson, C. A., “Finite-Element Calculations of Very High Rayleigh Number Thermal Convection,” Geophysical Journal of the Royal Astronomical Society, 80, 1985, pp. 575-601.

18Frick, H., Busse, F. H. & Clever, R. M., “Steady 3-Dimensional Convection at High Prandtl Numbers,” Journal of Fluid Mechanics, 127, 1983,  pp. 141-153.