From 056994564bab91ea85cd84ff6eae97263b35dc11 Mon Sep 17 00:00:00 2001 From: Scott Date: Fri, 29 Sep 2017 13:36:28 -0600 Subject: [PATCH] Initial commit; minimal working code after Montek's summer work with Irina --- .gitignore | 159 ++++++++ oldercodewithViZ/Python_tha_test.py | 611 ++++++++++++++++++++++++++++ oldercodewithViZ/radin_dailyavg.csv | 1 + radin_dailyavg.csv | 1 + setup.py | 10 + thalake/bmithalake.py | 146 +++++++ thalake/radin_dailyavg.csv | 1 + thalake/thalake.py | 553 +++++++++++++++++++++++++ 8 files changed, 1482 insertions(+) create mode 100755 .gitignore create mode 100644 oldercodewithViZ/Python_tha_test.py create mode 100644 oldercodewithViZ/radin_dailyavg.csv create mode 100644 radin_dailyavg.csv create mode 100644 setup.py create mode 100644 thalake/bmithalake.py create mode 100644 thalake/radin_dailyavg.csv create mode 100644 thalake/thalake.py diff --git a/.gitignore b/.gitignore new file mode 100755 index 0000000..a6a8f42 --- /dev/null +++ b/.gitignore @@ -0,0 +1,159 @@ +# ---------------------------------------------------------------------------- # +### Added from: +### https://github.com/github/gitignore/blob/master/Python.gitignore +### on 10/7/2016 + +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +env/ +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +*.egg-info/ +.installed.cfg +*.egg + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*,cover +.hypothesis/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +target/ + +# IPython Notebook +.ipynb_checkpoints + +# pyenv +.python-version + +# celery beat schedule file +celerybeat-schedule + +# dotenv +.env + +# virtualenv +.venv/ +venv/ +ENV/ + +# Spyder project settings +.spyderproject + +# Rope project settings +.ropeproject + +# ---------------------------------------------------------------------------- # +### Added from: +### https://github.com/github/gitignore/blob/master/Global/Vim.gitignore +### on 10/7/2016 +# swap +[._]*.s[a-w][a-z] +[._]s[a-w][a-z] +# session +Session.vim +# temporary +.netrwhist +*~ +# auto-generated tag files +tags + +# ---------------------------------------------------------------------------- # +### Added from: +### https://github.com/github/gitignore/blob/master/Global/Emacs.gitignore +### on 10/7/2016 +# ---------------------------------------------------------------------------- # + +# -*- mode: gitignore; -*- +*~ +\#*\# +/.emacs.desktop +/.emacs.desktop.lock +*.elc +auto-save-list +tramp +.\#* + +# Org-mode +.org-id-locations +*_archive + +# flymake-mode +*_flymake.* + +# eshell files +/eshell/history +/eshell/lastdir + +# elpa packages +/elpa/ + +# reftex files +*.rel + +# AUCTeX auto folder +/auto/ + +# cask packages +.cask/ +dist/ + +# Flycheck +flycheck_*.el + +# server auth directory +/server/ + +# projectiles files +.projectile diff --git a/oldercodewithViZ/Python_tha_test.py b/oldercodewithViZ/Python_tha_test.py new file mode 100644 index 0000000..51f124b --- /dev/null +++ b/oldercodewithViZ/Python_tha_test.py @@ -0,0 +1,611 @@ +# Thaw Lake Model-1D + +# This model a 1-D numerical model of permafrost and subsidence processes. +# It aims to investigate the subsurface thermal impact of thaw lakes of various depths, +# and to evaluate how this impact might change in a warming climate. + +# Key paper: Matell, N., Anderson, R.S., Overeem, I., Wobus, C., Urban, F., +# Clow, G., in review 2011. Modeling the subsurface thermal impact of Arctic thaw lakes in a warming climate. Computers and Geosciences. + +# Copyright (C) <2011> + +# Developer can be contacted by irina.overeem@colorado.edu + +# Dr. Irina Overeem +# CSDMS Community Surface Dynamics Modeling System +# INSTAAR, University of Colorado at Boulder +# PO Box 450, 80309-0450 +# Boulder, CO, USA + + +# This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. +# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +# You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +# This model a 1-D numerical model of permafrost and subsidence processes. +# It aims to investigate the subsurface thermal impact of thaw lakes of various depths, +# and to evaluate how this impact might change in a warming climate. + +# The model was designed for the Alaskan Arctic Coastal Plain +# The model uses average and amplitude of temperature from observed permafrost temperatures at 5cm in the subsurface at Drew Point data (USGS) +# Clow, G.D., 2008a. Continued permafrost warming in northern Alaska, 2008 update. NOAA/ESRL Global Monitoring Annual Conference, Boulder. +# LINE 32 CLOW___________ 24 p. error + +# The model uses observed heat flow for thermal gradient from: +# Lachenbruch, A.H., Sass, J.H., Marshall, B.V., Moses, T.H., Jr., 1982. Temperatures, heat flow, and the geothermal regime at Prudhoe By, Alaska. Journal of Geophysical Research 87, 9301-9316. + + +# in 1-D, subsurface temperature changes through time and space of a semi-infinite +# half-space, then plots result. The model code incorporates phase change by using +# the apparent heat capacity scheme for temperatures within the +# "phase-change envelope". + +# lake-permafrost model - lake freezes and thaws, permafrost changes temperature, +# when permafrost thaw it subsides, assuming that if excess ice melts all +# the water leaves the permafrost and thus the subsurface volume decreases. + +# boundary condition of ice bottom and lake water top = Tf +# if TsTf, top +# boundary condition --> Tf and excess energy is used to melt ice. When +# ice melted from the top, all other boundary just moved up by equivelant +# amount, if moved up a full control volume then bottom control volume +# added that is same temp as (end-1) control volume. Assumes that when ice +# free, lake is completely mixed. (ie Liston & Hall) + +import math +import numpy as np +import csv +import operator +import copy +import scipy.io as sio +import matplotlib.pyplot as plt + +# set up 1-d grid +dz0 = 0.05; # cell size in top 10 m +dz2 = 1; # cell size below 10 m +deptht = 100; # model depth [m] + +z_range1 = np.arange(0,10.05,dz0); +z_range2 = np.arange(11,deptht+1,dz2); +z = np.append(z_range1, z_range2); # vertical grid (distance below surface) [m] + +dz = np.diff(z); +zbetween = z[1:] + np.diff(z)/2; +dzbetween = np.diff(zbetween); +numcells = len(z) # number of grid cells +iceinit = 0.0001; # initial ice depth [m] +celltype = np.int16(1*np.ones(numcells)); # substrate type - 1=regular soil with excess ice, 2=compacted soil without excess ice, 8=water, 9=ice +celltype[np.where(z<6.7)] =2; # depth of subsided cells, this is a function of the initial lake depth (seed lake depth/excess ice content) +celltype[np.where(z<2.0)] = 8; # seed lake depth +celltype[np.where(z<=iceinit)] = 9; # ice thickness +celltype[np.where(z>=50)] = 2; # underlying cells with no excess ice +thawed = np.zeros(numcells); +thawedspec = np.zeros(numcells); + +deptht = max(z); +depthsubside = 2.0; # seed lake depth, has to match with ln 79, because the lake contains water +depthtalik = 0; # usually unkown, so the model will trend toward equilibrium in the first 100-200 years of a simulation +depthi = iceinit; +print (depthi) + # set time steps and simulation duration +periodyear = 3600*24*365.; # period (1 year) [s] +dt = 3600*24; # length of timestep [s] +years = 1000; # simulation duration +nt = int(((periodyear/dt)*years)+1); # number of timesteps +t = np.arange(0,dt*nt+1,dt); # time at each timestep [s] +tday = [x/(3600*24) for x in t]; +tyear = [x/(3600*24*365) for x in t]; +tyearsonly = np.arange(0,years); +yearnum = 0; + +# declare zarrays + +depthin = np.zeros(nt); +depthwn = np.zeros(nt); +depthsubsiden = np.zeros(nt+1); +depthtalikn = np.zeros(nt+1); +icethickness = np.zeros(nt+1); +icegrowth = np.zeros(nt+1); +soilthickness = np.zeros(nt+1); +waterthickness = np.zeros(nt+1); +Tbot = np.zeros(nt+1); +T3m = np.zeros(nt+1); +T5m = np.zeros(nt+1); +T10m = np.zeros(nt+1); +T25m = np.zeros(nt+1); +T50m = np.zeros(nt+1); +T100m = np.zeros(nt+1); +Tbotavgann = np.zeros(years); +T3mavgann = np.zeros(years); +T5mavgann = np.zeros(years); +T10mavgann = np.zeros(years); +T25mavgann = np.zeros(years); +T50mavgann = np.zeros(years); +T100mavgann = np.zeros(years); +Trecord = np.zeros((365,numcells)); +watts = np.zeros(np.size(t)); + +# define constants +# Model input and boundary conditions (compiled from field and modeling studies by West & Plug, 2008; Ling and Zhang, 2003;..... +# Ling & Zhang, 2004; Zhou and Huang, 2004; Romanovsky and Osterkamp, 2000; Hinzman, 1998; and French, 2007). + +kelvin = 273.15; + +rhoi = 917; # ice density [kg/m3] +ci = 2108; # ice specific heat capacity [J/kg/K] +ki = 2.18; # ice thermal conductivity [J/s/m/K] +kappai = ki/(ci*rhoi); +rhow = 1000; # water density [kg/m3] +cw = 4210; # water specific heat capacity [J/kg/K] +kw = 0.58; # water thermal conductivity [J/s/m/K] +kappaw = kw/(cw*rhow); +rhorock = 1200; # mineral soil density [kg/m3] +crock = 1000; # mineral soil specific heat capacity [J/kg/degree C] +krock = 1.5; # rock thermal conductivity [J/s/m/degree C] + +totalice = 0.65; # ice content of original substrate, from bulk samples in top 5m at Drew Point +excessice = 0.3; # excess ice content - gone once melted +porewater = totalice-excessice; # porewater ice - still there if refreezes +W = totalice; # total water content percent of material by mass +Wu = totalice*0.05; # unfrozen percent water content of material by mass at + #temperature T = Tpc-pcenv +Ws = porewater; +Wus = porewater*0.05; +totalice2 = totalice-Wu; +porewater2 = porewater-Wus; + +cf = (ci*totalice2)+(cw*Wu)+(crock*(1-totalice)); +cu = (cw*totalice)+(crock*(1-totalice)); +cfs = (ci*porewater2)+(cw*Wu)+(crock*(1-porewater)); +cus = (cw*porewater)+(crock*(1-porewater)); +rhof = (rhoi*totalice2)+(rhow*Wu)+(rhorock*(1-totalice)); +rhou = (rhow*totalice)+(rhorock*(1-totalice)); +rhofs = (rhoi*porewater2)+(rhow*Wus)+(rhorock*(1-porewater)); +rhous = (rhow*porewater)+(rhorock*(1-porewater)); +Cf = cf*rhof; # frozen volumetric heat capacity [J/m3/degree C] +Cu = cu*rhou; # thawed volumetric heat capacity [J/m3/degree C] +Cfs = cfs*rhofs; +Cus = cus*rhous; +kf = pow(ki,totalice2)*pow(kw,Wu)*pow(krock,(1-totalice)); +ku = pow(kw,totalice)*pow(krock,(1-totalice)); +kfs = pow(ki,porewater2)*pow(kw,Wus)*pow(krock,(1-porewater)); +kus = pow(kw,porewater)*pow(krock,(1-porewater)); + +L = 334*1000; # latent heat of fusion for water [J/kg] +hw = 0.56; # convective transfer coefficient =[J/s/m2/K] +Kh = 0.0; # turbulent diffusivity for heat +Tpc = kelvin+0; # freezing temperature [K] +pcenv = 1; # width of phase change envelope [m] + +# set up surface and initial temperatures, including geothermal gradient +# theoretical warming scenario is set here +Tbar = kelvin-11; # MAAT [K] +warming = 5; +Tbarplus = Tbar+(float(warming)/36500)*np.arange(1,36501); +Tamp = 17.5; # amplitude of air temperature fluctuations [K] + +Tsurface = [Tbar-Tamp*math.sin(2*math.pi*x/periodyear) for x in t]; +Tsurface[182499:218999] = Tbarplus[0:]-Tamp*np.sin(2.*np.pi)*t[182499:218999]/periodyear; +Tsurface[218999:] = Tbar+warming-Tamp*np.sin(2*np.pi)*t[218999:]/periodyear; + +q = 0.056; # mantle heat flow [J/m2/s] +dTdzbase = q/kf; +dTdzbase2 = q/kfs; +Tgrad_first_half = [np.multiply(dTdzbase,z[np.where(z<=10)])] +Tgrad_second_half= [np.multiply(dTdzbase2,z[np.where(z>10)])]; +Tgrad = np.append(Tgrad_first_half, Tgrad_second_half); + +Twi = Tsurface[0]; # start up water temperature [K] +Tsi = Tbar; # start up permafrost temperature [K] +Tinit = np.ones((np.size(z))); # initial temperature grid +Tinitlake_first_half = [Tpc*Tinit[np.where(celltype == 9)]] +Tinitlake_second_half= [Twi*Tinit[np.where(celltype == 8)]]; +Tinitlake = np.append(Tinitlake_first_half, Tinitlake_second_half); +water = len(Tinitlake); +for n in range(0,numcells): + if (celltype is 9): + Tinit[n] = Tpc; + elif (celltype is 8): + Tinit[n] = Twi; + else: + Tinit[n] = Tsi+Tgrad[n]; + +# load daily average radiation, function from Drew Point meteorological data +mat_file = np.genfromtxt('radin_dailyavg.csv', delimiter=','); + +#Tinitperm = [Tsi*Tinit(find(celltype==1))+Tgrad(find(celltype==1)),Tsi*Tinit(find(celltype==2))+Tgrad(find(celltype==2))]; +#Tinit = [Tinitlake, Tinitperm]; + +T = Tinit; +zstarthaw = math.sqrt((ku/Cu)*periodyear/(math.pi)); +zstarfrozen = math.sqrt((kfs/Cf)*periodyear/(math.pi)); +Tbase = Tbar+(dTdzbase2*deptht)+np.multiply(Tamp,np.exp(-deptht/zstarthaw))*np.sin((2*math.pi*np.divide(t,periodyear))-(deptht/zstarthaw)); +Tright = -kelvin+Tbar+np.multiply(dTdzbase2,z)+np.multiply(Tamp,np.exp(-z/zstarthaw)); #outer edges of the funnel +Trightf = -kelvin+Tbar+np.multiply(dTdzbase2,z)+np.multiply(Tamp,np.exp(-z/zstarfrozen)); #outer edges of the funnel +Tleftf = -kelvin+Tbar+np.multiply(dTdzbase2,z)-np.multiply(Tamp,np.exp(-z/zstarfrozen)); +Tleft = -kelvin+Tbar+np.multiply(dTdzbase2,z)-np.multiply(Tamp,np.exp(-z/zstarthaw)); +#Ts0 = Tsurface[1]; + +count = 1; +daynum = 0; +numplots = 50; +countprint = dt*nt/numplots; +nplot=0; +fig_one = plt.figure(1) +zero = np.zeros((np.size(z))); +plt.plot(zero,-z, linewidth = 2); +plt.plot(Tright,-z,Tleft,-z,Trightf,-z,Tleftf,-z,) +plt.plot(Tinit-kelvin,-z) +plt.axis([Tbar-Tamp-kelvin, Tbar+Tamp-kelvin, -deptht, 0]); +plt.xlabel('Temperature (C)',fontsize = 18) +plt.ylabel('Depth (m)',fontsize = 18) +#fig_one = plt.show(); + + +# model equations +print ('starting now') +loopcount = 0; +print (nt); +for n in range(0,nt): + loopcount = loopcount+1; + print ('____LOOPCOUNT ' , loopcount); + thawedspec = np.zeros(numcells); + thawed = np.zeros(numcells); + count = count+1; + time = (n+1)*dt; # time into run [s] + error = 1; + Told = T; # remember temperatures from last timestep... + Ts0 = Tsurface[n]; + T[0]= Ts0; + + if depthi>0 and Tsurface[n]>Tpc: + Ts0 = Tpc; + + water = -1; + ice = -1 + subsided = -1; + regsoil = -1; + + for i in range(0,numcells): + if celltype[i] == 8: + water = water+1; + elif (celltype[i] == 9): + ice = ice+1; + #print 'updated ice is ', ice; + elif (celltype[i] == 2): + subsided = subsided+1; + elif (celltype[i] == 1): + regsoil = regsoil+1; + + print ('ice is =', ice, ' ', 'Water is = ', water); + + day = (n % 365); + solarrad = mat_file[day]; + #print "water on top _ ", water; + #compute temperatures in water layer if no ice exists + #water = 0; + #ice = 38; + if celltype[0] == 8: + print ('in if loop_________') + Twater = T[0:water]; # T(find(celltype==8)); + Tmix = np.mean(Twater); #mix again + sextinc = 0.6; # solar extinction constant + albedo = 0.06; + qrad = (1-albedo)*solarrad*np.exp(-sextinc*zbetween[0:water]); + watts[count] = solarrad; + error = 1; + count = 0 + while error>0.0001: + Tcalc = [0] + Tit = Twater; + for i in range (1,water-1): + coeff1 = kw/dz[i]; + coeff2 = kw/dz[i-1]; + coeff3 = cw*rhow*dzbetween[i-1]/dt; + coeff4 = coeff1+coeff2+coeff3; + Tcalc = np.append(Tcalc, (coeff1*T[i+1]+coeff2*T[i-1]+coeff3*Told[i]+qrad[i-1]-qrad[i])/coeff4); + Twater_Partial = np.append(Ts0,Tcalc[1:]); + Twater = np.append(Twater_Partial, Tcalc[-1]); + error = max(abs(Twater-Tit)); + Twater2 = Twater[0:-1]; + Tmix = np.mean(Twater2); # thoroughly mix lake + T[1:water] = [Tmix*np.ones(np.size(x))for x in Twater2]; + # compute temperatures in ice layer if completely frozen to bottom + + elif water == 0: + print ('___________________in elif loop___________________') + Tice = (T[np.where(celltype==9)]); + #print 'Starting in the loop ', Tice + while error>0.0001: + Tcalc = [0] + #Tcalc= np.zeros(ice); + Tit = copy.deepcopy(Tice); + for i in range(1,ice): + coeff1 = ki/dz[i]; + coeff2 = ki/dz[i-1]; + coeff3 = ci*rhoi*dzbetween[i-1]/dt; + coeff4 = coeff1+coeff2+coeff3; + Tcalc = np.append(Tcalc, (coeff1*T[i+1]+coeff2*T[i-1]+coeff3*Told[i])/coeff4); + #print ' Tcalc is ', Tcalc; + Tice = np.append(Ts0,Tcalc[0:]); + #print 'Ts0 is ', Ts0; + #print "The final___ TICE ", Tice; + #print 'Tice length is ', len(Tice); #(Tice is 38) + #print 'Tit length is ', len(Tit); #(Tit is 39) + error = max(abs(Tice-Tit)); + #print error; + + T[np.where(celltype==9)] = Tice; + #compute temperatures in ice and water layers if both exist + + else: + print ('__in else loop') + print (ice); + if Tsurface[n]>=Tpc: + T[0:ice] = Tpc; + elif ice<3: + T[0] = Tsurface[n]; + T[ice] = Tpc; + print ("whatss elif"); + else: + Ti = T[0:ice]; + error = 1; + while error>0.0001: + Tcalci = [0] + #Tcalci= np.zeros(ice); + Titi = copy.deepcopy(Ti); + for i in range(1,ice-1): + coeff1 = ki/dz[i]; + coeff2 = ki/dz[i-1]; + coeff3 = ci*rhoi*dzbetween[i-1]/dt; + coeff4 = coeff1+coeff2+coeff3; + Tcalci = np.append(Tcalci, (coeff1*T[i+1]+coeff2*T[i-1]+coeff3*Told[i])/coeff4); + Ti_partial=np.append(Ts0,Tcalci[1:]); + Ti = np.append(Ti_partial,Tpc); + error = max(abs(Ti-Titi)); + T[0:ice] = Ti; + + error = 1; + if water==1: + T[np.where(celltype==8)] = Tpc; + else: + print ("this loop maybe"); + error = 1; + while error>0.0001: + #Tcalcw = [0] + Tit = copy.deepcopy(T); + Tcalcw= np.zeros(ice); + for i in range (ice,ice+water+1): + coeff1 = kw/dz[i]; + coeff2 = kw/dz[i-1]; + coeff3 = cw*rhow*dzbetween[i-1]/dt; + coeff4 = coeff1+coeff2+coeff3; + Tcalcw = np.append(Tcalcw, (coeff1*T[i+1]+coeff2*T[i-1]+coeff3*Told[i])/coeff4); + Tw = np.append(Tpc, Tcalcw[ice:]); #problem appears to be in Tcalcw[ice:] + T[ice:ice+water+2] = Tw; + error = max(abs(T-Tit)); + + # Calculate temperatures in underlying permafrost + #clear 'z2' 'zbetween2' 'Told2' 'T2' 'celltype2' + z2 = z[ice+water-1:]; + zbetween2 = zbetween[ice+water-1:]; + Told2 = Told[ice+water-1:]; + T2 = T[ice+water-1:]; + numcells2 = len(T2); + celltype2 = celltype[ice+water-1:]; + C = np.zeros(numcells2); + k = np.zeros(numcells2); + Tcalc = 0; + for i in range (0,numcells2): + if celltype2[i] == 8: + C[i] = cw*rhow; + k[i] = kw; + elif celltype2[i] == 9: + C[i] = ci*rhoi; + k[i] = ki; + elif celltype2[i]==1 and T2[i] < (Tpc-pcenv): + C[i] = Cf; + k[i] = kf; + #print('real permafrost, frozen') + elif celltype2[i]==1 and T2[i]>=(Tpc-pcenv) and T2[i]<=Tpc: + C[i] = Cf+L*rhof*((W-Wu)/pcenv); + k[i] = kf+((ku-kf)/pcenv)*(T2(i)-(Tpc-pcenv)); + print('real permafrost, thawing') + elif celltype2[i]==1 and T2[i]>Tpc: + C[i] = Cu; + k[i] = ku; + print('real permafrost, thawed') + elif celltype2[i]==2 and T2[i]<(Tpc-pcenv): + C[i] = Cfs; + k[i] = kfs; + elif celltype2[i]==2 and T2[i]>=(Tpc-pcenv) and T2[i]<=Tpc: + C[i] = Cfs+L*rhofs*((Ws-Wus)/pcenv); + k[i] = kfs+((kus-kfs)/pcenv)*(T2[i]-(Tpc-pcenv)); + elif celltype2[i]==2 and T2[i]>Tpc: + C[i] = Cus; + k[i] = kus; + + kbetween = k[1:]-(np.diff(k)/2); + dz2 = np.diff(z2); + dzbetween2 = np.diff(zbetween2); + error = 1; + #Tcalcperm= np.zeros(251) + while error>0.0001: + Tcalcperm = [0] + Tit = copy.deepcopy(T2); + for i in range(1,numcells2-1): + coeff1 = kbetween[i]/dz2[i]; + coeff2 = kbetween[i-1]/dz2[i-1]; + coeff3 = C[i-1]*dzbetween2[i-1]/dt; + coeff4 = coeff1+coeff2+coeff3; + Tcalcperm = np.append(Tcalcperm, ((coeff1*T2[i+1]+coeff2*T2[i-1]+coeff3*Told2[i])/coeff4)); + #Tcalcperm[i] = (coeff1*T2[i+1]+coeff2*T2[i-1]+coeff3*Told2[i])/coeff4; + T2_partial = np.append(T[ice+water],Tcalcperm[1:]) + T2 = np.append(Tcalcperm,Tbase[n]); + error = max(abs(T2-Tit)); + T[ice+water-1:] = T2; + + # subside permafrost + print (len(range(water+ice,numcells))); + for i in range (water+ice,numcells): + if T[i]Tpc and celltype[i] == 1: + thawedspec[i] = 1; + thawed[i] = 1; + celltype[i] = 2; + #print 'is it here' + elif T[i]>Tpc and celltype[i] == 2: + thawed[i] = 1; + #print 'where is it' + +##[Tbar-Tamp*math.sin(2*math.pi*x/periodyear) for x in t]; + numthawedspec = np.sum(thawedspec)+1; + depthsubsidenew = z[numthawedspec]*excessice; + depthsubside = depthsubside+depthsubsidenew; + depthsubsiden[n] = depthsubside; + cellsubside = np.int16(math.floor((depthsubside)/dz0)); + maxcellsubside = cellsubside; + #if depthsubside>0: #depthsubsidenew>0 + # print celltype; + # for i in range(ice,maxcellsubside): + # celltype[i] = 8; + + if (Tsurface[n]ice+water: + icecheckdepth = ice+water-1; + Twater = T[icecheckdepth]; + if water==0: + Twater = T[np.where(celltype==8)]; + Twater = Twater[0]; + print ("here"); + if Twater 0: + ddepthidt_partial = ((Tpc-Ts0)*pow((depthi/ki),(-1))-hw*(Twater-Tpc)); + ddepthidt = np.divide(ddepthidt_partial, rhoi*L) + ddepthi = ddepthidt*dt; + #### ADJUSTMENT HERE + if ddepthi>0.05: + ddepthi = 1.05; + print('__________________had to adjust________________'); + icegrowth[n] = ddepthi; + depthi = depthi+ddepthi; + + if depthi<0.00001: + depthi = 0; + + depthw = depthsubside-depthi; + if depthi>depthsubside: + depthi = depthsubside; + depthw = 0; + + # check if Ts>Tpc (if there is ice there already) + elif Tsurface[n]>Tpc and depthi>0: + Qm = (Tsurface[n]-Tpc)*(pow((dz0/ki),(-1))); + dMdt = Qm/(L*rhoi); + dM = dMdt*dt; + icegrowth[n] = -dM; + depthi = depthi-dM; + if depthi<0: + depthi = 0; + + maxcelli = np.int16(round(depthi/dz0)); + celltype[0:ice+water] = 8; + + celltype[0:maxcelli] = 9; + icethickness[n] = depthi; + soilthickness[n] = deptht-depthsubside; + waterthickness[n] = deptht-icethickness[n]-soilthickness[n]; + Tbot[n] = T[ice+water]-kelvin; + T3m[n] = T[59]-kelvin; + T5m[n] = T[99]-kelvin; + T10m[n] = T[199]-kelvin; + T25m[n] = T[214]-kelvin; + T50m[n] = T[239]-kelvin; + T100m[n] = T[289]-kelvin; + if (tday[n]%365)==0 and tday[n]>364: + yearnum = yearnum+1; + Tbotavgann[yearnum] = np.mean(Tbot[n-364:n]); + T3mavgann[yearnum] = np.mean(T3m[n-364:n]); + T5mavgann[yearnum] = np.mean(T5m[n-364:n]); + T10mavgann[yearnum] = np.mean(T10m[n-364:n]); + T25mavgann[yearnum] = np.mean(T25m[n-364:n]); + T50mavgann[yearnum] = np.mean(T50m[n-364:n]); + T100mavgann[yearnum] = np.mean(T100m[n-364:n]); + + if n>=nt-364: + daynum = daynum+1; + Trecord[daynum, 0:] = T; + + if((n % (365*50))==0): + # nplot = nplot+1 + fig_one = plt.figure(1) + # ice + icethickness[n] + # water + waterthickness[n] + + lakedepthx = np.arange(-20,25,1); + lakedepth = icethickness[n]+waterthickness[n]*np.ones(np.size(lakedepthx)); + Tplot = T-kelvin; + plt.plot(Tplot,-z,'m') + timeplot = tday[n]/365 #will stamp time in years on your screen + plt.xlabel('temperature',fontsize = 18) + plt.ylabel('depth', fontsize = 18) + # + #figure(2) + #clf(2) + icedepth = icethickness[n]*np.ones(np.size(lakedepthx)); + plt.plot(zero,-z, linewidth = 2); + plt.plot(lakedepthx,-lakedepth, lakedepthx,-icedepth,Tplot,-z) + plt.xlabel('temperature', fontsize = 18) + plt.ylabel('depth', fontsize = 18) + #axis([-30 25 -5 0]) + +# Plot output parameters after run is completed +print ('after the loop') +'''# plot surface temperature fluctuations +fig_two = plt.figure(2) +TsK = 100;#Tsurface-kelvin; +plt.plot(tday,TsK) +plt.xlabel('time(days)',fontsize = 18) +plt.ylabel('surface temperature(degrees C)',fontsize = 18) + +# plot thickness +fig_three = figure(3) +plt.plot(tyear[1:end-1],icethickness[1:end-1],tyear[1:end-1],waterthickness[1:end-1],tyear[1:end-1],depthsubsiden[1:end-1]) +plt.xlabel('time(years)','fontname','arial','fontsize',18) +plt.ylabel('thickness (m)','fontname','arial','fontsize',18) +plt.legend('ice thickness','water thickness','depth subside') + +# plot temps at various depths +fig_four = figure(4) +T0C = np.zeros(size(tyear)); +plt.plot(tyearsonly,Tbotavgann,tyearsonly,T3mavgann,tyearsonly,T5mavgann,tyearsonly,T10mavgann,tyearsonly,T25mavgann, tyearsonly,T50mavgann,tyearsonly,T100mavgann,tyear[1:end-1],T0C[1:end-1]) +plt.xlabel('time(years)','fontname','arial','fontsize',18) +plt.ylabel('temperature (C)','fontname','arial','fontsize',18) +plt.legend('at lake bottom','at 3m depth','at 5m depth','at 10m depth','at 25m depth','at 50m depth','at 100m depth','0C') + +# plot temps at lake bottom +fig_five = figure(5) +T0C = np.zeros(size(tyear)); +plt.plot(tyear[1:end-1],Tbot[1:end-1],tyearsonly,Tbotavgann,tyear[1:end-1],T0C[1:end-1]) +plt.xlabel('time(years)','fontname','arial','fontsize',18) +plt.ylabel('temperature (C)','fontname','arial','fontsize',18) +plt.legend('at lake bottom','lake bottom average annual temp','0C') + +fig_two = plt.show(); +fig_three = plt.show(); +fig_four = plt.show(); +fig_five = plt.show(); + +input() +''' diff --git a/oldercodewithViZ/radin_dailyavg.csv b/oldercodewithViZ/radin_dailyavg.csv new file mode 100644 index 0000000..613391f --- /dev/null +++ b/oldercodewithViZ/radin_dailyavg.csv @@ -0,0 +1 @@ +6.286257077690969908e+01,6.094700543309402718e+01,5.902912919992142093e+01,5.710961000726149450e+01,5.518911863055893008e+01,5.326832822565901182e+01,5.134791386111293576e+01,4.942855204884011044e+01,4.751092027402129503e+01,4.559569652509747328e+01,4.368355882474016738e+01,4.177518476265366587e+01,3.990277293672754411e+01,3.831239339362680596e+01,3.676791799077611955e+01,3.522837210237122463e+01,3.369430585545826062e+01,3.216626720050992816e+01,3.064480155131714412e+01,2.913045143102182166e+01,2.762375612485181975e+01,2.612525134012183514e+01,2.463546887401730601e+01,2.315493628966762074e+01,2.171288715035241879e+01,2.055301198811046959e+01,1.942724853294134846e+01,1.831007557459968993e+01,1.720188784544575711e+01,1.610307562389593627e+01,1.501402452872270921e+01,1.393511532187155133e+01,1.286672371998985120e+01,1.180922021482493989e+01,1.076296990263483444e+01,9.728332322717138325e+00,8.937502651117220864e+00,8.216318407411662861e+00,7.504082120964423730e+00,6.801036929238594375e+00,6.107421582998850518e+00,5.423470365354699396e+00,4.749413017312665630e+00,4.085474669784134427e+00,3.431875781970633810e+00,2.805127455544981441e+00,2.420881179340481193e+00,2.048159507583730221e+00,1.682114724363046321e+00,1.322866093415396849e+00,9.705299341290687387e-01,6.252196071109451525e-01,2.918960039338246437e-01,1.771454084358596581e-01,6.928392916400623325e-02,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,6.928392916400623325e-02,1.771454084358596581e-01,2.918960039340988133e-01,6.252196071109380471e-01,9.705299341290672954e-01,1.322866093415394184e+00,1.682114724362788749e+00,2.048159507583477090e+00,2.420881179340474088e+00,2.805127455545747051e+00,3.431875781970863404e+00,4.085474669784111335e+00,4.749413017312162033e+00,5.423470365354933875e+00,6.107421582999201348e+00,6.801036929238826190e+00,7.504082120964650215e+00,8.216318407411169034e+00,8.937502651117322117e+00,9.728332322717884395e+00,1.076296990263472075e+01,1.180922021482483153e+01,1.286672371998961495e+01,1.393511532187169344e+01,1.501402452872334514e+01,1.610307562389588654e+01,1.720188784544571448e+01,1.831007557459963309e+01,1.942724853294126675e+01,2.055301198811109487e+01,2.171288715035418804e+01,2.315493628966740047e+01,2.463546887401664165e+01,2.612525134012043893e+01,2.762375612485159948e+01,2.913045143102206325e+01,3.064480155131664674e+01,3.216626720050921762e+01,3.369430585545777035e+01,3.522837210236924221e+01,3.676791799077442846e+01,3.831239339362482355e+01,3.990277293672902914e+01,4.177518476265455405e+01,4.368355882474136820e+01,4.559569652509853199e+01,4.751092027402186346e+01,4.942855204883975517e+01,5.134791386111260181e+01,5.326832822565912551e+01,5.518911863055910061e+01,5.710961000725983894e+01,5.902912919992107987e+01,6.094700543309508589e+01,6.286257077691075068e+01,6.504408142923783487e+01,6.731150555251983292e+01,6.957441968809095556e+01,7.183204804714833358e+01,7.408362017742892647e+01,7.632837147931614652e+01,7.856554371291085204e+01,8.079438549505744049e+01,8.301415278545449894e+01,8.522410936093632472e+01,8.742352727703548965e+01,8.964320922166764660e+01,9.212783985038407764e+01,9.463991702357323277e+01,9.713710436664312908e+01,9.961860333783822341e+01,1.020836271066133776e+02,1.045314009761038250e+02,1.069611627851832054e+02,1.093721632894912403e+02,1.117636665208914906e+02,1.141349501247561733e+02,1.164853056746731710e+02,1.188427495135249075e+02,1.214497742963900464e+02,1.240571397607725856e+02,1.266381670955139924e+02,1.291921142963542479e+02,1.317182562071990048e+02,1.342158847170373974e+02,1.366843089299077008e+02,1.391228553079557742e+02,1.415308677876249988e+02,1.439077078692092471e+02,1.462527546799526874e+02,1.487972463569419119e+02,1.513660869080079863e+02,1.538976689973949590e+02,1.563913496292526020e+02,1.588465052042980403e+02,1.612625314572597972e+02,1.636388433684695940e+02,1.659748750503042913e+02,1.682700796092566407e+02,1.705402243539989229e+02,1.730002404051292899e+02,1.754184532106346524e+02,1.777893410209722731e+02,1.801123719239223817e+02,1.823870336813989752e+02,1.846128334527778918e+02,1.867941480007042685e+02,1.891370012481262108e+02,1.914291683874612033e+02,1.937018180404253371e+02,1.959851575209921464e+02,1.982094242104809894e+02,2.003741594103126715e+02,2.024789240166101365e+02,2.045232981138738353e+02,2.065068805601038378e+02,2.084292885649008156e+02,2.102901572620931745e+02,2.120891392784248808e+02,2.138259042997732422e+02,2.155001386365036069e+02,2.171115447893790815e+02,2.186598410175286062e+02,2.201447609099212457e+02,2.215660529617904047e+02,2.229234801573986715e+02,2.242168195604715493e+02,2.254458619136769073e+02,2.266104112484254074e+02,2.277102845062142933e+02,2.287453111727197665e+02,2.297153329258348151e+02,2.306202032986815311e+02,2.314597873586975254e+02,2.322339614037742592e+02,2.329426126763821969e+02,2.335856390965545870e+02,2.341629490145511738e+02,2.346744609839448401e+02,2.351201035558342767e+02,2.354998150947743341e+02,2.358135436170063883e+02,2.360612466514553205e+02,2.362428911239031777e+02,2.363584532647145124e+02,2.364079185403368228e+02,2.363912816088215152e+02,2.363085462994700947e+02,2.361597256166645309e+02,2.359448417678810017e+02,2.356639262157777068e+02,2.353170197541903690e+02,2.349041726078544627e+02,2.344254445554753374e+02,2.338809050758206070e+02,2.332706335163791209e+02,2.325947192840671107e+02,2.318532620573982967e+02,2.310463720194735799e+02,2.301741701110701399e+02,2.292367883030632356e+02,2.282343698873253004e+02,2.271670697851909324e+02,2.260350548725836859e+02,2.248385043207103138e+02,2.235776099513061865e+02,2.222525766052952179e+02,2.208636225236912196e+02,2.194109797395206272e+02,2.178948944795100999e+02,2.163156275742284436e+02,2.146734548753608181e+02,2.129686676787072201e+02,2.112015731515328980e+02,2.093724947628086568e+02,2.074817727148911501e+02,2.055297643751804344e+02,2.035168447062063706e+02,2.014434066926854712e+02,1.993098617639869872e+02,1.971166402105051532e+02,1.948641915923882664e+02,1.926222690682654672e+02,1.903606555467356145e+02,1.880481723164200787e+02,1.858970337891098552e+02,1.837012053469015598e+02,1.814562992408918944e+02,1.791627970429301797e+02,1.768211994494532462e+02,1.744320265686427831e+02,1.720007880949516448e+02,1.697591462864160974e+02,1.674921575025541927e+02,1.651840417026052990e+02,1.628353328502733746e+02,1.604465837208960863e+02,1.580183660552347931e+02,1.555512706895913766e+02,1.530459076614164644e+02,1.505029062895320635e+02,1.479838601359597590e+02,1.456638513409760378e+02,1.433116278570334430e+02,1.409277966072513095e+02,1.385129821228730407e+02,1.360678265003785725e+02,1.335929893327690365e+02,1.310891476142284091e+02,1.285569956182047520e+02,1.259972447487754863e+02,1.234106233649290090e+02,1.208237967608994978e+02,1.184859918317059737e+02,1.161547027258366569e+02,1.138019000959190947e+02,1.114282795770612324e+02,1.090345509664079060e+02,1.066214379290218517e+02,1.041896776807610649e+02,1.017400206486066594e+02,9.927323010911993606e+01,9.679008180511546300e+01,9.429136354183235369e+01,9.181940132609368277e+01,8.961168731599777004e+01,8.742352727703655546e+01,8.522410936093739053e+01,8.301415278545991328e+01,8.079438549504895661e+01,7.856554371291851169e+01,7.632837147930952426e+01,7.408362017743618821e+01,7.183204804714088709e+01,6.957441968809786204e+01,6.731150555251622336e+01,6.504408142924167180e+01 diff --git a/radin_dailyavg.csv b/radin_dailyavg.csv new file mode 100644 index 0000000..613391f --- /dev/null +++ b/radin_dailyavg.csv @@ -0,0 +1 @@ +6.286257077690969908e+01,6.094700543309402718e+01,5.902912919992142093e+01,5.710961000726149450e+01,5.518911863055893008e+01,5.326832822565901182e+01,5.134791386111293576e+01,4.942855204884011044e+01,4.751092027402129503e+01,4.559569652509747328e+01,4.368355882474016738e+01,4.177518476265366587e+01,3.990277293672754411e+01,3.831239339362680596e+01,3.676791799077611955e+01,3.522837210237122463e+01,3.369430585545826062e+01,3.216626720050992816e+01,3.064480155131714412e+01,2.913045143102182166e+01,2.762375612485181975e+01,2.612525134012183514e+01,2.463546887401730601e+01,2.315493628966762074e+01,2.171288715035241879e+01,2.055301198811046959e+01,1.942724853294134846e+01,1.831007557459968993e+01,1.720188784544575711e+01,1.610307562389593627e+01,1.501402452872270921e+01,1.393511532187155133e+01,1.286672371998985120e+01,1.180922021482493989e+01,1.076296990263483444e+01,9.728332322717138325e+00,8.937502651117220864e+00,8.216318407411662861e+00,7.504082120964423730e+00,6.801036929238594375e+00,6.107421582998850518e+00,5.423470365354699396e+00,4.749413017312665630e+00,4.085474669784134427e+00,3.431875781970633810e+00,2.805127455544981441e+00,2.420881179340481193e+00,2.048159507583730221e+00,1.682114724363046321e+00,1.322866093415396849e+00,9.705299341290687387e-01,6.252196071109451525e-01,2.918960039338246437e-01,1.771454084358596581e-01,6.928392916400623325e-02,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,6.928392916400623325e-02,1.771454084358596581e-01,2.918960039340988133e-01,6.252196071109380471e-01,9.705299341290672954e-01,1.322866093415394184e+00,1.682114724362788749e+00,2.048159507583477090e+00,2.420881179340474088e+00,2.805127455545747051e+00,3.431875781970863404e+00,4.085474669784111335e+00,4.749413017312162033e+00,5.423470365354933875e+00,6.107421582999201348e+00,6.801036929238826190e+00,7.504082120964650215e+00,8.216318407411169034e+00,8.937502651117322117e+00,9.728332322717884395e+00,1.076296990263472075e+01,1.180922021482483153e+01,1.286672371998961495e+01,1.393511532187169344e+01,1.501402452872334514e+01,1.610307562389588654e+01,1.720188784544571448e+01,1.831007557459963309e+01,1.942724853294126675e+01,2.055301198811109487e+01,2.171288715035418804e+01,2.315493628966740047e+01,2.463546887401664165e+01,2.612525134012043893e+01,2.762375612485159948e+01,2.913045143102206325e+01,3.064480155131664674e+01,3.216626720050921762e+01,3.369430585545777035e+01,3.522837210236924221e+01,3.676791799077442846e+01,3.831239339362482355e+01,3.990277293672902914e+01,4.177518476265455405e+01,4.368355882474136820e+01,4.559569652509853199e+01,4.751092027402186346e+01,4.942855204883975517e+01,5.134791386111260181e+01,5.326832822565912551e+01,5.518911863055910061e+01,5.710961000725983894e+01,5.902912919992107987e+01,6.094700543309508589e+01,6.286257077691075068e+01,6.504408142923783487e+01,6.731150555251983292e+01,6.957441968809095556e+01,7.183204804714833358e+01,7.408362017742892647e+01,7.632837147931614652e+01,7.856554371291085204e+01,8.079438549505744049e+01,8.301415278545449894e+01,8.522410936093632472e+01,8.742352727703548965e+01,8.964320922166764660e+01,9.212783985038407764e+01,9.463991702357323277e+01,9.713710436664312908e+01,9.961860333783822341e+01,1.020836271066133776e+02,1.045314009761038250e+02,1.069611627851832054e+02,1.093721632894912403e+02,1.117636665208914906e+02,1.141349501247561733e+02,1.164853056746731710e+02,1.188427495135249075e+02,1.214497742963900464e+02,1.240571397607725856e+02,1.266381670955139924e+02,1.291921142963542479e+02,1.317182562071990048e+02,1.342158847170373974e+02,1.366843089299077008e+02,1.391228553079557742e+02,1.415308677876249988e+02,1.439077078692092471e+02,1.462527546799526874e+02,1.487972463569419119e+02,1.513660869080079863e+02,1.538976689973949590e+02,1.563913496292526020e+02,1.588465052042980403e+02,1.612625314572597972e+02,1.636388433684695940e+02,1.659748750503042913e+02,1.682700796092566407e+02,1.705402243539989229e+02,1.730002404051292899e+02,1.754184532106346524e+02,1.777893410209722731e+02,1.801123719239223817e+02,1.823870336813989752e+02,1.846128334527778918e+02,1.867941480007042685e+02,1.891370012481262108e+02,1.914291683874612033e+02,1.937018180404253371e+02,1.959851575209921464e+02,1.982094242104809894e+02,2.003741594103126715e+02,2.024789240166101365e+02,2.045232981138738353e+02,2.065068805601038378e+02,2.084292885649008156e+02,2.102901572620931745e+02,2.120891392784248808e+02,2.138259042997732422e+02,2.155001386365036069e+02,2.171115447893790815e+02,2.186598410175286062e+02,2.201447609099212457e+02,2.215660529617904047e+02,2.229234801573986715e+02,2.242168195604715493e+02,2.254458619136769073e+02,2.266104112484254074e+02,2.277102845062142933e+02,2.287453111727197665e+02,2.297153329258348151e+02,2.306202032986815311e+02,2.314597873586975254e+02,2.322339614037742592e+02,2.329426126763821969e+02,2.335856390965545870e+02,2.341629490145511738e+02,2.346744609839448401e+02,2.351201035558342767e+02,2.354998150947743341e+02,2.358135436170063883e+02,2.360612466514553205e+02,2.362428911239031777e+02,2.363584532647145124e+02,2.364079185403368228e+02,2.363912816088215152e+02,2.363085462994700947e+02,2.361597256166645309e+02,2.359448417678810017e+02,2.356639262157777068e+02,2.353170197541903690e+02,2.349041726078544627e+02,2.344254445554753374e+02,2.338809050758206070e+02,2.332706335163791209e+02,2.325947192840671107e+02,2.318532620573982967e+02,2.310463720194735799e+02,2.301741701110701399e+02,2.292367883030632356e+02,2.282343698873253004e+02,2.271670697851909324e+02,2.260350548725836859e+02,2.248385043207103138e+02,2.235776099513061865e+02,2.222525766052952179e+02,2.208636225236912196e+02,2.194109797395206272e+02,2.178948944795100999e+02,2.163156275742284436e+02,2.146734548753608181e+02,2.129686676787072201e+02,2.112015731515328980e+02,2.093724947628086568e+02,2.074817727148911501e+02,2.055297643751804344e+02,2.035168447062063706e+02,2.014434066926854712e+02,1.993098617639869872e+02,1.971166402105051532e+02,1.948641915923882664e+02,1.926222690682654672e+02,1.903606555467356145e+02,1.880481723164200787e+02,1.858970337891098552e+02,1.837012053469015598e+02,1.814562992408918944e+02,1.791627970429301797e+02,1.768211994494532462e+02,1.744320265686427831e+02,1.720007880949516448e+02,1.697591462864160974e+02,1.674921575025541927e+02,1.651840417026052990e+02,1.628353328502733746e+02,1.604465837208960863e+02,1.580183660552347931e+02,1.555512706895913766e+02,1.530459076614164644e+02,1.505029062895320635e+02,1.479838601359597590e+02,1.456638513409760378e+02,1.433116278570334430e+02,1.409277966072513095e+02,1.385129821228730407e+02,1.360678265003785725e+02,1.335929893327690365e+02,1.310891476142284091e+02,1.285569956182047520e+02,1.259972447487754863e+02,1.234106233649290090e+02,1.208237967608994978e+02,1.184859918317059737e+02,1.161547027258366569e+02,1.138019000959190947e+02,1.114282795770612324e+02,1.090345509664079060e+02,1.066214379290218517e+02,1.041896776807610649e+02,1.017400206486066594e+02,9.927323010911993606e+01,9.679008180511546300e+01,9.429136354183235369e+01,9.181940132609368277e+01,8.961168731599777004e+01,8.742352727703655546e+01,8.522410936093739053e+01,8.301415278545991328e+01,8.079438549504895661e+01,7.856554371291851169e+01,7.632837147930952426e+01,7.408362017743618821e+01,7.183204804714088709e+01,6.957441968809786204e+01,6.731150555251622336e+01,6.504408142924167180e+01 diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..003b4b1 --- /dev/null +++ b/setup.py @@ -0,0 +1,10 @@ +from setuptools import setup, find_packages + +setup(name='thalake', + version='0.1', + author='Montek Thind', + author_email='montek.thind@colorado.edu', + license='MIT', + description='python version of tha-lake', + packages=find_packages(), +) diff --git a/thalake/bmithalake.py b/thalake/bmithalake.py new file mode 100644 index 0000000..b8e7443 --- /dev/null +++ b/thalake/bmithalake.py @@ -0,0 +1,146 @@ +"""Basic Model Interface (BMI) for the Diffusion model.""" + +import numpy as np +from basic_modeling_interface import Bmi +from thalake import thaLakeModel + +class bmithaLakeModel(Bmi): + """BMI for the Diffusion model.""" + + _name = 'thalake permafrost Component"' + _input_var_names = ('depth_in_meters',) + _output_var_names = ('depth_in_meters',) + + def __init__(self): + """Create a Diffusion model that's ready for initialization.""" + self._model = None + self._values = {} + self._var_units = {} + self._grids = {} + self._grid_type = {} + self.ngrids = 0 + + + def initialize(self, filename=None): + self._model = thaLakeModel() + #self._model.initialize(filename=cfg_file) + #self._model = thaLakeModel(config_file=filename) + #self._values = { + # 'depth_in_meters': self._model.depth, + #} + self._var_units = { + 'depth_in_meters': 'm' + } + self._grids = { + 0: ['plate_surface__temperature'] + } + self._grid_type = { + 0: 'uniform_rectilinear_grid' + } + print ("opening file") + + def update(self): + """Advance model by one time step.""" + self._model.advance() + + def update_frac(self, time_frac): + time_step = self.get_time_step() + self._model.dt = time_frac * time_step + self.update() + self._model.dt = time_step + + def update_until(self, then): + n_steps = (then - self.get_current_time()) / self.get_time_step() + + for _ in range(int(n_steps)): + self.update() + + if (n_steps - int(n_steps)) > 0.0: + self.update_frac(n_steps - int(n_steps)) + + def finalize(self): + """Finalize model.""" + self._model = None + #assert(self._model.status == 'initialized') + + def get_var_type(self, var_name): + return str(self.get_value(var_name).dtype) + + def get_var_units(self, var_name): + return self._var_units[var_name] + + def get_var_nbytes(self, var_name): + return self.get_value(var_name).nbytes + + def get_var_grid(self, var_name): + for grid_id, var_name_list in self._grids.items(): + if var_name in var_name_list: + return grid_id + + def get_grid_rank(self, grid_id): + return len(self.get_grid_shape(grid_id)) + + def get_grid_size(self, grid_id): + size = self.bmi.get_grid_size(grid) + return str(size) + + def get_value_ref(self, var_name): + return self._values[var_name].reshape(-1) + + def get_value(self, var_name): + return self.get_value_ref(var_name).copy() + + + def set_value(self, var_name, src): + val = self.get_value_ref(var_name) + val[:] = src + + def get_component_name(self): + """Name of the component.""" + return self._name + + def get_input_var_names(self): + """Get names of input variables.""" + return self._input_var_names + + def get_output_var_names(self): + """Get names of output variables.""" + return self._output_var_names + + def get_grid_shape(self, grid_id): + """Number of columns and rows of uniform rectilinear grid.""" + return (self._model.ny, self._model.nx) + + def get_grid_spacing(self, grid_id): + """Spacing of columns and rows of uniform rectilinear grid.""" + assert_true(grid_id < self.ngrids) + return np.array([1, 1], dtype='float32') + #return (self._model.dy, self._model.dx) + + def get_grid_origin(self, grid_id): + """Origin of uniform rectilinear grid.""" + return (0.0, 0.0) + + def get_grid_type(self, grid_id): + """Type of grid.""" + return self._grid_type[grid_id] + + def get_start_time(self): + """Start time of model.""" + return 0.0 + + def get_end_time(self): + """End time of model.""" + return np.finfo('d').max + + def get_current_time(self): + """Current time of model.""" + return self._model.time + + def get_time_step(self): + """Time step of model.""" + return self._model.dt + + def get_time_units(self): + """Time units of model.""" + return '-' diff --git a/thalake/radin_dailyavg.csv b/thalake/radin_dailyavg.csv new file mode 100644 index 0000000..613391f --- /dev/null +++ b/thalake/radin_dailyavg.csv @@ -0,0 +1 @@ +6.286257077690969908e+01,6.094700543309402718e+01,5.902912919992142093e+01,5.710961000726149450e+01,5.518911863055893008e+01,5.326832822565901182e+01,5.134791386111293576e+01,4.942855204884011044e+01,4.751092027402129503e+01,4.559569652509747328e+01,4.368355882474016738e+01,4.177518476265366587e+01,3.990277293672754411e+01,3.831239339362680596e+01,3.676791799077611955e+01,3.522837210237122463e+01,3.369430585545826062e+01,3.216626720050992816e+01,3.064480155131714412e+01,2.913045143102182166e+01,2.762375612485181975e+01,2.612525134012183514e+01,2.463546887401730601e+01,2.315493628966762074e+01,2.171288715035241879e+01,2.055301198811046959e+01,1.942724853294134846e+01,1.831007557459968993e+01,1.720188784544575711e+01,1.610307562389593627e+01,1.501402452872270921e+01,1.393511532187155133e+01,1.286672371998985120e+01,1.180922021482493989e+01,1.076296990263483444e+01,9.728332322717138325e+00,8.937502651117220864e+00,8.216318407411662861e+00,7.504082120964423730e+00,6.801036929238594375e+00,6.107421582998850518e+00,5.423470365354699396e+00,4.749413017312665630e+00,4.085474669784134427e+00,3.431875781970633810e+00,2.805127455544981441e+00,2.420881179340481193e+00,2.048159507583730221e+00,1.682114724363046321e+00,1.322866093415396849e+00,9.705299341290687387e-01,6.252196071109451525e-01,2.918960039338246437e-01,1.771454084358596581e-01,6.928392916400623325e-02,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,6.928392916400623325e-02,1.771454084358596581e-01,2.918960039340988133e-01,6.252196071109380471e-01,9.705299341290672954e-01,1.322866093415394184e+00,1.682114724362788749e+00,2.048159507583477090e+00,2.420881179340474088e+00,2.805127455545747051e+00,3.431875781970863404e+00,4.085474669784111335e+00,4.749413017312162033e+00,5.423470365354933875e+00,6.107421582999201348e+00,6.801036929238826190e+00,7.504082120964650215e+00,8.216318407411169034e+00,8.937502651117322117e+00,9.728332322717884395e+00,1.076296990263472075e+01,1.180922021482483153e+01,1.286672371998961495e+01,1.393511532187169344e+01,1.501402452872334514e+01,1.610307562389588654e+01,1.720188784544571448e+01,1.831007557459963309e+01,1.942724853294126675e+01,2.055301198811109487e+01,2.171288715035418804e+01,2.315493628966740047e+01,2.463546887401664165e+01,2.612525134012043893e+01,2.762375612485159948e+01,2.913045143102206325e+01,3.064480155131664674e+01,3.216626720050921762e+01,3.369430585545777035e+01,3.522837210236924221e+01,3.676791799077442846e+01,3.831239339362482355e+01,3.990277293672902914e+01,4.177518476265455405e+01,4.368355882474136820e+01,4.559569652509853199e+01,4.751092027402186346e+01,4.942855204883975517e+01,5.134791386111260181e+01,5.326832822565912551e+01,5.518911863055910061e+01,5.710961000725983894e+01,5.902912919992107987e+01,6.094700543309508589e+01,6.286257077691075068e+01,6.504408142923783487e+01,6.731150555251983292e+01,6.957441968809095556e+01,7.183204804714833358e+01,7.408362017742892647e+01,7.632837147931614652e+01,7.856554371291085204e+01,8.079438549505744049e+01,8.301415278545449894e+01,8.522410936093632472e+01,8.742352727703548965e+01,8.964320922166764660e+01,9.212783985038407764e+01,9.463991702357323277e+01,9.713710436664312908e+01,9.961860333783822341e+01,1.020836271066133776e+02,1.045314009761038250e+02,1.069611627851832054e+02,1.093721632894912403e+02,1.117636665208914906e+02,1.141349501247561733e+02,1.164853056746731710e+02,1.188427495135249075e+02,1.214497742963900464e+02,1.240571397607725856e+02,1.266381670955139924e+02,1.291921142963542479e+02,1.317182562071990048e+02,1.342158847170373974e+02,1.366843089299077008e+02,1.391228553079557742e+02,1.415308677876249988e+02,1.439077078692092471e+02,1.462527546799526874e+02,1.487972463569419119e+02,1.513660869080079863e+02,1.538976689973949590e+02,1.563913496292526020e+02,1.588465052042980403e+02,1.612625314572597972e+02,1.636388433684695940e+02,1.659748750503042913e+02,1.682700796092566407e+02,1.705402243539989229e+02,1.730002404051292899e+02,1.754184532106346524e+02,1.777893410209722731e+02,1.801123719239223817e+02,1.823870336813989752e+02,1.846128334527778918e+02,1.867941480007042685e+02,1.891370012481262108e+02,1.914291683874612033e+02,1.937018180404253371e+02,1.959851575209921464e+02,1.982094242104809894e+02,2.003741594103126715e+02,2.024789240166101365e+02,2.045232981138738353e+02,2.065068805601038378e+02,2.084292885649008156e+02,2.102901572620931745e+02,2.120891392784248808e+02,2.138259042997732422e+02,2.155001386365036069e+02,2.171115447893790815e+02,2.186598410175286062e+02,2.201447609099212457e+02,2.215660529617904047e+02,2.229234801573986715e+02,2.242168195604715493e+02,2.254458619136769073e+02,2.266104112484254074e+02,2.277102845062142933e+02,2.287453111727197665e+02,2.297153329258348151e+02,2.306202032986815311e+02,2.314597873586975254e+02,2.322339614037742592e+02,2.329426126763821969e+02,2.335856390965545870e+02,2.341629490145511738e+02,2.346744609839448401e+02,2.351201035558342767e+02,2.354998150947743341e+02,2.358135436170063883e+02,2.360612466514553205e+02,2.362428911239031777e+02,2.363584532647145124e+02,2.364079185403368228e+02,2.363912816088215152e+02,2.363085462994700947e+02,2.361597256166645309e+02,2.359448417678810017e+02,2.356639262157777068e+02,2.353170197541903690e+02,2.349041726078544627e+02,2.344254445554753374e+02,2.338809050758206070e+02,2.332706335163791209e+02,2.325947192840671107e+02,2.318532620573982967e+02,2.310463720194735799e+02,2.301741701110701399e+02,2.292367883030632356e+02,2.282343698873253004e+02,2.271670697851909324e+02,2.260350548725836859e+02,2.248385043207103138e+02,2.235776099513061865e+02,2.222525766052952179e+02,2.208636225236912196e+02,2.194109797395206272e+02,2.178948944795100999e+02,2.163156275742284436e+02,2.146734548753608181e+02,2.129686676787072201e+02,2.112015731515328980e+02,2.093724947628086568e+02,2.074817727148911501e+02,2.055297643751804344e+02,2.035168447062063706e+02,2.014434066926854712e+02,1.993098617639869872e+02,1.971166402105051532e+02,1.948641915923882664e+02,1.926222690682654672e+02,1.903606555467356145e+02,1.880481723164200787e+02,1.858970337891098552e+02,1.837012053469015598e+02,1.814562992408918944e+02,1.791627970429301797e+02,1.768211994494532462e+02,1.744320265686427831e+02,1.720007880949516448e+02,1.697591462864160974e+02,1.674921575025541927e+02,1.651840417026052990e+02,1.628353328502733746e+02,1.604465837208960863e+02,1.580183660552347931e+02,1.555512706895913766e+02,1.530459076614164644e+02,1.505029062895320635e+02,1.479838601359597590e+02,1.456638513409760378e+02,1.433116278570334430e+02,1.409277966072513095e+02,1.385129821228730407e+02,1.360678265003785725e+02,1.335929893327690365e+02,1.310891476142284091e+02,1.285569956182047520e+02,1.259972447487754863e+02,1.234106233649290090e+02,1.208237967608994978e+02,1.184859918317059737e+02,1.161547027258366569e+02,1.138019000959190947e+02,1.114282795770612324e+02,1.090345509664079060e+02,1.066214379290218517e+02,1.041896776807610649e+02,1.017400206486066594e+02,9.927323010911993606e+01,9.679008180511546300e+01,9.429136354183235369e+01,9.181940132609368277e+01,8.961168731599777004e+01,8.742352727703655546e+01,8.522410936093739053e+01,8.301415278545991328e+01,8.079438549504895661e+01,7.856554371291851169e+01,7.632837147930952426e+01,7.408362017743618821e+01,7.183204804714088709e+01,6.957441968809786204e+01,6.731150555251622336e+01,6.504408142924167180e+01 diff --git a/thalake/thalake.py b/thalake/thalake.py new file mode 100644 index 0000000..811f4a3 --- /dev/null +++ b/thalake/thalake.py @@ -0,0 +1,553 @@ +# Thaw Lake Model-1D + +# This model a 1-D numerical model of permafrost and subsidence processes, it is a Python version. +# It aims to investigate the subsurface thermal impact of thaw lakes of various depths, +# and to evaluate how this impact might change in a warming climate. + +# Key paper: Matell, N., Anderson, R.S., Overeem, I., Wobus, C., Urban, F., +# Clow, G., 2013 Modeling the subsurface thermal impact of Arctic thaw lakes in a warming climate. Computers and Geosciences. + +# Originates from earlier code by Nora Matell and co-authors. +# Copyright to Python version (C) <2017> + +# Developer can be contacted by irina.overeem@colorado.edu + +# Dr. Irina Overeem +# CSDMS Community Surface Dynamics Modeling System +# INSTAAR, University of Colorado at Boulder +# PO Box 450, 80309-0450 +# Boulder, CO, USA + + +# This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. +# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. +# You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +# This model a 1-D numerical model of permafrost and subsidence processes. +# It aims to investigate the subsurface thermal impact of thaw lakes of various depths, +# and to evaluate how this impact might change in a warming climate. + +# The model was designed for the Alaskan Arctic Coastal Plain +# The model uses average and amplitude of temperature from observed permafrost temperatures at 5cm in the subsurface at Drew Point data (USGS) +# Clow, G.D., 2008a. Continued permafrost warming in northern Alaska, 2008 update. NOAA/ESRL Global Monitoring Annual Conference, Boulder. + +# The model uses observed heat flow for thermal gradient from: +# Lachenbruch, A.H., Sass, J.H., Marshall, B.V., Moses, T.H., Jr., 1982. Temperatures, heat flow, and the geothermal regime at Prudhoe Bay, Alaska. Journal of Geophysical Research 87, 9301-9316. + + +# in 1-D, subsurface temperature changes through time and space of a semi-infinite +# half-space, then plots result. The model code incorporates phase change by using +# the apparent heat capacity scheme for temperatures within the +# "phase-change envelope". + +# lake-permafrost model - lake freezes and thaws, permafrost changes temperature, +# when permafrost thaw it subsides, assuming that if excess ice melts all +# the water leaves the permafrost and thus the subsurface volume decreases. + +# boundary condition of ice bottom and lake water top = Tf +# if TsTf, top +# boundary condition --> Tf and excess energy is used to melt ice. When +# ice melted from the top, all other boundary just moved up by equivelant +# amount, if moved up a full control volume then bottom control volume +# added that is same temp as (end-1) control volume. Assumes that when ice +# free, lake is completely mixed. (ie Liston & Hall) + +import math +import numpy as np +import csv +import operator +import copy +import scipy.io as sio +import matplotlib.pyplot as plt +#from permamodel.utils import model_input +#from permamodel.components import perma_base + +class thaLakeModel(): + print ("hello maybe this is working") + def __init__(self): + # set up 1-d grid + self.dz0 = 0.05; # cell size in top 10 m + self.dz2 = 1; # cell size below 10 m + self.deptht = 100; # model depth [m] + self.z_range1 = np.arange(0,10.05,self.dz0); + self.z_range2 = np.arange(11,self.deptht+1,self.dz2); + self.iceinit = 0.0001; # initial ice depth [m] + self.periodyear = 3600*24*365.; # period (1 year) [s] + self.dt = 3600*24; # length of timestep [s] + self.years = 1000; # simulation duration + self.depthsubside = 2.0; # seed lake depth, has to match with ln 79, because the lake contains water + self.depthtalik = 0; # usually unkown, so the model will trend toward equilibrium in the first 100-200 years of a simulation + self.yearnum = 0; + self.kelvin = 273.15; + self.rhoi = 917; # ice density [kg/m3] + self.ci = 2108; # ice specific heat capacity [J/kg/K] + self.ki = 2.18; # ice thermal conductivity [J/s/m/K] + self.rhow = 1000; # water density [kg/m3] + self.cw = 4210; # water specific heat capacity [J/kg/K] + self.kw = 0.58; # water thermal conductivity [J/s/m/K] + self.rhorock = 1200; # mineral soil density [kg/m3] + self.crock = 1000; # mineral soil specific heat capacity [J/kg/degree C] + self.krock = 1.5; # rock thermal conductivity [J/s/m/degree C] + self.totalice = 0.65; # ice content of original substrate, from bulk samples in top 5m at Drew Point + self.excessice = 0.3; # excess ice content - gone once melted + self.L = 334*1000; # latent heat of fusion for water [J/kg] + self.hw = 0.56; # convective transfer coefficient =[J/s/m2/K] + self.Kh = 0.0; # turbulent diffusivity for heat + self.pcenv = 1; # width of phase change envelope [m] + self.warming = 5; + self.Tamp = 17.5; # amplitude of air temperature fluctuations [K] + self.q = 0.056; # mantle heat flow [J/m2/s] + self.count = 1; + self.daynum = 0; + self.numplots = 50; + + def initialize(self): + self.z = np.append(self.z_range1, self.z_range2); # vertical grid (distance below surface) [m] + self.dz = np.diff(self.z); + self.zbetween = self.z[1:] + np.diff(self.z)/2; + self.dzbetween = np.diff(self.zbetween); + self.numcells = len(self.z) # number of grid cells + self.celltype = np.int16(1*np.ones(self.numcells)); # substrate type - 1=regular soil with excess ice, 2=compacted soil without excess ice, 8=water, 9=ice + self.celltype[np.where(self.z<6.7)] =2; # depth of subsided cells, this is a function of the initial lake depth (seed lake depth/excess ice content) + self.celltype[np.where(self.z<2.0)] = 8; # seed lake depth + self.celltype[np.where(self.z<=self.iceinit)] = 9; # ice thickness + self.celltype[np.where(self.z>=50)] = 2; # underlying cells with no excess ice + self.thawed = np.zeros(self.numcells); + self.thawedspec = np.zeros(self.numcells); + + self.deptht = max(self.z); + self.depthi = self.iceinit; + + # set time steps and simulation duration + + self.nt = int(((self.periodyear/self.dt)*self.years)+1); # number of timesteps + self.t = np.arange(0,self.dt*self.nt+1,self.dt); # time at each timestep [s] + self.tday = [x/(3600*24) for x in self.t]; + self.tyear = [x/(3600*24*365) for x in self.t]; + self.tyearsonly = np.arange(0,self.years); + + # declare zarrays + + self.depthin = np.zeros(self.nt); + self.depthwn = np.zeros(self.nt); + self.depthsubsiden = np.zeros(self.nt+1); + self.depthtalikn = np.zeros(self.nt+1); + self.icethickness = np.zeros(self.nt+1); + self.icegrowth = np.zeros(self.nt+1); + self.soilthickness = np.zeros(self.nt+1); + self.waterthickness = np.zeros(self.nt+1); + self.Tbot = np.zeros(self.nt+1); + self.T3m = np.zeros(self.nt+1); + self.T5m = np.zeros(self.nt+1); + self.T10m = np.zeros(self.nt+1); + self.T25m = np.zeros(self.nt+1); + self.T50m = np.zeros(self.nt+1); + self.T100m = np.zeros(self.nt+1); + self.Tbotavgann = np.zeros(self.years); + self.T3mavgann = np.zeros(self.years); + self.T5mavgann = np.zeros(self.years); + self.T10mavgann = np.zeros(self.years); + self.T25mavgann = np.zeros(self.years); + self.T50mavgann = np.zeros(self.years); + self.T100mavgann = np.zeros(self.years); + self.Trecord = np.zeros((365,self.numcells)); + self.watts = np.zeros(np.size(self.t)); + + + # define constants + # Model input and boundary conditions (compiled from field and modeling studies by West & Plug, 2008; Ling and Zhang, 2003;..... + # Ling & Zhang, 2004; Zhou and Huang, 2004; Romanovsky and Osterkamp, 2000; Hinzman, 1998; and French, 2007). + + self.kappai = self.ki/(self.ci*self.rhoi); + self.kappaw = self.kw/(self.cw*self.rhow); + + self.porewater = self.totalice-self.excessice; # porewater ice - still there if refreezes + self.W = self.totalice; # total water content percent of material by mass + self.Wu = self.totalice*0.05; # unfrozen percent water content of material by mass at + #temperature T = Tpc-pcenv + self.Ws = self.porewater; + self.Wus = self.porewater*0.05; + self.totalice2 = self.totalice-self.Wu; + self.porewater2 = self.porewater-self.Wus; + + self.cf = (self.ci*self.totalice2)+(self.cw*self.Wu)+(self.crock*(1-self.totalice)); + self.cu = (self.cw*self.totalice)+(self.crock*(1-self.totalice)); + self.cfs = (self.ci*self.porewater2)+(self.cw*self.Wu)+(self.crock*(1-self.porewater)); + self.cus = (self.cw*self.porewater)+(self.crock*(1-self.porewater)); + self.rhof = (self.rhoi*self.totalice2)+(self.rhow*self.Wu)+(self.rhorock*(1-self.totalice)); + self.rhou = (self.rhow*self.totalice)+(self.rhorock*(1-self.totalice)); + self.rhofs = (self.rhoi*self.porewater2)+(self.rhow*self.Wus)+(self.rhorock*(1-self.porewater)); + self.rhous = (self.rhow*self.porewater)+(self.rhorock*(1-self.porewater)); + self.Cf = self.cf*self.rhof; # frozen volumetric heat capacity [J/m3/degree C] + self.Cu = self.cu*self.rhou; # thawed volumetric heat capacity [J/m3/degree C] + self.Cfs = self.cfs*self.rhofs; + self.Cus = self.cus*self.rhous; + self.kf = pow(self.ki,self.totalice2)*pow(self.kw,self.Wu)*pow(self.krock,(1-self.totalice)); + self.ku = pow(self.kw,self.totalice)*pow(self.krock,(1-self.totalice)); + self.kfs = pow(self.ki,self.porewater2)*pow(self.kw,self.Wus)*pow(self.krock,(1-self.porewater)); + self.kus = pow(self.kw,self.porewater)*pow(self.krock,(1-self.porewater)); + + self.Tpc = self.kelvin+0; # freezing temperature [K] + + # set up surface and initial temperatures, including geothermal gradient + # theoretical warming scenario is set here + self.Tbar = self.kelvin-11; # MAAT [K] + + self.Tbarplus = self.Tbar+(float(self.warming)/36500)*np.arange(1,36501); + + + self.Tsurface = [self.Tbar-self.Tamp*math.sin(2*math.pi*x/self.periodyear) for x in self.t]; + self.Tsurface[182499:218999] = self.Tbarplus[0:]-self.Tamp*np.sin(2.*np.pi)*self.t[182499:218999]/self.periodyear; + self.Tsurface[218999:] = self.Tbar+self.warming-self.Tamp*np.sin(2*np.pi)*self.t[218999:]/self.periodyear; + + self.dTdzbase = self.q/self.kf; + self.dTdzbase2 = self.q/self.kfs; + self.Tgrad_first_half = [np.multiply(self.dTdzbase,self.z[np.where(self.z<=10)])] + self.Tgrad_second_half= [np.multiply(self.dTdzbase2,self.z[np.where(self.z>10)])]; + self.Tgrad = np.append(self.Tgrad_first_half, self.Tgrad_second_half); + + self.Twi = self.Tsurface[0]; # start up water temperature [K] + self.Tsi = self.Tbar; # start up permafrost temperature [K] + self.Tinit = np.ones((np.size(self.z))); # initial temperature grid + self.Tinitlake_first_half = [self.Tpc*self.Tinit[np.where(self.celltype == 9)]] + self.Tinitlake_second_half= [self.Twi*self.Tinit[np.where(self.celltype == 8)]]; + self.Tinitlake = np.append(self.Tinitlake_first_half, self.Tinitlake_second_half); + self.water = len(self.Tinitlake); + for n in range(0,self.numcells): + if (self.celltype is 9): + self.Tinit[n] = self.Tpc; + elif (self.celltype is 8): + self.Tinit[n] = self.Twi; + else: + self.Tinit[n] = self.Tsi+self.Tgrad[n]; + + # load daily average radiation, function from Drew Point meteorological data + self.mat_file = np.genfromtxt('radin_dailyavg.csv', delimiter=','); + + #Tinitperm = [Tsi*Tinit(find(celltype==1))+Tgrad(find(celltype==1)),Tsi*Tinit(find(celltype==2))+Tgrad(find(celltype==2))]; + #Tinit = [Tinitlake, Tinitperm]; + + self.T = self.Tinit; + self.zstarthaw = math.sqrt((self.ku/self.Cu)*self.periodyear/(math.pi)); + self.zstarfrozen = math.sqrt((self.kfs/self.Cf)*self.periodyear/(math.pi)); + self.Tbase = self.Tbar+(self.dTdzbase2*self.deptht)+np.multiply(self.Tamp,np.exp(-self.deptht/self.zstarthaw))*np.sin((2*math.pi*np.divide(self.t,self.periodyear))-(self.deptht/self.zstarthaw)); + self.Tright = -self.kelvin+self.Tbar+np.multiply(self.dTdzbase2,self.z)+np.multiply(self.Tamp,np.exp(-self.z/self.zstarthaw)); #outer edges of the funnel + self.Trightf = -self.kelvin+self.Tbar+np.multiply(self.dTdzbase2,self.z)+np.multiply(self.Tamp,np.exp(-self.z/self.zstarfrozen)); #outer edges of the funnel + self.Tleftf = -self.kelvin+self.Tbar+np.multiply(self.dTdzbase2,self.z)-np.multiply(self.Tamp,np.exp(-self.z/self.zstarfrozen)); + self.Tleft = -self.kelvin+self.Tbar+np.multiply(self.dTdzbase2,self.z)-np.multiply(self.Tamp,np.exp(-self.z/self.zstarthaw)); + #Ts0 = Tsurface[1]; + + + self.countprint = self.dt*self.nt/self.numplots; + self.nplot=0; + self.Loopcount = 0; + + def updateModel(self): + for n in range(0,5): + self.Loopcount = self.Loopcount+1; # this is a total loop count for debugging a.o. + print (self.Loopcount); + self.thawedspec = np.zeros(self.numcells); + self.thawed = np.zeros(self.numcells); + self.count = self.count+1; + time = (n+1)*self.dt; # time into run [s] + error = 1; + Told = self.T; # remember temperatures from last timestep... + self.Ts0 = self.Tsurface[n]; + self.T[0]= self.Ts0; + + if self.depthi>0 and self.Tsurface[n]>self.Tpc: + self.Ts0 = self.Tpc; + + self.water = -1; + self.ice = -1; + self.subsided = -1; + self.regsoil = -1; + + for i in range(0,self.numcells): + if self.celltype[i] == 8: + self.water = self.water+1; + elif (self.celltype[i] == 9): + self.ice = self.ice+1; + elif (self.celltype[i] == 2): + self.subsided = self.subsided+1; + elif (self.celltype[i] == 1): + self.regsoil = self.regsoil+1; + + self.day = (n % 365); + self.solarrad = self.mat_file[self.day]; + + if self.celltype[0] == 8: + print ('in if loop') + Twater = T[0:water]; # T(find(celltype==8)); + Tmix = np.mean(Twater); #mix again + sextinc = 0.6; # solar extinction constant + albedo = 0.06; + qrad = (1-albedo)*solarrad*np.exp(-sextinc*zbetween[0:water]); + watts[count] = solarrad; + error = 1; + count = 0 + while error>0.0001: + self.Tcalc = [0] + self.Tit = Twater; + for i in range (1,self.water-1): + self.coeff1 = self.kw/self.dz[i]; + self.coeff2 = self.kw/self.dz[i-1]; + self.coeff3 = self.cw*self.rhow*self.dzbetween[i-1]/self.dt; + self.coeff4 = self.coeff1+self.coeff2+self.coeff3; + self.Tcalc = np.append(self.Tcalc, ((self.coeff1*self.T[i+1]+self.coeff2*T[i-1]+coeff3*Told[i]+qrad[i-1]-qrad[i])/coeff4)); + self.Twater_Partial = np.append(self.Ts0,self.Tcalc[1:]); + self.Twater = np.append(self.Twater_Partial, self.Tcalc[-1]); + error = max(abs(self.Twater-self.Tit)); + self.Twater2 = self.Twater[0:-1]; + self.Tmix = np.mean(self.Twater2); # thoroughly mix lake + self.T[1:self.water] = [self.Tmix*np.ones(np.size(x))for x in self.Twater2]; + # compute temperatures in ice layer if completely frozen to bottom + + elif self.water == 0: + print ('in elif loop') + self.Tice = (self.T[np.where(self.celltype==9)]); + #print 'Starting in the loop ', Tice + while error>0.0001: + self.Tcalc = [0] + Tit = copy.deepcopy(Tice); + for i in range(1,self.ice): + self.coeff1 = self.ki/self.dz[i]; + self.coeff2 = self.ki/self.dz[i-1]; + self.coeff3 = self.ci*self.rhoi*self.dzbetween[i-1]/self.dt; + self.coeff4 = self.coeff1+self.coeff2+self.coeff3; + self.Tcalc = np.append(self.Tcalc, (self.coeff1*self.T[i+1]+self.coeff2*self.T[i-1]+self.coeff3*self.Told[i])/self.coeff4); + self.Tice = np.append(self.Ts0,self.Tcalc[1:]); + #print "The final___ TICE ", Tice; + error = max(abs(self.Tice-self.Tit)); + + self.T[np.where(self.celltype==9)] = self.Tice; + #compute temperatures in ice and water layers if both exist + + else: + print ("else loop") + if self.Tsurface[n]>=self.Tpc: + self.T[0:self.ice] = self.Tpc; + print ("Else -If") + elif self.ice<3: + self.T[0] = self.Tsurface[n]; + self.T[self.ice-1] = self.Tpc; + print ("Maybe Elif it is") + else: + print ('Maybe in this loop') + self.Ti = self.T[0:self.ice]; + error = 1; + while error>0.0001: + self.Tcalci = [0] + #Tcalci= np.zeros(ice); + self.Titi = copy.deepcopy(Ti); + for i in range(1,self.ice-1): + coeff1 = self.ki/dz[i]; + coeff2 = self.ki/self.dz[i-1]; + self.coeff3 = self.ci*self.rhoi*self.dzbetween[i-1]/self.dt; + self.coeff4 = self.coeff1+self.coeff2+self.coeff3; + self.Tcalci = np.append(self.Tcalci, (self.coeff1*self.T[i+1]+self.coeff2*self.T[i-1]+self.coeff3*self.Told[i])/self.coeff4); + self.Ti_partial=np.append(self.Ts0,self.Tcalci[1:]); + self.Ti = np.append(self.Ti_partial,self.Tpc); + error = max(abs(self.Ti-self.Titi)); + self.T[0:self.ice] = self.Ti; + error = 1; + if self.water==1: + self.T[np.where(self.celltype==8)] = self.Tpc; + else: + error = 1; + while error>0.0001: + Tcalcw= np.zeros(self.ice); + self.Tit = copy.deepcopy(self.T); + for i in range (self.ice,self.ice+self.water+1): + coeff1 = self.kw/self.dz[i]; + coeff2 = self.kw/self.dz[i-1]; + coeff3 = self.cw*self.rhow*self.dzbetween[i-1]/self.dt; + coeff4 = coeff1+coeff2+coeff3; + Tcalcw = np.append(Tcalcw, (coeff1*self.T[i+1]+coeff2*self.T[i-1]+coeff3*Told[i])/coeff4); + Tw = np.append(self.Tpc, Tcalcw[self.ice:]); + self.T[self.ice:self.ice+self.water+2] = Tw; + error = max(abs(self.T-self.Tit)); + + # Calculate temperatures in underlying permafrost + #clear 'z2' 'zbetween2' 'Told2' 'T2' 'celltype2' + #clear 'z2' 'zbetween2' 'Told2' 'T2' 'celltype2' + z2 = self.z[self.ice+self.water-1:]; + zbetween2 = self.zbetween[self.ice+self.water-1:]; + Told2 = Told[self.ice+self.water-1:]; + T2 = self.T[self.ice+self.water-1:]; + numcells2 = len(T2); + celltype2 = self.celltype[self.ice+self.water-1:]; + C = np.zeros(numcells2); + k = np.zeros(numcells2); + Tcalc = 0; + for i in range (0,numcells2): + if celltype2[i] == 8: + C[i] = self.cw*self.rhow; + k[i] = self.kw; + elif celltype2[i] == 9: + C[i] = self.ci*self.rhoi; + k[i] = self.ki; + elif celltype2[i]==1 and T2[i] < (self.Tpc-self.pcenv): + C[i] = self.Cf; + k[i] = self.kf; + #print('real permafrost, frozen') + elif celltype2[i]==1 and T2[i]>=(self.Tpc-self.pcenv) and T2[i]<=self.Tpc: + C[i] = self.Cf+self.L*self.rhof*((self.W-self.Wu)/self.pcenv); + k[i] = self.kf+((self.ku-self.kf)/self.pcenv)*(T2[i]-(self.Tpc-self.pcenv)); + print('real permafrost, thawing') + elif celltype2[i]==1 and T2[i]>self.Tpc: + C[i] = self.Cu; + k[i] = self.ku; + print('real permafrost, thawed') + elif celltype2[i]==2 and T2[i]<(self.Tpc-self.pcenv): + C[i] = self.Cfs; + k[i] = self.kfs; + elif celltype2[i]==2 and T2[i]>=(self.Tpc-self.pcenv) and T2[i]<=self.Tpc: + C[i] = self.Cfs+self.L*self.rhofs*((self.Ws-self.Wus)/self.pcenv); + k[i] = self.kfs+((self.kus-self.kfs)/self.pcenv)*(T2[i]-(self.Tpc-self.pcenv)); + elif celltype2[i]==2 and T2[i]>self.Tpc: + C[i] = self.Cus; + k[i] = self.kus; + + kbetween = k[1:]-(np.diff(k)/2); + dz2 = np.diff(z2); + dzbetween2 = np.diff(zbetween2); + error = 1; + #Tcalcperm= np.zeros(251) + while error>0.0001: + Tcalcperm = [0] + Tit = copy.deepcopy(T2); + for i in range(1,numcells2-1): + coeff1 = kbetween[i]/dz2[i]; + coeff2 = kbetween[i-1]/dz2[i-1]; + coeff3 = C[i-1]*dzbetween2[i-1]/self.dt; + coeff4 = coeff1+coeff2+coeff3; + Tcalcperm = np.append(Tcalcperm, ((coeff1*T2[i+1]+coeff2*T2[i-1]+coeff3*Told2[i])/coeff4)); + #Tcalcperm[i] = (coeff1*T2[i+1]+coeff2*T2[i-1]+coeff3*Told2[i])/coeff4; + T2_partial = np.append(self.T[self.ice+self.water],Tcalcperm[1:]) + T2 = np.append(Tcalcperm,self.Tbase[n]); + error = max(abs(T2-Tit)); + self.T[self.ice+self.water-1:] = T2; + + # subside permafrost + print (len(range(self.water+self.ice,self.numcells))); + + for i in range (self.water+self.ice, self.numcells): + if self.T[i]self.Tpc and celltype[i] == 1: + thawedspec[i] = 1; + thawed[i] = 1; + celltype[i] = 2; + #print 'is it here' + elif T[i]>Tpc and celltype[i] == 2: + thawed[i] = 1; + #print 'where is it' + + ##[Tbar-Tamp*math.sin(2*math.pi*x/periodyear) for x in t]; + numthawedspec = np.sum(self.thawedspec)+1; + depthsubsidenew = self.z[numthawedspec]*self.excessice; + #depthsubsidenew = self.z*self.excessice; + depthsubside = self.depthsubside+depthsubsidenew; + self.depthsubsiden[n] = self.depthsubside; + cellsubside = np.int16(math.floor((self.depthsubside)/self.dz0)); + maxcellsubside = cellsubside; + #print range(self.ice,maxcellsubside) + if self.depthsubside>0: #depthsubsidenew>0 + for i in range(self.ice,maxcellsubside): + self.celltype[i] = 8; + #print 'p' + + if (self.Tsurface[n]self.ice+self.water: + icecheckdepth = self.ice+self.water-1; + Twater = self.T[icecheckdepth]; + if self.water==1: + Twater = T[np.where(self.celltype==8)]; + Twater = Twater[0]; + if Twater0: + ddepthidt_partial = ((self.Tpc-self.Ts0)*pow((depthi/self.ki),(-1))-self.hw*(Twater-self.Tpc)); + ddepthidt = np.divide(ddepthidt_partial, self.rhoi*self.L) + ddepthi = ddepthidt*self.dt; + #### ADJUSTMENT HERE + if ddepthi>0.05: + ddepthi = 0.05; + print('had to adjust'); + self.icegrowth[n] = ddepthi; + depthi = depthi+ddepthi; + if depthi<0.00001: + depthi = 0; + depthw = self.depthsubside-depthi; + + if depthi>self.depthsubside: + depthi = self.depthsubside; + depthw = 0; + + # check if Ts>Tpc (if there is ice there already) + elif self.Tsurface[n]>Tpc and depthi>0: + Qm = (self.Tsurface[n]-Tpc)*(pow((self.dz0/self.ki),(-1))); + dMdt = Qm/(L*rhoi); + dM = dMdt*self.dt; + icegrowth[n] = -dM; + depthi = depthi-dM; + if depthi<0: + depthi = 0; + maxcelli = np.int16(round(depthi/self.dz0)); + self.celltype[0:self.ice+self.water] = 8; + self.celltype[0:maxcelli] = 9; + self.icethickness[n] = depthi; + self.soilthickness[n] = self.deptht-self.depthsubside; + self.waterthickness[n] = self.deptht-self.icethickness[n]-self.soilthickness[n]; + + self.Tbot[n] = self.T[self.ice+self.water]-self.kelvin; + self.T3m[n] = self.T[59]-self.kelvin; + self.T5m[n] = self.T[99]-self.kelvin; + self.T10m[n] = self.T[199]-self.kelvin; + self.T25m[n] = self.T[214]-self.kelvin; + self.T50m[n] = self.T[239]-self.kelvin; + self.T100m[n] = self.T[289]-self.kelvin; + if (self.tday[n]%365)==0 and self.tday[n]>364: + yearnum = yearnum+1; + Tbotavgann[yearnum] = np.mean(self.Tbot[n-364:n]); + T3mavgann[yearnum] = np.mean(T3m[n-364:n]); + T5mavgann[yearnum] = np.mean(T5m[n-364:n]); + T10mavgann[yearnum] = np.mean(T10m[n-364:n]); + T25mavgann[yearnum] = np.mean(T25m[n-364:n]); + T50mavgann[yearnum] = np.mean(T50m[n-364:n]); + T100mavgann[yearnum] = np.mean(T100m[n-364:n]); + + if n>=self.nt-364: + daynum = daynum+1; + Trecord[daynum, 0:] = T; + + if((n % (365*50))==0): + # nplot = nplot+1 + fig_one = plt.figure(1) + # ice + self.icethickness[n] + # water + self.waterthickness[n] + + lakedepthx = np.arange(-20,25,1); + lakedepth = self.icethickness[n]+self.waterthickness[n]*np.ones(np.size(lakedepthx)); + Tplot = self.T-self.kelvin; + timeplot = self.tday[n]/365 #will stamp time in years on your screen) + + icedepth = self.icethickness[n]*np.ones(np.size(lakedepthx)); + + #compute temperatures in water layer if no ice exist + +def printSomething(): + print ('This is working'); + +if __name__ == '__main__': + classFunc = thaLakeModel(); + classFunc.initialize(); + classFunc.updateModel(); + print ("Loop ended")