Mississippi Embayment Project: Smith and Bronson
From GeoMod
Contents |
Determine model domain
To determine the domain of the model, the parameters that to are to be set must be dependent upon the initial calculations of the average temperature of the crust (the "z" coordinate), the temperature gradient of the crust (the heat transfer), the heat supplied to the crust (the source) and making certain heat gain will be equal to loss.
Governing equations
Temperature Diffusion:
Thermal Expansion of Granite and Sandstones: 9*10^-6
Discretized equations
Thermal Expansion: Exp = (temp_new - T[1:nrows-1,1:ncols-1]) * 9. * 0.0000001 * dz_crust[1:nrows-1,1:ncols-1]
Temperature Diffusion:
Model parameters
In the parameters of the model, the results of running the model should demonstrate the possible effects of the moving matle over the plume within this particular grid. This approach is based upon the erosion portion of the uplift of the model and showing the removal of the uplift. The model is based upon the movement of the plume in millions of years or this could be reduced in division of tens.
Initial conditions
Utilizing the diffusion to show how energy is moved from one cell to the surrounding cells and so forth. The edge cells will be used to limit the diffusing of the input energy in order to stop the movement of the energy to give it an end. The movement of the earth's crust over the plume was done in an opposite manner where it appeared the plume was moving under the crust of the earth.
Programming
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 are by default constant head '''
''' set dimensions of grid'''
nrows = 41
ncols = 41
'''set 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 = 2200.0 + 273.0# Hot spot temp. Kelvin
dz =32000.0 #meters
A = (dx * dx)
dT = NT-TT
dt = 100.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
density = 2600.0 * ones((nrows,ncols)) # density of granite
dz_crust = dz * ones((nrows,ncols)) #thickness of crust
'''Temperature array'''
T = AT * ones((nrows, ncols)) #intial conditons
'''boundary temp'''
T[:,0]=AT # bottom row
T[:,nrows-1]=AT # top row
T[0,:]=AT #left col
T[ncols- 1,:]=AT #right col
print T
Tmap_base = raster_map(zeros((nrows,ncols)), dx=dx)
Tmap = raster_map(T, dx=dx)
color_scale = color_map(300.0, 400.0, 1.0, colmin=vector(1,0,0), colmax=vector(0,1,0))
Tmap.line_3d(scale=1.0, color_map = color_scale)
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), 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)
##
##'''Create embayment as a layer_raster with the cr 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))
##color_scale = color_map(50.0, 62.0, 1.0, colmin=vector(1,0,0), colmax=vector(0,1,0))
##'''draw lines for embayment top with the color_map'''
##crust.top.line_3d(color_map=color_scale)
##'''draw layer surfaces with the color_map (may not be stable'''
###aquifer.draw_layer_surfaces(color_scale)
##'''draw whole layer with colors determined by the elevation of the top layer'''
embayment.draw_layer(color_map=color_scale)
'''molding a ball'''
ball = sphere(pos=(0,-20*dx,-dz-10000), color=(0,0,1),radius=50000.0)
'''make pos Q = ball.pos'''
ball.pos = (dx,-20*dx,-dz-10000)
(r,c)=cr_top.get_rc(ball.pos)
'''moving ball'''
vel = 0.0
t = 0
xvel = .03 / (365.25 * 24.0 * 60.0 * 60.0) ##meters per year
##if ball.pos.x > - 1:
##
## vel += xvel * dt
## ball.pos.x += vel * dt
'''set boundary conditions'''
Q = 0.0 * ones((nrows,ncols)) #hotspot
Q[r,c] = K * ((NT-AT)/dz) * A
nsteps = 100000
for nt in range(1, nsteps+1):
M = dx * dx * dz_crust[1:nrows-1, 1:ncols-1] * density[1:nrows-1, 1:ncols-1]
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
#print Exp.shape
#print dz_crust.shape
dz_crust[1:nrows-1,1:ncols-1] = dz_crust[1:nrows-1,1:ncols-1] \
+ Exp * dz_crust[1:nrows-1,1:ncols-1]
'''update'''
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'''
if nt % 100 == 0:
print "exp=", dz_crust[21,:]
print "position", ball.pos, r, c
print "Q", NT, T[r,c], Q[r,c]*dt/(C*M[r,c])
print nt, T[21, :]
Tmap.update_line_3d()
##plot_line(dz_crust[21,:])
embayment.update_layer()
'''update image'''
Tmap.update_line_3d()
Results of model simulations
1. The first result was from the input of the values from the diffusion equation. This would show if the equation(s) were stable and if not, how, and what was needed to keep the program stable.
2. A grid was put into the program to relate to the expanse of the uplift and to give a ratio in distance between the appilachian and the rocky mountains.
3. Time steps had to be separated out of the running program in order to determine if the program was running with visual results and if not what could be done to demonstrate that without the program becoming unstable.
4. a ball program was inserted into the program in order to represent the plume, the movemnt of the plume, and the heat dissipated throught the heat from the plume.
