Mississippi Embayment Project: Class
From GeoMod
Contents |
Abstract
To simulate the formation of the Mississippi Embayment, a model for thermal expansion (uplift) and a model for erosion were created using the Python code. These two models were then combined into one model used to simulate the uplift, erosion, and subsidence of the Mississippi Embayment.
Introduction
The Mississippi Embayment was formed when the faulted crust passed over the Bermuda hotspot in Mid-Cretaceous (Cox & VanArsdale, 2002). The Bermuda plume is thought to have caused the Mississippi Valley Graben complex to uplift 1 - 3 km (Cox & VanArsdale, 2002). This uplift was then eroded and as the crust cooled it subsided, creating a depression called the Mississippi Embayment.
Using the Pyhton code, a model for the Bermuda hotspot was created. This model shows the hotspot moving underneath the crust, and the uplift of the crust due to thermal expansion. The details of this model can be found at Mississippi Embayment Project: Smith and Bronson.
Using the Python code a model for the erosion of the crust was made. This model shows the erosion of the crust (a topographic high). The equations used and the detailed workings of the erosion model can be found at Mississippi Embayment Project: Bradshaw.
The two models were then combined and altered slightly to create a single model simulating the complete formation of the Mississippi Embayment.
- Some of the changes include:
- Changing units and timesteps so that are uniform throughout the model
- Changing the temperature model boundaries (east and west boundaries set equal to temperature in neighbor cell)
- Incorporating isostatic rebound in the model (producing subsidence)
Isostatic Rebound
Governing Equation: The mass of the crust equals the mass of the mantle displaced: Mc=Mmd
Discretized Equation:
After dz_i is calculated, it is added to Z_new to get the elevation of the crust after erosion and isostatic rebound have occurred.
Model Code
from visual import *
from raster_map import *
from visual.graph import *
def plot_line(data_array):
glist = []
ct = -1
for i in data_array:
ct+=1
glist.append((ct,i))
c = gcurve(pos=glist)
'''Note that the boundary conditions for erosion are by default constant head '''
''' Dimensions of grid'''
nrows = 41
ncols = 41
'''Model Parameters'''
K = 3.0 # thermal conductivty (W/mK)
dx = 100000.0 #meters
TT = 15.0 + 273.0#assumption: may need to take average of cell
NT = 8400.0 + 273.0# Hot spot temp. Kelvin
dz =32000.0 #meters
A = (dx * dx)
dT = NT-TT
dt = 500.0* 365.25 * 24.0 * 60.0 * 60.0
BT = 335.0 + 273.0 #Base Temperature of crust
q = K * (dT/dz)
C = 790 # W/mK specific heat capacity of granite
S = K*dt / (dx * dx * C)
AT = (TT + BT) / 2.0 #average temp cell: boundary
'''Erosion parameters'''
dt_E = dt
K_E = 50000.0/ (365.25 * 24.0 * 60.0 * 60.0)
S_E = K_E * dt_E / (dx*dx)
'''Plume = ball'''
vel = 0.0
t = 0
xvel = .3 / (365.25 * 24.0 * 60.0 * 60.0) ##meters per year
density_crust = 2600.0 * ones((nrows,ncols)) # density of granite
density_mantle = 3300 * ones((nrows,ncols)) #density of mantle
dz_crust = dz * ones((nrows,ncols)) #thickness of crust
'''Temperature array'''
T = AT * ones((nrows, ncols)) #intial conditons
Tmap_base = raster_map(zeros((nrows,ncols)), dx=dx)
Tmap = raster_map(T, dx=dx)
tcolor_scale = color_map(300.0, 1000.0, 1.0, colmin=vector(1,0,0), colmax=vector(0,1,0))
Tmap.line_3d(scale=1.0, color_map = tcolor_scale, center=1)
'''Elevation array'''
Z = raster_map(dz_crust, dx=dx)
scene.center = vector(dx*20.0,dx*20.0,4500.)
'''Create arrays and raster_map's for embayment with the top of
the model being the crust (cr_top), and the bottom of the model
being at 0.0 (cr_base).'''
cr_top = raster_map(dz_crust, dx=dx)
cr_base = raster_map(0.01 * ones((nrows,ncols)), dx=dx)
'''Set Z = cr_top'''
Z = cr_top
'''Create embayment as a layer_raster with the cr_top and cr_base'''
embayment = layer_raster(cr_top, cr_base)
'''Draw map of embayment'''
color_scale = color_map(31990.0, 32050.0, 5.0, colmin=vector(1,0,0), colmax=vector(0,1,0))
embayment.draw_layer(color_map=color_scale)
'''Molding a ball'''
ball = sphere(pos=(dx,-20*dx,dz), color=(0,0,1),radius=dx*.25)
'''make pos Q = ball.pos'''
##ball.pos = (dx,-20*dx,-dz-10000)
(r,c)=cr_top.get_rc(ball.pos)
'''Boundary conditions'''
Q = 0.0 * ones((nrows,ncols)) #hotspot
Q[r,c] = K * ((NT-AT)/dz) * A
nt = 0
runtime = 20000000.0 * 365.25 * 24.0 * 60.0 * 60.0
tt = 0.0
##for nt in range(1, nsteps+1):
while tt < runtime:
tt += dt
nt += 1
M = dx * dx * dz_crust[1:nrows-1, 1:ncols-1] * density[1:nrows-1, 1:ncols-1]
'''Temperature boundaries'''
T[:,0]=T[:,1] # Left column
T[:,ncols-1]=T[:,ncols-2] # Right column
T[0,:]=AT #top row
T[nrows- 1,:]=AT#bottom col
'''Temperature Diffusion Equation'''
temp_new = T[1:nrows-1, 1:ncols-1] + (Q[1:nrows-1, 1:ncols-1]*dt/(C*M))+ \
(S/(density[1:nrows-1, 1:ncols-1]*dz_crust[1:nrows-1, 1:ncols-1]))* \
((T[2:nrows, 1:ncols-1]- T[1:nrows-1, 1:ncols-1])*(dx*(dz_crust[1:nrows-1, 1:ncols-1]+dz_crust[2:nrows, 1:ncols-1])/2.0)- \
(T[1:nrows-1, 1:ncols-1] - T[0:nrows-2, 1:ncols-1])*(dx*(dz_crust[1:nrows-1, 1:ncols-1]+dz_crust[0:nrows-2, 1:ncols-1])/2.0)+ \
(T[1:nrows-1, 2:ncols]-T[1:nrows-1, 1:ncols-1])*(dx*(dz_crust[1:nrows-1, 1:ncols-1]+dz_crust[1:nrows-1, 2:ncols])/2.0)- \
(T[1:nrows-1, 1:ncols-1]-T[1:nrows-1, 0:ncols-2])*(dx*(dz_crust[1:nrows-1, 1:ncols-1]+dz_crust[1:nrows-1, 0:ncols-2])/2.0))
'''Thermal Expansion'''
Exp = (temp_new - T[1:nrows-1,1:ncols-1]) * 9. * 0.0000001 * dz_crust[1:nrows-1,1:ncols-1]
#print Exp.shape
#print dz_crust.shape
dz_crust[1:nrows-1,1:ncols-1] = dz_crust[1:nrows-1,1:ncols-1] + Exp
'''Erosion equation'''
Z_new = Z.data[1:nrows-1, 1:ncols-1] + S_E * (Z.data[2:nrows, 1:ncols-1] - (4 * Z.data[1:nrows-1, 1:ncols-1]) +
Z.data[0:nrows-2, 1:ncols-1] + Z.data[1:nrows-1, 2:ncols] +
Z.data[1:nrows-1, 0:ncols-2])
ze = Z.data[1:nrows-1,1:ncols-1] - Z_new
'''Isostatic rebound'''
dz_i = Z_new * (1.0 - density[1:nrows-1, 1:ncols-1]/density_mantle[1:nrows-1, 1:ncols-1]) - \
Z.data[1:nrows-1,1:ncols-1] * (1.0 - density[1:nrows-1, 1:ncols-1]/density_mantle[1:nrows-1, 1:ncols-1])
'''Update erosion'''
Z.data[1:nrows-1,1:ncols-1] = Z_new + dz_i
Z.data[1:nrows-1,1:ncols-1] = Z.data[1:nrows-1,1:ncols-1] + Exp
'''Update temperature'''
T[1:nrows-1,1:ncols-1] = temp_new
'''Move plume'''
ball.pos = ball.pos + vector(xvel*dt,0.0,0.0)
(r, c) = cr_top.get_rc(ball.pos)
'''Update plume heat flux'''
Q = 0.0 * Q
Q[r,c] = K * ((NT-T[r,c])/dz_crust[r,c]) * dx * dx *5000.0
'''Output'''
topo_cross = cr_top.data[21,:]-(32000. * ones((ncols)))
plot_line(topo_cross)
if nt % 100 == 0:
print Exp, ze
print "exp=", dz_crust[21,:]
print "position", ball.pos, r, c
print cr_top.data[21,:]
print "Q", NT, T[r,c], Q[r,c]*dt/(C*M[r,c])
print nt, T[21, :], ze[20, :]
Tmap.update_line_3d()
##plot_line(dz_crust[21,:])
embayment.update_layer()
'''Update image'''
Tmap.update_line_3d()
Model Assumptions
- Initial topography of the crust is flat (dz = 32000 m)
- Uniform erosion rate for the modeled area
- No effect from water runoff or streams on the erosion rate
- All erosional boundaries are no flow
- East and west temperature boundaries are set equal to the temperature of the cell next to it
- North and south temperature boundaries are no flow
- Density of the crust is equal to the density of granite
Model Simulation
The Python code written for the model was then displayed visually using VPython. The result of the simulation is shown in the following four figures, showing the elevation of the area as the model is running.
- Python model showing the elevation of the area, due to thermal expansion, erosion, and isostatic rebound, at the beginning of the model. Green is higher elevation, red is lower elevation, and the plume is the blue dot.
- Python model showing the elevation of the area, due to thermal expansion, erosion, and isostatic rebound, almost half of the way through the model. Green is higher elevation, red is lower elevation, and the plume is the blue dot. Note that on the left side of the figure, you can see subsidence starting to occur. Subsidence is marked by the lighter red color.
- Python model showing the elevation of the area, due to thermal expansion, erosion, and isostatic rebound, just over half of the way through the model. Green is higher elevation, red is lower elevation, and the plume is the blue dot. Note that on the left side of the figure, you can see the Mississippi Embayment forming. This is marked by the lighter red color.
- Python model showing the elevation of the area, due to thermal expansion, erosion, and isostatic rebound, at the end of the model. Green is higher elevation, red is lower elevation, and the plume is the blue dot. Note that the Mississippi Embayment has grown to about half of the modeled area. The Mississippi Embayment will continue to form as the crust cools and subsides.
Model Limitations
When the model is run it is clear that the heat spreads further laterally. This is due to the size of the model embayment.
The effect of uplift caused by the force of the hotspot is not included in this model.
The diffusion eguation used for temperature is only 2-D.
Future Work
To better improve the model the embayment should be enlarged and the force of uplift caused by the hotspot should be included. Since it is hard to measure the force of uplift caused by a hotspot this model could be used to predict the force by seeing home much up force would need to be added in addition to the rise in elevation created by thermal expansion alone.
References
Bishop, 2007: Bishop, P. 2007. Long-term landscape evolution: linking tectonics and surface processes. Earth Surf. Process. Landforms. Wiley Interscience: (www.interscience.wiley.com)
Cox & VanArsdale, 2002: Cox,R., VanArsdale,R. 2002. The Misssissippi Embayment, North America: a first order continental structure generated by the Cretaceous superplume mantle event. Journal of Geodynamics, Vol. 34. Pergamon. pp. 163-176
Densmore, Ellis, & Anderson, 1998: Densmore, A.L., Ellis, M.A., and Anderson, R.S. 1998. Landsliding and the evolution of normal faultbounded ranges, J. Geophys. Res., in press, 1998.
Guillou-Frottier et. al.: Plume head–lithosphere interactions near intra-continental plate boundaries
Martin & Church, 2004: Martin, Y. and Church, C. 2004. Numerical modeling of landscape evolution: Progress in Physical Geography 28,3. pp. 317–339
