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.
Introduction
Several evolving methods of modeling planetary convection processes exist today. Until recently, none of these models has incorporated a temperature dependent viscosity in a threedimensional spherical space. In the current study, a finitevolume multigrid numerical simulation called ŒDIPUS^{1} 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 H_{2}O layer with estimates of its thickness ranging up to 120 km^{2}. 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 life^{3}, 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 Earthlike 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 noslip 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 Earthlike 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 temperaturedependence of viscosity. Finally, where previous scaling laws do not exist, the study provides qualitative explanations for the observed behavior of convectionrelated 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). T_{avg} in this layer is adiabatic, though in ŒDIPUS the effects of hydrostatic pressure are not incorporated and thus T_{avg} 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.
Methods
ŒDIPUS
Œ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 finitevolume grid with 128x128x128 elements to numerically integrate the NavierStokes equations. The Europan model has a 2x2x1 aspect ratio in the x,y,z directions, respectively, while the width in the xydirection for the spherical model is related to the curvature of the layer, defined using the ratio r_{c}/r_{p}. 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 r_{c}/r_{p} value and did not run successfully. We used multiple runs to test freeslip and noslip boundary conditions for both the spherical and Cartesian models. Ideally, one might want to set a noslip bottom boundary (representing the mantle/core interface) and a freeslip 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.
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:
(1).
Unfortunately, this equation leads to large viscosity differences on the scale of 10^{15} Pa s, which cause divergences in numerical simulations. Therefore, ŒDIPUS utilizes the following approximation:
(2).
where h_{0} is the viscosity at T=T_{0}, and avis is a constant that is related to H through the viscous temperature scale, defined in the following equations:
(3).
The nonNewtonian effects on viscosity have previously been investigated, but their effect is much less than that of temperature^{4,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 14. 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 nondimensional value that, in Cartesian coordinates, represents the ratio of convective heat flux (q_{conv}, the sum of advective and diffusive fluxes) to what the conductive heat flux would be if there were no advection {q_{cond}=k(T_{1}T_{0})/d}. It is defined as
(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 nondimensional value that can be compared to heat flux. In this case it is defined as
(6).
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 r^{2} difference in surface area from the coremantle boundary (CMB) to the surface, where r=r_{c}/r_{p}, it is necessary to multiply Nu_{top} by r^{2} before comparing it to Nu_{bot}. 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:
(7).
For each value of avis, h_{0} is not defined directly. Rather h_{1}, the bottom viscosity, is defined for all avis values. For Venus, h_{1}= 10^{21} Pa s. This value is based loosely on an Earth analog^{6}. The Europan runs are performed with h_{1} = 10^{14}, based on findings in Ref. 7. Equation (2) is then inverted to determine h_{0}. Ra varies when avis is changed due to the change in h_{0}.
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 nondimensionalized using the factors in Table 1. To obtain velocities in the center of the elements (so that u_{x}, u_{y}, and u_{z} 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 (q_{conv}/k) extended toward the interior intersects with T_{m}^{8}. 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 temperaturedependence 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/1000^{th} 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.
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,
(8).
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 t_{visc} should not be significantly affected.
Both t_{x,z} and t_{y,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: ¶u_{x}/¶z, ¶u_{z}/¶x, ¶u_{y}/¶z, and ¶u_{z}/¶y. These derivatives are calculated as follows:
(9,10).
The ¶u_{y}/¶z, and ¶u_{z}/¶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 (t_{visc}) for each point by finding the sum of the two strain and the two stress vectors:
(11).
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):
(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 T_{avg} 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 midlayer 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:
(13).
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,
(14),
where z is depth. Values of C_{0} range from 25 to 100 MPa for silicates. For ice, according to laboratory experiments^{11}, C_{0}=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 C_{0} value for silicates of 80 MPa and values of m from 0.05 and 0.2 are used. The C_{0} 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 Ra_{d,hot} and Ra_{d,cold}, to Ra_{conv} (Figure 9), where
(15,16).
The DT_{TBL} value is the temperature difference across the appropriate thermal boundary layer and h is the viscosity of the interior (constanttemperature) layer. The second compares DT_{h} to DT_{TBL,cold} and DT_{TBL,hot}. (Figure 6). The third plots T_{m} and Ra_{conv} vs. avis (Figure 7), while the fourth compares Nu to Ra (Figure 8). The fifth and final plot compares observed values of q_{conv} 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):
(17).
This equation is applicable to the current study because it has been determined that, for a given temperaturedependence of viscosity, Ra_{d} is constant over a wide range of Ra values^{13}.
Results
Scaling law comparison
Much of the scaling law analysis conducted in this study involves the data from the Venusbased 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.
(18).
As can be seen in Figure 6, as avis grows, the DT_{TBL,hot}/DT_{h} ratio seems to approach 1.43. Specifically, at avis=5, DT_{TBL,hot} = 1.43DT_{h}. This seems to support the previous conclusion; however, the experiments in this study were based on noslip boundary conditions, which should make the thermal boundary layer more stable. This means it would require a higher DT_{TBL} to form plumes than would be the case with freeslip 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 DT_{TBL,hot} = 1.43DT_{h} line at values of avis higher than 5. This would correspond to a DT_{TBL,hot }/DT_{h} ratio >1.43, which is predicted for noslip boundary conditions. Also in agreement with Ref. 8 is the finding that DT_{TBL,cold} does not seem to show a strong proportionality to DT_{h}, though again this finding is based on a very limited dataset.First, we compare DT_{h} to DT_{TBL} for both thermal boundary layers (Figure 6). In Ref. 8, it is observed that for a 2D fluid heated from below with highly temperaturedependent viscosity and freeslip boundary conditions, one obtains the following scaling law:
We next compare the mean convectivelayer temperature and the convective Ra value with avis. The results are shown in Figure 7. For both T_{m} and Ra_{conv}, 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 T_{m} can be understood by considering that increasing avis corresponds to increasing the temperaturedependence of the viscosity (Equation 2). Since h_{bot }is kept constant, h_{surf} 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 Ra_{conv} 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, Ra_{conv}. h(T_{m}) may increase or decrease because, although T_{m} will grow, the increased temperaturedependence will cause h(T_{m}) to decrease for a given T_{m}. Regardless, the third degree of dependence on convective layer thickness overcomes any possible decrease in h(T_{m}) and Ra_{conv} decreases with increasing avis.
Next, Nu is compared to Ra_{surf} (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:
(19).
Using Eqs. (4,17), it can be seen that
(20).
and it follows that
(21).
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 Ra_{d} 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 fluids^{8,16,17}. For temperaturedependent viscosities, the heat flux is controlled by the instability of the hot thermal boundary layer^{8}. In the current study we are varying the temperaturedependence of viscosity and are consequently decreasing Ra_{d,hot} as we increase avis (Figure 9). Thus, because varying Ra implies varying avis, and thus varying Ra_{d,hot}, L, and T_{c}, 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 R^{2} value of 0.9956 for the regression. Those two values have 95% confidence ranges of 0.70321.837 and 0.12950.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 convectivelayer viscosity with increasing avis values. Indeed, when Nu_{conv} is compared with Ra_{conv} 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 convectivelayer viscosity does not exist. The noslip boundaries employed in this study may also contribute to this low value. In previous noslipboundary experiments with isoviscous fluids, a value for b slightly less than 1/4 has been obtained^{8,18}, rather than the predicted 1/3 that is confirmed in simulations of isoviscous fluids with freeslip boundary conditions.
Finally, the study examines the prediction of q_{conv} 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 q_{conv,bot}agrees well with the observed data, but the predicted q_{conv,top}is 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 d_{cold}, DT_{TBL,cold}, and, consequently, the predicted q_{conv,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 DT_{TBL} to DT_{h} 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 Venusbased simulations. First, it analyzes the difference between these properties in runs with noslip and with freeslip 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 noslip boundaries and the other freeslip conditions. The fifth run uses the aforementioned parameters that give the lowest stress values (freeslip, 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  s_{max}, Mpa  Radius, km  ô_{yield} at this depth, MPa 
R_{big}  1140  Noslip  0.05  299.2  6018.0  3.77E015  3181.0  1033.3  6018.0  132.8 
R_{small}  840  Noslip  0.05  213.7  6018.0  3.14E015  3181.0  650.6  6018.0  132.8 
F_{big}  1140  Freeslip  0.05  274.1  5791.0  4.37E015  3793.0  1064.5  5813.7  449.9 
F_{small}  840  Freeslip  0.05  199.3  5768.3  3.50E015  3589.6  678.0  5791.0  485.1 
m_{0.2}  840  Freeslip  0.2  199.3  5768.3  3.50E015  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 noslip vs. freeslip boundary conditions (Figure 12)
Freeslip simulations lead to similar values for maximum viscous stress (τ_{visc,max}) and horizontal normal stress (s_{max}) as those obtained in noslipboundary simulations. More precisely, τ_{visc,max} values differ by a factor <0.085 between the two boundary conditions cases, and s_{max} values differ by a factor of <0.045. The main difference is in the location of the maximum stresses. Noslip 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 noslip boundaries cause the radius of the maximum stress values to be equal to the planetary radius, freeslip 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 noslip, on the other hand, velocity must approach zero. This leads to larger ¶u_{xy}/¶z 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 freeslip and the noslip 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 noslip boundary regime than it does in the freeslip regime. Again, this is due to the fact that velocity must reach zero at the boundaries in the noslip 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.
Freeslip conditions also lead to a thinner conductive lid, by approximately 120140 km. This is again due to the fact that velocities do not approach zero at the surface in the freeslip regime. Because of this, the lid’s boundary, which is defined as the point where horizontally averaged maximum u_{z} is equal to 1/1000^{th} of the global maximum, occurs at a higher point in the freeslip case.
Considering all of the abovementioned factors, it is clear that at most points in the conductive lid, especially those near the surface where t_{yield} is small, s_{free}£ s_{noslip}. Still, for m=0.05 and for either value of DT, s > t_{yield} at some point in the lid for both the freeslip and noslip regimes (Table 2). This implies lid fracture, which conflicts with Venus’ current stagnant lid regime.
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 noslip/freeslip comparison, both values chosen for DT cause sto exceed τ_{yield} at a point in the lid.
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, freeslip boundaries) were used. The only value that changes with a higher m is the rate of increase of t_{yield} 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, freeslip 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.
Conclusion
This study presents an initial evaluation of data produced by 3D spherical and Cartesian numerical simulations of a fluid with temperaturedependent 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, Europabased 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, Venusbased 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 domain^{15}). We also determine that mean convectivelayer temperature increases while the convective Rayleigh number decreases with increased temperaturedependence, and these phenomena are qualitatively explained. We explain that the relationship between Nusselt number and Rayleigh number is not NuµRa^{1/3}, as observed in isoviscous models, because a changing Rayleigh number in this study implies a changing temperaturedependence, which leads to a changing DT_{TBL} (see Eqs. 20, 21). Next, it is explained why increasing avis values lead to smaller DT_{TBL} 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. Freeslip 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 lowstress parameters (low DT, freeslip boundaries) and a friction coefficient equal to that obtained in laboratory experiments with Earthlike 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 performed^{15}, 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 100200 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.
Acknowledgements
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.
Appendix
Parameter 
Europa 
Venus 
h_{bottom} 
10^{14} Pa s 
10^{21} Pa s 
a 
5×10^{5} K^{1} 
3×10^{5} K^{1} 
g 
1.2 m s^{2} 
8.87 m s^{2} 
r 
920 kg m^{3} 
3500 kg m^{3} 
DT 
173 K 
840, 1140 K 
T_{0} 
100 K 
735 K 
T_{1} 
273 K 
1175, 1575 K 
r_{p} 
1570 km^{19} 
6052 km 
r_{c} 
1520 km 
3120 km 
d 
50 km 
2932 km 
k 
1.3×10^{6} m^{2} s^{1} 
10^{6} m^{2} s^{1} 
Appendix A. Input parameters used in all Europa and Venusbased simulations
Planet 
avis 
DT 
Boundary Conditions 
T_{m}, K 
δ_{hot}, km 
ΔT_{TBL,hot}, K 
δ_{cold}, km 
ΔT_{TBL,cold}, K 
L, km 
T_{c} – T_{0}, K 
DT_{h}, K 
Nu_{bot} 
Nu_{top} 
20 
840 
NS 
1438.8 
129.01 
136.20 
664.04 
341.31 
990.92 
362.50 
42 
3.86 
1.07 

VENUS 
20 
840 
FS 
1463.3 
85.26 
111.73 
620.52 
364.35 
850.81 
363.91 
42 
4.79 
1.28 
20 
1140 
NS 
1685.2 
118.35 
189.81 
641.95 
481.18 
873.10 
469.01 
57 
4.31 
1.18 

20 
1140 
FS 
1727.8 
75.23 
147.16 
612.50 
525.04 
750.20 
467.80 
57 
5.30 
1.40 

1 
173 
NS 
195.38 
1.51 
77.62 
1.63 
89.11 
0.12 
6.26 
173 
14.91 
15.06 

EUROPA 
3 
173 
NS 
210.37 
1.49 
62.63 
2.20 
92.78 
0.42 
17.59 
58 
12.21 
12.10 
4 
173 
NS 
217.41 
1.54 
55.59 
2.64 
96.40 
0.59 
21.01 
43 
10.45 
10.39 

5 
173 
NS 
223.46 
1.60 
49.54 
3.09 
96.86 
0.87 
26.60 
35 
8.95 
8.89 
Planet 
avis 
DT 
Boundary Conditions 
Ra 
Ra_{conv} 
Ra_{δ,hot} 
Ra_{δ,cold} 
Predicted q_{conv,hot}, W/m^{2} 
Predicted q_{conv,cold}, W/m^{2} 
Observed q_{conv,hot}, W/m^{2} 
Observed q_{conv,cold}, W/m^{2} 
20 
840 
NS 
4.06E02 
1.27E+05 
10.6 
3635.1 
1.06E03 
5.14E04 
1.10E03 
3.06E04 

VENUS 
20 
840 
FS 
4.06E02 
2.80E+05 
4.5 
5669.7 
1.31E03 
5.87E04 
1.37E03 
3.67E04 
20 
1140 
NS 
5.52E02 
1.95E+05 
10.5 
4243.6 
1.60E03 
7.50E04 
1.68E03 
4.60E04 

20 
1140 
FS 
5.52E02 
4.92E+05 
4.4 
8499.8 
1.96E03 
8.57E04 
2.06E03 
5.45E04 

1 
173 
NS 
3.38E+06 
5.61E+06 
73.0 
103.9 
5.13E02 
5.48E02 
5.16E02 
5.21E02 

EUROPA 
3 
173 
NS 
4.57E+05 
2.71E+06 
29.8 
141.3 
4.20E02 
4.22E02 
4.22E02 
4.19E02 
4 
173 
NS 
1.68E+05 
2.15E+06 
24.0 
209.5 
3.60E02 
3.64E02 
3.62E02 
3.59E02 

5 
173 
NS 
6.19E+04 
1.76E+06 
20.4 
290.6 
3.10E02 
3.13E02 
3.10E02 
3.08E02 
Appendix B. Table of relevant results from Europabased Cartesian simulations and Venusbased spherical simulations.
References
^{1}Choblet, 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. 930.
^{2}Billings, 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.
^{3}Hand, K.P., “Unmasking Europa: The Search for Life on Jupiter’s Ocean Moon,” Nature, 457, 2009, pp. 384385.
^{4}Dumoulin, C., Doin, M. P. & Fleitout, L., “Heat transport in stagnant lid convection with temperature and pressuredependent Newtonian or nonNewtonian Rheology,” Journal of Geophysical ResearchSolid Earth, 104, 1999, pp. 1275912777.
^{5}Solomatov, V.S., and Moresi, L.N., “Scaling of timedependent stagnant lid convection: Application to smallscale convection on Earth and other terrestrial planets,” Journal of Geophysical. Research, 105, 21, 2000, pp. 79521,817.
^{6}Steinberger, B. & Calderwood, A. R., “Models of largescale viscous flow in the Earth’s mantle with constraints from mineral physics and surface observations,” Geophysical Journal International , 167, 2006, pp. 14611481.
^{7}Sotin, 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.
^{8}Deschamps, F. & Sotin, C., “Inversion of twodimensional numerical convection experiments for a fluid with a strongly temperaturedependent viscosity,” Geophysical Journal International, 143, 2000, pp. 204218.
^{9}Valencia, D., O’Connell, R. J. & Sasselov, D. D., “Inevitability of plate tectonics on superearths,” Astrophysical Journal, 670, 2007, pp. L45L48.
^{10}Jaeger, J.C., and N.G.W. Cook. Fundamentals of Rock Mechanics, Chapman and Hall, New York, 1979, pp. 593.
^{11}Beeman, M., Durham, W. B. & Kirby, S. H., “Friction of Ice,” Journal of Geophysical ResearchSolid Earth and Planets, 93, 1988, pp. 76257633.
^{12}Foster, A. & Nimmo, F., “Comparisons between the rift systems of East Africa, Earth and Beta Regio, Venus,” Earth and Planetary Science Letters, 143, 1996, pp. 183195.
^{13}Sotin, C. & Labrosse, S., “Threedimensional thermal convection in an isoviscous, 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. 171190.
^{14}Parmentier, E. M. & Sotin, C., “Threedimensional 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. 609617.
^{15}Taylor, H., “Applications of 3D Numerical Simulations for Mantle Convection Within Venus,” Jet Propulsion Laboratory, National Aeronautics and Space Administration, Los Angeles, CA, 2009 (unpublished).
^{16}Hansen, U. & Yuen, D. A., “High Rayleigh Number Regime of TemperatureDependent Viscosity Convection and the Earth’s Early Thermal History,” Geophysical Research Letters, 20, 1993, pp. 21912194.
^{17}Schubert, G. & Anderson, C. A., “FiniteElement Calculations of Very High Rayleigh Number Thermal Convection,” Geophysical Journal of the Royal Astronomical Society, 80, 1985, pp. 575601.
^{18}Frick, H., Busse, F. H. & Clever, R. M., “Steady 3Dimensional Convection at High Prandtl Numbers,” Journal of Fluid Mechanics, 127, 1983, pp. 141153.