Heat flux in FiPy with zero flux

zeroflux.py

##CREATING DATAFRAMES FOR SURFACE TEMPERATURE AND LITHOLOGY PROPERTIES 


#SURFACE TEMPERATURE 

import numpy as np 
import pandas as pd 

m = 10 #multiplyer; used so that if we want to run the model over a longer time only one value has to be changed 
time = 800 * m #in real time should be over 800.000 years 

surface_temp = pd.DataFrame(np.arange(time), columns=['time'])
surface_temp['temperature']=0
surface_temp['value']=0


#THIS LIST SHOULD PROBABLY BE CONVERTED AGAIN (FROM OLD T=0 TO CURRENT SURFACE TEMP)
surface_temp.temperature.iloc[0*m:10*m] = -5         #define surface temperatures
surface_temp.temperature.iloc[10*m:25*m] = 0
surface_temp.temperature.iloc[25*m:50*m] = -15
surface_temp.temperature.iloc[50*m:90*m] = 0
surface_temp.temperature.iloc[90*m:100*m] = -15
surface_temp.temperature.iloc[100*m:120*m] = 0
surface_temp.temperature.iloc[120*m:130*m] = -15
surface_temp.temperature.iloc[130*m:150*m] = 0
surface_temp.temperature.iloc[150*m:800*m] = -5

surface_temp.value[surface_temp.temperature == 0] = 1      #give the different surface temperatures a value
surface_temp.value[surface_temp.temperature == -5] = 2     #for later reading 
surface_temp.value[surface_temp.temperature== -15] = 3


#LITHOLOGY PROPERTIES  


df =pd.DataFrame(np.arange(6001), columns=['lithology'])
df['k'] = 0
df['dens'] = 0
df['cp'] = 0

toyr = 60 * 60 * 24 * 365 #how may seconds in a year
#something like this could be read from a csv file or petrel??
#DH4 example 
df.lithology.iloc[0:190] = 'sandstone'
df.lithology.iloc[190:540] = 'shale/siltstone'
df.lithology.iloc[540:590] = 'sandstone'
df.lithology.iloc[590:670] = 'shale/siltstone'
df.lithology.iloc[670:700] = 'sandstone'
df.lithology.iloc[700:770] = 'shale/siltstone'
df.lithology.iloc[770:800] = 'sandstone'
df.lithology.iloc[800:1200] = 'shale/siltstone'
df.lithology.iloc[1200:1230] = 'ign_intrusion'
df.lithology.iloc[1230:1500] = 'shale'
df.lithology.iloc[1500:6001] = 'sandstone' #actually not sandstone but equal properties



df.k[df.lithology=='sandstone'] = 3 * toyr #where df.lithology =='sandstone', a k value is assigned accordingly
df.k[df.lithology=='shale/siltstone'] = 2 * toyr 
df.k[df.lithology=='shale'] = 2 * toyr 
df.k[df.lithology=='ign_intrusion'] = 2 * toyr 

df.dens[df.lithology=='sandstone'] = 2200
df.dens[df.lithology=='shale/siltstone'] = 2400
df.dens[df.lithology=='shale'] = 2600
df.dens[df.lithology=='ign_intrusion'] = 2700

df.cp[df.lithology=='sandstone'] = 850
df.cp[df.lithology=='shale/siltstone'] = 850
df.cp[df.lithology=='shale'] = 850
df.cp[df.lithology=='ign_intrusion'] = 850

df['d'] = df.k / (df.dens * df.cp)

#HEATFLOW MODEL 

from fipy import *

nx = 6000        #amount of steps 
dx = 1.          #stepsize 
L = int(nx * dx) #depth of well 

mesh = Grid1D(nx=nx, dx=dx)

initial_temp = 180.


startingvalue_1 = (mesh.cellCenters[0]) / (nx/initial_temp)
startingvalue_2 = 0. #choose a starting value; defines with what kind of geothermal gradient you start

T = CellVariable(name="DH4", mesh=mesh, value= (startingvalue_1))


X = mesh.faceCenters[0]


#variables into FaceVariables
k= FaceVariable(mesh=mesh) #thermal conductivity
dens= FaceVariable(mesh=mesh) #density
cp = FaceVariable(mesh=mesh) #heat capacity



k.setValue(0)
dens.setValue(0)
cp.setValue(0)


for ii in range(len(df.d)): 
    if df.lithology.iloc[ii] == 'sandstone':
        k.setValue(3 * toyr, where=(X == ii))
        dens.setValue(2200, where=(X == ii))
        cp.setValue(850, where=(X == ii))
    elif df.lithology.iloc[ii] == 'shale/siltstone':
        k.setValue(2 * toyr, where=(X == ii))
        dens.setValue(2400, where=(X == ii))
        cp.setValue(850, where=(X == ii))
    elif df.lithology.iloc[ii] == 'ign_intrusion':
        k.setValue(2 * toyr, where=(X == ii))
        dens.setValue(2700, where=(X == ii))
        cp.setValue(850, where=(X == ii))
    elif df.lithology.iloc[ii] == 'shale':
        k.setValue(2 * toyr, where=(X == ii))
        dens.setValue(2600, where=(X == ii))
        cp.setValue(850, where=(X == ii))
#what happens here is that based on the lithologies, values are assigned to the lithologies 

D = k / (dens * cp)
#diffussion coeff is then calculated based on these properties per meter. 

time = Variable()
dt = 10 #time step

#T.constrain(300., mesh.facesRight)
#T.faceGrad.constrain([fluxRight], mesh.facesRight)

fluxRight = Variable([7.]) * dt #(J/yr) / m^2; according to CMR-report this should be 1893(J/yr)/m^2 or 70mW/m^2
fluxRight = fluxRight / (dens * cp) # m * K /yr

D = FaceVariable(mesh=mesh, value=D)
D.setValue(0., where=mesh.facesRight)

eqI = (TransientTerm() == DiffusionTerm(D) + (mesh.facesRight * fluxRight).divergence)
#diffusion equation + implemented heatflow 





phiAnalytical = CellVariable(name="analytical value", mesh=mesh, value=(mesh.cellCenters[0]) / (nx/initial_temp))
if __name__ == '__main__':
    viewer = Viewer(vars=(T, phiAnalytical), limits={'xmin': 0, 'xmax': 7000, 'ymin':-15, 'ymax':400}, title="temp vs depth")
    viewer.plot()



while time() < 800 * m:      #when choosing new m remember to re-run the first cell! 
    if surface_temp.value.iloc[int(time())] == 1:
        valueLeft = 0
        T.constrain(valueLeft, mesh.facesLeft)  
    elif surface_temp.value.iloc[int(time())] == 2:
        valueLeft = -5
        T.constrain(valueLeft, mesh.facesLeft) 
    elif surface_temp.value.iloc[int(time())] == 3: 
        valueLeft = -15
        T.constrain(valueLeft, mesh.facesLeft) 

    #here add the temperature of -6 which is the CMR approximation for the 
    #surface temp 150000 - 800000 yrs ago
  
    time.setValue(time() + dt)
    eqI.solve(var=T, dt=dt)
    if __name__ == '__main__':
        viewer.plot()