erfa | Fachverein Erdwissenschaften

Numerical modeling of the interaction
between a mantle plume and a moving
oceanic lithosphere
Valentin Gischig
WS 2005/06
Contents
1 Abstract
2
2 Introduction
2
3 Model description
3.1 Numerical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Physical parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
3
6
7
4 Results
9
4.1 Buoyancy flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
4.2 Background viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4.3 Arrhenius activation energy . . . . . . . . . . . . . . . . . . . . . . . . . . 13
5 Discussion and conclusion
15
Literature
17
1
1 Abstract
The topographic swell above a mantle plume below an oceanic lithosphere is partly due
to plume underplating and partly due to heating up and thinning of the lithosphere by
the plume. The relationship between lithospheric erosion by plume heating and mantle
and plume properties is examined with a numerical model. The empirical relationships,
so-called scaling laws, are ∆ ∝ F 0.0912 , ∆ ∝ η0−0.2029 and ∆ ∝ E −0.4198 , where ∆ is the
maximum distance (in km) between the 1250◦ C-isotherms of a thinned lithosphere and
an undisturbed lithosphere, F the buoyancy flux (in kg/s), η0 the background viscosity
(in P a · s), and E the Arrhenius activation energy (in kJ/mol). Only the first two scaling
laws seem to be reasonable with respect to physical expectations for the vertical thinning.
Better results could be expected if the same calculations were done for the whole volume
of eroded lithosphere instead of the maximal vertical thinning.
2 Introduction
For more than 40 years now the pattern of volcanism of the Hawaiian Islands has been
thought to be due to a mantle plume meeting the lithosphere underneath Hawaii. Another
phenomena observed around this area is a large region of shallow sea floor, the so-called
“swell” of the lithosphere with around 1km height and 1000km width. Several models for
the plate-plume-interaction have been proposed to explain this swell [Ribe, 2004] (Figure
1). A first model (Figure 1a) suggests that the uplift of the lithosphere can be explained
with the heat of the plume that heats and thins the lithosphere which is thus pushed
upwards, because of the lower density of the warmer lithosphere material [Detrick and
Crough, 1978]. But calculations showed that the time scale for heating the lithosphere
by diffusion is much larger than observed. A second model (Figure 1b) explains the uplift
of the lithosphere with a mantle plume underplating the lithosphere and being dragged
with the moving plate in the direction of the plate motion [Olson, 1990]. The buoyancy
of this pancake-shaped plume, again less dense than the lithosphere and the surrounding
mantle material, applies a upward force on the lithosphere. In this model there is no
heat transferred to the lithosphere and its thickness stays constant. Seismic studies of
the Hawaiian swell by Li et al. [2004] suggest a certain thinning of the lithosphere in this
area. Li et al. [2004] assume that the best model is a mixture of the first two models
(Figure 1c). An explanation for the heat transfer faster than by diffusion is small scale
convection in the mantle plume that could transfer heat fast enough to the overlying
plate to make some thinning of the plate possible [Moore et al., 1998].
Most previous numerical calculations have considered only the underplating model (no
heat transfer by diffusion) up to now. But some numerical models, that came up in the
last years, showed some degree of plate-thinning. It’s obvious that the dimensions of the
2
Figure 1: Three different models for the formation of the lithospheric swell produced by
a mantle plume. Explanations in the text. From [Ribe, 2004].
swell are dependent on several mantle parameters such as viscosity, plume dimensions,
buoyancy flux, excess temperature of the plume, Arrhenius activation energy, etc. Empirical relationships between these parameters and the swell dimensions, so-called scaling
laws, help better understanding the physical processes responsible for the swell. This
has already been done for width and height of the swell and the buoyancy flux of the
plume [van Hunen and Zhong, 2003]. The aim of this semester work is to find out about
the (empirical) relationship between the thinning of the lithosphere and buoyancy flux,
background viscosity and Arrhenius activation energy. This is done by using a numerical
model of a very general situation of plume-plate-interaction. Finding scaling laws for all
the other important parameters would go beyond the scope of this semesterwork. Since
for many plumes on Earth the plume head is no longer on top of the plume tail and is
very often already subducted, just the thinning by the plume tail is studied here.
The structure of this semester work is as follows: In the next section the setup of the numerical model (code, equations, geometry, boundary conditions, etc.), the post-processing
steps and the meaning of the model parameters, that have been varied, will be described.
In section 4 the results of the calculations for each parameter will be presented separately.
These results will be put into a physical context and are compared with previous studies
from literature in section 5.
3 Model description
3.1 Numerical model
To simulate the plate-plume-interaction the parallel Cartesian 3-D code “citcom” is used.
This numerical finite element model, written in C-language, is adapted for plume-plate
interaction as in [Zhong and Watts, 2002] using the setup from [Ribe and Christensen,
3
1994]. The dimensional equations solved by the code are:
Continuity equation:
∇·u=0
(1)
−∇P + ∇ · η ∇u + ∇T u + ∆ρg = 0
(2)
Stokes equation:
Energy equation:
ρcp
Parameter
u
P
η
ρ
∆ρ
g
cp
T
t
k
∂T
+ u · ∇T
∂t
= ∇ · k∇T
Description
Velocity vector
Pressure
Viscosity
Density
Density defect
Gravity constant
Specific heat
Absolute temperature
Time
Thermal conductivity
(3)
Units/values
m/s
Pa
Pa · s
kg/m3
kg/m3
m/s2
1250kJ/kg
K
s
W/mK
Table 1: Explanation of the parameters as used in equations 1 to 3.
For the mantle rheology diffusion creep is assumed using the linear Arrhenius viscosity
(see below). The model consists of a box of 400km depth, 1600km width and 3200km
length and and has a 5km element size. The boundary conditions are as follows(Figure
2):
• At the surface: T=0◦ C, the velocity v is the value estimated for the Pacific plate
around Hawaii relative to the mantle below (8.6cm/yr) and points in the x-direction.
• At the bottom: T=1350◦ C, v=0.
• At the y-z-planes: Temperature profile is the error function (cooling halfspace
temperature profile, [Turcotte and Schubert, 2002]), velocity grows linearly from 0
at the bottom to the plate velocity at the surface, such that the shear stress σxz is
contstant.
• At the x-z-planes: Symmetry and free slip boundary conditions: vy = 0, σxy =
σzy = 0 and ∂T /∂y = 0.
4
Figure 2: Sketch of the geometry and the boundary conditions of the models as used for
this study.
Since the model is symmetric to the y=0-plane the calculations have the be run only for
one half of the plume. The initial model consists of a plume positioned at the x-axis at
◦
800km as a hot circular region with a excess temperature ∆T
p = 300 C in the middle
2
2
of the plume and a temperature of ∆T = ∆Tp exp −r /R at the distance r from the
plume center. The cross section parallel to the plate motion through this initial situation
is shown in figure 3.
Figure 3: Cross section at y=0 of the initial model. The black line is the 1250◦ C-isotherm,
which is used as the lithosphere-asthenosphere boundary in this work. The hot
region at the bottom is the initial plume. Note that the scaling of the x- and
the z-axis is not equal to make features better visible.
Starting from this initial plume the plume head rises up, meets the lithosphere and
flows along the lithosphere, until just the plume tail remains underneath the plate. The
5
numerical model is run for 5000 to 15000 time steps to reach this steady state situation.
After that the shape of the plume does not change significantly anymore and the thinning
of the lithosphere stays constant. The model time between each time step is not constant
and depends on the velocity of the mantle material. This means that the time between
two time steps is adjusted such that material does not move more than half an element
size in order to make the numerical timestepping scheme stable. For example for high
mantle viscosities, 5000 time steps give a time span of around 50 million years. The
buoyancy flux (defined below) has to be constant during the whole run of a model in
these calculations. Since temperature and the temperature dependent viscosity change
in the material surrounding the plume during the uprising of the plume the flux would
also change during the run of one model. To make sure that the flux stays constant
in time, the numerical code measures the flux after a certain number of time steps and
changes the radius of the plume automatically to fix the flux to the desired value. The
calculations were made on the Linux-cluster named “stella”, mostly on two processors,
and took 10 to 20 hours of computation time each, depending on the number of time
steps.
3.2 Post-processing
To measure the thinning of the lithosphere, the total volume of lithosphere that has been
eroded away (i.e. heated up by small scale convection), should be considered. Since
this is difficult to measure, alternatively just the maximum thinning of the lithosphere is
measured, which is the maximum distance between the eroded lithosphere and a lithosphere that has not been altered by a plume. For measuring the maximum thinning
of the lithosphere it is enough to look at a 2-D cross-section of the whole 3-D model
where the thinning through the point where thinning is maximal. This is the vertical
section through the x-axis and containing the plume center and the direction of the
velocity of the plume. For comparison a reference model with no plume at all is also
calculated (Figure 4). The lower boundary of the lithosphere is defined thermally and
is the 1250◦ C-isotherm. The maximum thinning is then plotted against the parameters
that are varied in order to find a relationship between these values. The relationships are
set up using linear regression. Correlation coefficient, the rms-value of the residuals and
distribution of the residuals are used to decide how good the fitted functions are. This
post-processing is done with “matlab”.
6
Figure 4: Cross section at y=0 for a plume in its steady state (top) and an undisturbed
mantle (below), i.e. with no plume at all. The lithosphere thickens as it moves
from left to right and gets older.
3.3 Physical parameters
Buoyancy flux: The buoyancy flux is the amount of density defect that rises up through
the plume in a certain time interval:
Z ~ · ∆ρdS
B=
−K
(4)
S
~ is the velocity of the material penetrating through a certain surface increment,
where K
∆ρ the density difference between the hot plume and the background mantle and S the
horizontal surface through which the flux is measured. B depends on several parameters
as viscosity of the plume material, excess temperature and the radius of the plume.
Numerical models are run for half flux values between 1000kg/s and 3000kg/s (because
only one half of the plume is modeled). For the Hawaiian plume a flux of 6000kg/s has
been estimated [Davies, 1988, Sleep, 1990], which is a half flux of 3000kg/s.
Background viscosity: For the mantle a linear viscous Arrhenius rheology for diffusion
creep is assumed, which means that the mantle viscosity is
E 1
1
η = η0 exp
−
(5)
R T
Tm
where T is the absolute temperature, Tm is the absolute temperature of the mantle, R is
the gas constant and E is the background activation energy. The background viscosity
η0 is referred to as the viscosity of the mantle material and is thus a realistic value for
the average asthenosphere. This parameter is changed here between 1019 and 1020 P a · s.
7
These values are thought to be appropriate for the mantle although for the asthenosphere
even lower values are suggested.
Arrhenius activation energy: According to the exponential function in equation 5,
mantle viscosity changes slowly with T at temperatures around Tm and increases very
fast for lower temperatures. This means that at the transition to the lithosphere the
viscosity increases very abruptly and material gets stiff. How fast the stiffness of the
lithosphere increases at the bottom depends on the Arrhenius activation energy E. Since
lithosphere material can only be eroded away if it is weak (i.e. has a low viscosity),
the activation energy E determines how much of the lithosphere is enough mobile to
be eroded and has a strong influence on lithosphere thinning. For low values of E the
stiffness of the lithosphere increases slower and it becomes easier to thin the lithosphere.
Mantle rheology is usually thought to be best described by dislocation creep viscosity.
For dislocation creep an activation energy of 350kJ/mol is estimated. But according to
[Christensen, 1984] dislocation energy can be simulated with a diffusion creep with 2
to 3 times lower Arrhenius activation energy, which is much easier to implement in a
numerical model and takes much less computation time. Thus in this models diffusion
creep with Arrhenius activation energies between 120kJ/mol and 300kJ/mol have been
used. Values higher than 200kJ/mol can be considered as rather high for the mantle in
this respect.
8
4 Results
4.1 Buoyancy flux
The numerical model is run for buoyancy half fluxes between 1000kg/s and 3000kg/s in
steps of 500kg/s. The cross section containing the symmetry axis (i.e. at y=0) with
the plume in its steady state is shown in figure 5 for all values. Especially for higher
fluxes one can see at the left side of the plume head areas of cold material that seem
to be dragged down to the mantle and can be considered as small scale convection. As
expected higher fluxes thin the lithosphere more since there is more heat brought to the
lithosphere.
In order to decide whether the model reached its steady state for all 500 time steps the
maximum thinning of the lithosphere is plotted (Figure 6). For all different values the
maximum thinning first increases very much until it reaches a maximum value. This is
due to the head of the plume reaching the lithosphere. The lithosphere then increases
its thickness again and reaches a steady state. In this case no longer the plume head but
the tail of the plume causes the lithosphere thinning. This is the situation we meet at
most plumes of the earth. Steady state is reached for all buoyancy fluxes within 5000
time steps. Since the maximum thinning stays constant at steady state, the values at
the 5000th time step is used for setting up a scaling law. These data points and fitted
curve is shown in figure 7.
Here the data points are tried to be fitted with a power-law, which means that first the
logarithm has been calculated for all data points to which then a linear regression is
applied. It turns out that the logarithmic data is best fitted with the following linear
function:
log(∆) = 0.0912 · log(F ) + 1.3307
(6)
where ∆ is the maximum thinning in km and F the buoyancy flux in kg/s. The regression has a correlation coefficient of 0.9874 and a rms-value of the residuals of 0.0025,
which means that the following empirical power-law relation can explain the relationship
between lithosphere thinning and the flux of the mantle plume quite well:
∆ = 21.4 · F 0.0912
(7)
Although the residuals in figure 7 are very small, they have a certain dependency on the
flux (i.e. not normally distributed), which means that the scaling law can not explain all
of the relationship between thinning and flux.
9
Figure 5: Cross section at y=0 for different values of buoyancy flux. In order to make
details more visible the scale of the x- and z-axis are not the same. The plots
are vertically enlarged.
Figure 6: Maximum lithosphere thinning ∆ at all 500 time steps. The steady state is
reached as soon the line becomes flat.
10
Figure 7: Buoyancy flux vs. lithosphere thinning. The black crosses are the data points
coming out of the numerical calculations, the red line is a power-law fitted into
this data.
4.2 Background viscosity
In figure 8 again the cross section through the plume in its steady state is shown, this
time for different values of the background viscosity. Small scale convection becomes
stronger and more cold material is dragged down into the mantle with lower viscosities.
To reach a steady state plume this time up to 15000 time steps have been necessary (Figure 9), because the mantle material is more mobile at lower viscosities and thus moves
more vigorously before it reaches a steady state. For very low viscosities the steady state
is slightly oscillating, so that not the maximum thinning value of the last time step but
the average of the last few time steps are used in further calculations.
The scaling law for viscosity is again a power-law (Figure 10). The data is fitted using
linear regression in the same way as in the previous paragraph with a correlation coefficient of -0.9989 and a rms-value of the residuals of 0.0035. The result of the linear
regression is
log(∆) = −0.2029 · log(η0 ) + 5.7031
(8)
what gives a power-law relationship between maximum thinning and background viscosity
of the mantle:
∆ = 5.05 · 105 · (η0 )−0.2029
(9)
11
Figure 8: Cross section at y=0 for different values of background viscosity. Again the
scale is not equal vertically and horizontally.
Figure 9: Maximum lithosphere thinning ∆ at all 500 time steps. For low viscosities an
oscillating steady state is reached.
12
Figure 10: Background viscosity vs. lithosphere thinning.
4.3 Arrhenius activation energy
Figure 11 shows the result of calculations with different values of the Arrhenius activation
energy. For values higher than 200kJ/mol it was necessary to run the model for up to
10000 time steps (Figure 12). Since the material moves faster at higher values of E in the
high-temperature region of the plume (because of very low viscosity inside the plume),
time steps are getting smaller and there are more timesteps necessary to reach a steady
state.
The data coming out of these calculations are again best fitted with a power-law (Figure
13), which has a correlation coefficient of -0.9962 and a rms-value of 0.0062:
log(∆) = −0.4198 · log(E) + 2.5150
(10)
whereas the residuals are again not normally distributed. This gives
∆ = 327 · E −0.4198
13
(11)
Figure 11: Cross section at y=0 for different values of Arrhenius activation energy. (Vertical axis enlarged)
Figure 12: Maximum lithosphere thinning ∆ at all 500 time steps.
14
Figure 13: Activation energy vs. lithosphere thinning.
5 Discussion and conclusion
Numerical results and the obtained scaling laws are compared with other studies and
tried to be set into a physical background here. Looking at figures 7, 10 and 13 we
recognize that in general lithosphere is thinned by an uprising plume within a range
of 30km at very high activation energies and 70km for very low background viscosities.
The results of seismic studies made by [Li et al., 2004] showed that the undisturbed
90Ma to 100Ma year old pacific lithosphere being 100 to 110km thick is thinned out
by the Hawaiian plume to a thickness of 50 to 60km underneath the island Kauai. Although the lithosphere-asthenosphere-boundary is defined by a seismic discontinuity in
that study and not by the 1250◦ C-isotherm, as it is done here, the thinning of the Hawaiian lithosphere is comparable to the results of these numerical calculations. This can be
considered as evidence, that model 3 as explained in the introduction (Figure 1c) is the
most reasonable.
Obviously the heat eroding the lithosphere is provided by the plume underplating the
lithosphere. Just a small amount of heat carried by the plume is transferred to the
lithosphere since most of its heat stays in the mantle and is recycled by mantle convection
(i.e. in a subduction zone). Still the amount of heat transferred to lithosphere can
be expected to be proportional to the amount of material that is transported by the
plume and carries heat upwards. Thus also the volume of eroded lithosphere should
be proportional to the buoyancy flux. Unfortunately we are not looking at the volume
of eroded lithosphere, but at just the vertical dimension of this volume. If it would
be equally easy to erode the lithosphere in the horizontal directions as in the vertical
direction, the vertical erosion would be ∆ ∝ F 1/3 . But the lateral erosion must be
much easier because in the vertical direction the lithosphere becomes quickly stiffer for
shallower parts of the lithosphere. Thus for a scaling law for the buoyancy flux we can
15
expect ∆ ∝ F a where a is a positive value much lower than 1/3. Here I calculated
a = 0.0912, which satisfies quite well these expectations. This hypothesis, assuming that
the volume of the eroded plate is proportional to the buoyancy flux, could be tested only
by doing the same study for the volume of thinned lithosphere instead of the maximum
thinning.
Scaling laws set up in [van Hunen and Zhong, 2003] showed that the background viscosity
is related to the height of the elevated topography on a plume and the width of the swell
by a power-law. For the height of the swell this is H ∝ η −0.5 . Since the elevated
topography is partly caused by thinning of the lithosphere one could expect that the
scaling law for lithosphere thinning and the background viscosity is a similar power-law.
The scaling law found here is ∆ ∝ η −0.2 which satisfies these expectations quite well.
Comparing these two relationships suggests that there is probably no proportionality
between lithospheric thinning and the elevation of the swell, because H decreases faster
with higher viscosity than ∆. If there would be a proportional relationship between ∆
and H, both ∆ and H were likely to have a similar relation to η
For the Arrhenius activation energy E, van Hunen and Zhong [2003] found a exponential
function as scaling law that relates it to the topography height, which leads again to
the assumption that for lithosphere thinning an exponential function also best fits our
data. But an exponential function fits the data only with a rms-value of the residuals of
0.0397 which is one magnitude higher than for a power-law relation. Still the correlation
coefficient for the exponential function is -0.9717 which says that the fit is also acceptable.
A power-law would imply that the activation energy has not the same kind of influence
on thinning and elevation of the lithosphere. Doing the calculations for more than 6
values would probably be helpful for the decision whether a power-law or an exponential
relationship explains the process in a better way.
The scaling laws that express the relationships between lithospheric thinning and the parameters buoyancy flux, background viscosity and Arrhenius activation energy all turned
out to be best fitted by power-laws. Unfortunately not in all cases (especially for Arrhenius activation energy) these empirical laws can be very well explained physically.
As already mentioned, doing the same calculations with the total volume of eroded
lithosphere instead of the maximum thinning would probably provide data that can be
fitted better and give more reasonable scaling laws. Also for another definition of the
asthenosphere-lithosphere-boundary, e.g. by another isotherm or a certain value of viscosity, the results would change to some degree and give other information for a better
understanding of the physical processes responsible for the lithospheric swell.
16
References
U.R. Christensen. Convection with pressure and temperature dependent non-newtonian
rheology. Geoph. J. R. Astron. Soc., 77:242–284, 1984.
G.H. Davies. Ocean bathymetry and mantle convection, 1, large-scale flow and hotspots.
J. Geophys. Res., 83:10,467–10,480, 1988.
R.S. Detrick and S.T. Crough. Island subsidence, hot spots and lithospheric thinning. J.
Geophys. Res., 83:1236–1244, 1978.
X. Li, R. Kind, X. Yuan, I. Wölbern, and W. Hanka. Rejuvenation af the lithosphere by
the hawaiiian plume. Nature, 427:827–829, 2004.
W.B. Moore, G. Schubert, and P. Tackley. Three-dimensional simulations of plumelithosphere interaction at the hawaiian swell. Science, 279:1008–1011, 1998.
P. Olson. Hot spots, swells and mantle plumes. Magma Transport and Storage, edited
by M.P.Ryan:33–51, 1990.
N.M. Ribe. Through thick and thin. Nature, 427:793 – 795, 2004.
N.M. Ribe and U.R. Christensen. Three-dimensional modeling of plume-lithosphere interaction. J. Geophys. Res., 99(B1):669–682, 1994.
N.H. Sleep. Hot spots and mantle plumes: some phenomenology. J. Geophys. Res., 95:
6715–6736, 1990.
D. L. Turcotte and G. Schubert. Geodynamics. Cambridge Press, 2nd edition, April
2002.
J. van Hunen and S. Zhong. New insight in the hawaiian plume swell dynamics from
scaling law. Geophys. Res. Let., 30(15):1–4, 2003.
S. Zhong and A.B. Watts. Constraints on the dynamics of the mantle plumes from uplift
of the hawaiian islands. Earth Plan. Sci. Let., 203:105–116, 2002.
17

Download Report