updating to v1.1.11. Details in changelog.

CFM_main/CFMoutput_example/csv/example_csv.json

@@ -5,6 +5,7 @@
 "InputFileNameIso": "example_ISOTOPE.csv",
 "InputFileNamerho": "example_RHOS.csv",
 "InputFileNamemelt": "example_SMELT.csv",
+"InputFileNameStrain": "example_STRAIN.csv",
 "InputFileNameSublim": "example_SUBLIM.csv",
 "InputFileNameRain": "example_RAIN.csv",
 "resultsFolder": "CFMoutput_example/csv",
@@ -55,8 +56,10 @@
 "iso": ["d18O","dD"],
 "isoOptions":["d18O","dD","NoDiffusion"],
 "spacewriteint": 1,
-"strain": false,
-"du_dx": 1e-5,
+"horizontal_divergence": false,
+"strain_softening": false,
+"tuning_bias_correction": false,
+"residual_strain": 0.0002,
 "outputs": ["density", "depth", "temperature", "age", "DIP","LWC","meltoutputs","climate"],
 "outputs_options": ["density", "depth", "temperature", "age", "Dcon", "bdot_mean", "climate", "compaction", "grainsize", "temp_Hx", "isotopes", "BCO", "DIPc", "DIP", "LWC","gasses", "PLWC_mem", "viscosity", "runoff", "refrozen"],
 "resultsFileName": "CFMresults.hdf5",

CFM_main/CFMoutput_example/df/example_df.json

@@ -5,6 +5,7 @@
 "InputFileNameIso": "example_ISOTOPE.csv",
 "InputFileNamerho": "example_RHOS.csv",
 "InputFileNamemelt": "example_SMELT.csv",
+"InputFileNameStrain": "example_STRAIN.csv",
 "InputFileNameSublim": "example_SUBLIM.csv",
 "InputFileNameRain": "example_RAIN.csv",
 "resultsFolder": "CFMoutput_example/df",
@@ -55,8 +56,10 @@
 "iso": ["d18O","dD"],
 "isoOptions":["d18O","dD","NoDiffusion"],
 "spacewriteint": 1,
-"strain": false,
-"du_dx": 1e-5,
+"horizontal_divergence": false,
+"strain_softening": false,
+"tuning_bias_correction": false,
+"residual_strain": 0.0002,
 "outputs": ["density", "depth", "temperature", "age", "DIP","LWC","meltoutputs","climate"],
 "outputs_options": ["density", "depth", "temperature", "age", "Dcon", "bdot_mean", "climate", "compaction", "grainsize", "temp_Hx", "isotopes", "BCO", "DIPc", "DIP", "LWC","gasses", "PLWC_mem", "viscosity", "runoff", "refrozen"],
 "resultsFileName": "CFMresults.hdf5",

CFM_main/firn_density_nospin.py

@@ -445,13 +445,18 @@ def __init__(self, configName, climateTS = None, NewSpin = False):
                 self.iceout = self.c['iceout']
                 print('Ensure that your iceout value has units m ice eq. per year!')
             else:
-                self.iceout = np.mean(self.bdot+self.sublim) # this is the rate of ice flow advecting out of the column, units m I.E. per year.
-
+                if self.c['SUBLIM']:
+                    self.iceout = np.mean(self.bdot+self.sublim) # this is the rate of ice flow advecting out of the column, units m I.E. per year.
+                else:
+                    self.iceout = np.mean(self.bdot) # this is the rate of ice flow advecting out of the column, units m I.E. per year.
         except Exception:
             print('add field "manual_iceout" to .json file to set iceout value manually')
             self.iceout = np.mean(self.bdot) # this is the rate of ice flow advecting out of the column, units m I.E. per year.
 
-        self.w_firn         = np.mean(self.bdot+self.sublim) * RHO_I / self.rho 
+        if self.c['SUBLIM']:
+            self.w_firn         = np.mean(self.bdot+self.sublim) * RHO_I / self.rho 
+        else:
+            self.w_firn         = np.mean(self.bdot) * RHO_I / self.rho 
 
         if (np.any(self.bdotSec<0.0) and self.c['bdot_type']=='instant'):
             print('ERROR: bdot_type set to "instant" in .json and input')
@@ -946,13 +951,13 @@ def time_evolve(self):
                 if self.c['liquid'] == 'bucket':
 
                     if (self.snowmeltSec[iii]>0) or (np.any(self.LWC > 0.)) or (self.rainSec[iii] > 0.): #i.e. there is water
-                        self.rho, self.age, self.dz, self.Tz, self.r2, self.z, self.mass, self.dzn, self.LWC, meltgridtrack, self.refreeze, self.runoff = bucket(self,iii)
+                        self.rho, self.age, self.dz, self.Tz, self.r2, self.z, self.mass, self.dzn, self.LWC, meltgridtrack, self.refreeze, self.runoff, self.dh_melt = bucket(self,iii)
                         if self.doublegrid==True: # if we use doublegrid -> use the gridtrack corrected for melting
                             self.gridtrack = np.copy(meltgridtrack)
                         self.meltvol = self.snowmeltSec[iii]*S_PER_YEAR*0.917 #[m w.e.]
                     else: # Dry firn column and no input of meltwater                        
                         self.dzn     = self.dz[0:self.compboxes] # Not sure this is necessary
-                        self.refreeze, self.runoff, self.meltvol = 0.,0.,0.
+                        self.refreeze, self.runoff, self.meltvol, self.dh_melt = 0.,0.,0.,0.
                 ### end bucket ##################
 
                 elif self.c['liquid'] == 'darcy':
@@ -1032,6 +1037,7 @@ def time_evolve(self):
 
             else: # no melt, dz after compaction
                 self.dzn    = self.dz[0:self.compboxes]
+                self.dh_melt = 0
             ### end MELT #########
             ######################
 
@@ -1131,6 +1137,7 @@ def time_evolve(self):
             # MS 2/10/17: should double check that everything occurs in correct order in time step (e.g. adding new box on, calculating dz, etc.)               
                 self.age        = np.concatenate(([0], self.age[:-1])) + self.dt[iii]      
                 self.dzNew      = self.bdotSec[iii] * RHO_I / self.rhos0[iii] * S_PER_YEAR
+                self.dh_acc = self.dzNew
                 self.dz         = np.concatenate(([self.dzNew], self.dz[:-1]))
                 self.z          = self.dz.cumsum(axis = 0)
                 znew = np.copy(self.z) 
@@ -1157,6 +1164,7 @@ def time_evolve(self):
                 self.z          = self.dz.cumsum(axis=0)
                 self.z          = np.concatenate(([0],self.z[:-1]))
                 self.dzNew      = 0
+                self.dh_acc = 0
                 znew = np.copy(self.z)                             
                 self.compaction = (self.dz_old[0:self.compboxes]-self.dzn)
                 if not self.c['manualT']:
@@ -1167,9 +1175,10 @@ def time_evolve(self):
             if self.c['SUBLIM']:
                 if self.sublimSec[iii]<0:
                     self.mass_sum   = self.mass.cumsum(axis = 0)
-                    self.rho, self.age, self.dz, self.Tz, self.r2, self.z, self.mass, self.dzn, self.LWC, self.PLWC_mem, self.totwatersublim, sublgridtrack     = sublim(self,iii) # keeps track of sublimated water for mass conservation                    
+                    self.rho, self.age, self.dz, self.Tz, self.r2, self.z, self.mass, self.dzn, self.LWC, self.PLWC_mem, self.totwatersublim, sublgridtrack, dh_sub     = sublim(self,iii) # keeps track of sublimated water for mass conservation                    
                     self.compaction = (self.dz_old[0:self.compboxes]-self.dzn)
-                    self.dzNew      = 0
+                    # self.dzNew      = 0
+                    self.dh_acc = self.dh_acc + dh_sub #dh_sub is negative 
                     if self.doublegrid == True: # gridtrack corrected for sublimation
                         self.gridtrack = np.copy(sublgridtrack)
                     znew = np.copy(self.z)
@@ -1376,7 +1385,7 @@ def update_dH(self,iii):
         updates the surface elevation change
         '''
 
-        self.dH     = (self.sdz_new - self.sdz_old) + self.dzNew - (self.iceout*self.t[iii]) # iceout has units m ice/year, t is years per time step. 
+        self.dH     = (self.sdz_new - self.sdz_old) + self.dh_acc + self.dh_melt - (self.iceout*self.t[iii]) # iceout has units m ice/year, t is years per time step. 
         # self.dH2 = self.z[-1] - self.z_old[-1] #- (self.iceout*self.t) # alternative method. Should be the same?    
         self.dHAll.append(self.dH)
         self.dHtot  = np.sum(self.dHAll)
@@ -1386,7 +1395,7 @@ def update_dH(self,iii):
         ### domain and the 917 density horizon.
 
         iceout_corr = self.iceout*RHO_I/self.rho[-1]
-        self.dHcorr = (self.sdz_new - self.sdz_old) + self.dzNew - (iceout_corr*self.t[iii]) # iceout has units m ice/year, t is years per time step. 
+        self.dHcorr = (self.sdz_new - self.sdz_old) + self.dh_acc + self.dh_melt - (iceout_corr*self.t[iii]) # iceout has units m ice/year, t is years per time step. 
 
         self.dHAllcorr.append(self.dHcorr)
         self.dHtotcorr = np.sum(self.dHAllcorr)

CFM_main/main.py

@@ -22,7 +22,7 @@
 
 __author__ = "C. Max Stevens, Vincent Verjans, Brita Horlings, Annika Horlings, Jessica Lundin"
 __license__ = "MIT"
-__version__ = "1.1.10"
+__version__ = "1.1.11"
 __maintainer__ = "Max Stevens"
 __email__ = "[email protected]"
 __status__ = "Production"

CFM_main/melt.py

@@ -92,6 +92,12 @@ def bucket(self,iii):
     pm_lwc  = self.LWC[ind1]/self.dz[ind1]*pm_dz #LWC of the pm node
     pm_Tz   = T_MELT
 
+    dzo = self.dz.copy()
+    if ind1>0:
+        dh_melt = -1 * (np.sum(dzo[0:ind1]) + (dzo[ind1]-pm_dz))
+    else:
+        dh_melt = -1 * (dzo[ind1]-pm_dz)
+
     ### Liquid water input at the surface ###
     liq_in_mass = max(melt_mass + (np.sum(self.LWC[0:ind1+1]) - pm_lwc) * RHO_W_KGM, 0) #avoid negative lwcinput due to numerical round-off errors
     liq_in_vol  = liq_in_mass/RHO_W_KGM
@@ -421,7 +427,7 @@ def bucket(self,iii):
 
     self.rho[self.rho>RHO_I] = RHO_I
 
-    return self.rho, self.age, self.dz, self.Tz, self.r2, self.z, self.mass, self.dzn, self.LWC, meltgridtrack, refrozentot, runofftot
+    return self.rho, self.age, self.dz, self.Tz, self.r2, self.z, self.mass, self.dzn, self.LWC, meltgridtrack, refrozentot, runofftot, dh_melt
 
 #############
 

CFM_main/sublim.py

@@ -36,6 +36,7 @@ def sublim(self,iii):
         sys.exit()
 
     # ps is the partial sublimation (the model volume that has a portion sublimated away)   
+    dzo = self.dz.copy()
     if ind1 > 0: # if we sublimate the entire layers we retrieve the mass+lwc of these layers to the sublimated mass of the ps layer
         ps_sublim = sublim_mass - (np.cumsum(self.mass[ind1-1])+1000*np.cumsum(self.LWC[ind1-1])) 
     elif ind1 == 0: # if only the surface layer gets partially sublimated, all sublimation occurs in surface layer (logically)
@@ -46,6 +47,10 @@ def sublim(self,iii):
     ps_plwc = np.maximum(self.PLWC_mem[ind1] - ps_sublimlwc/1000, 0.) # assumes plwc is sublimated before mlwc
     ps_mass = np.maximum(1.0e-6,(self.mass[ind1] - ps_sublimmass)) # finally ice is sublimated if there is still mass to sublimate (<-> if ps_sublim<1000*self.LWC[ind1]), avoid rounding errors to cause negative mass
     ps_dz = ps_mass / self.rho[ind1] # remaining thickness [m]
+    if ind1>0:
+        dh_sub = -1 * (np.sum(dzo[0:ind1]) + (dzo[ind1]-ps_dz))
+    else:
+        dh_sub = -1 * (dzo[ind1]-ps_dz)
 
     ## Sublimated boxes are accomodated by just adding more (new) boxes at the bottom of the column
     ## Beware of this if you are not modeling to firn-ice transition depth.
@@ -89,6 +94,6 @@ def sublim(self,iii):
         print('setting to zero and continuing')
     self.LWC[self.LWC<0]=0.0
 
-    return self.rho, self.age, self.dz, self.Tz, self.r2, self.z, self.mass, self.dzn, self.LWC, self.PLWC_mem, self.totwatersublim,sublgridtrack
+    return self.rho, self.age, self.dz, self.Tz, self.r2, self.z, self.mass, self.dzn, self.LWC, self.PLWC_mem, self.totwatersublim, sublgridtrack, dh_sub
 
 

changelog.md

@@ -12,7 +12,7 @@ git push origin vX.Y.Z
 Then, on github do a release, which will trigger an updated DOI. 
 
 ## Current Version
-1.1.10
+1.1.11
 
 ## Full Documentation
 
@@ -31,6 +31,17 @@ https://communityfirnmodel.readthedocs.io/en/latest/
 	- Goujon physics work, but could possibly be implemented more elegantly (it would be nice to avoid globals)
 	- Not exactly in progress, but at some point adding a log file that gets saved in the results folder would be a good idea.
 
+## [1.1.11] 2022-12-13
+### Notes
+- This release fixes an issue in the estimated surface elevation change *dh*, which was not accounting properly for elevation change due to sublimation and melt processes (it was just considerine dh due to firn compaction and new snow acccumulation). The new code explicitly includes *dh_melt* and *dh_acc*, which are the elevation change (for that time step) due to melt and accumulation + sublimation. **The elevation change calculation should be considered to be in beta.** 
+
+### Changed
+- *firn_density_nospin.py, melt.py, sublim.py* firn_density_nospin's *update_dh* function is now:
+*dH = (sdz_new - sdz_old) + dh_acc + dh_melt - (iceout \* t[iii])*
+The sdz terms are the sum of the layer thicknesses before and after compaction (different is thus dh from firn compaction); dh_acc is the elevation change due to new snow accumulation minus the sublimated volume; dh_melt is the elevation decrease due to surface melt; iceout (rate of m ice e.q./year) times the time step size is the elevation change due to ice flow. Note that this assumes steady state, and that the ice flow is calculated using the spin up climate (iceout is set to be the mean ice-equivalent accumulation rate during the spin up), unless set explicitly in the config file. 
+
+>> *dh_melt* is calculated in *melt.py*, and *dh_acc* is calculated using dh_sub, which is now returned by *sublim.py*.
+
 ## [1.1.10]
 ### Notes
 - Version 1.1.10 is a minor update to fix an issue with the strain softening routine introduced in v1.1.9. (The code would throw an error if InputFileNamedudx was not in the .json).