Modification to the heat equation module

constants.py

@@ -72,6 +72,8 @@
 #albedo_mod_snow_aging = 6                      # effect of ageing on snow albedo [days] (Moelg et al. 2012, TC)
 #albedo_mod_snow_depth = 8                      # effect of snow depth on albedo [cm] (Moelg et al. 2012, TC)
 
+basal_heat_flux = 0.04                          # Basal/geothermal heat flux [Wm^-2]
+
 roughness_fresh_snow = 0.24                     # surface roughness length for fresh snow [mm] (Moelg et al. 2012, TC)
 roughness_ice = 1.7                             # surface roughness length for ice [mm] (Moelg et al. 2012, TC)
 roughness_firn = 4.0                            # surface roughness length for aged snow [mm] (Moelg et al. 2012, TC)

cosipy/modules/constants.py

@@ -0,0 +1,102 @@
+from config import WRF_X_CSPY
+"""
+    Declaration of constants
+    Do not modify unless you are absolutely sure what you are doing.
+"""
+
+' GENERAL INFORMATION ' 
+dt = 3600                                           # Time step in the input files [s]
+max_layers = 200                                    # Max. number of layers, just for the restart file
+z = 2.0                                             # Measurement height [m]
+
+' PARAMETERIZATIONS '
+stability_correction = 'Ri'                         # possibilities: 'Ri','MO'
+albedo_method = 'Oerlemans98'                       # possibilities: 'Oerlemans98'
+densification_method = 'Boone'                      # possibilities: 'Boone','empirical','constant' TODO: solve error Vionnet
+penetrating_method = 'Bintanja95'                   # possibilities: 'Bintanja95'
+roughness_method = 'Moelg12'                        # possibilities: 'Moelg12'
+saturation_water_vapour_method = 'Sonntag90'        # possibilities: 'Sonntag90'
+thermal_conductivity_method = 'bulk'		        # possibilities: 'bulk', 'empirical'
+sfc_temperature_method = 'SLSQP'                    # possibilities: 'L-BFGS-B', 'SLSQP'(faster), 'Newton' (Secant, fastest)'
+heat_equation_lower_boundary = 'basal_heat_flux'    # possibilities: 'basal_heat_flux','prescribed_temp'
+
+# WRF_X_CSPY: for efficiency and consistency
+if WRF_X_CSPY:
+    stability_correction = 'MO'
+    sfc_temperature_method = 'Newton'
+
+
+' INITIAL CONDITIONS '
+initial_snowheight_constant = 0.2                   # Initial snowheight
+initial_snow_layer_heights = 0.10                   # Initial thickness of snow layers
+initial_glacier_height = 40.0                       # Initial glacier height without snowlayers
+initial_glacier_layer_heights = 0.5                 # Initial thickness of glacier ice layers
+initial_top_density_snowpack = 300.0                # Top density for initial snowpack
+initial_bottom_density_snowpack = 600.0             # Bottom density for initial snowpack
+
+temperature_bottom = 270.16                         # Lower boundary condition for temperature profile (K) (initial temperature only for 'basal heat flux' method)
+const_init_temp = 0.1                               # constant for init temperature profile used in exponential function (exponential decay)
+
+zlt1 = 0.06					                        # First depth for temperature interpolation which is used for calculation of ground heat flux
+zlt2 = 0.1					                        # Second depth for temperature interpolation which is used for calculation of ground heat flux
+
+' MODEL CONSTANTS '
+center_snow_transfer_function = 1.0                 # center (50/50) when total precipitation is transferred to snow and rain
+spread_snow_transfer_function = 1                   # 1: +-2.5
+mult_factor_RRR = 1.0                               # multiplication factor for RRR
+
+minimum_snow_layer_height = 0.001                   # minimum layer height
+minimum_snowfall = 0.001                            # minimum snowfall per time step in m which is added as new snow
+
+
+' REMESHING OPTIONS'
+remesh_method = 'log_profile'                       # Remeshing (log_profile or adaptive_profile)
+first_layer_height = 0.01                           # The first layer will always have the defined height (m)
+layer_stretching = 1.20                             # Stretching factor used by the log_profile method (e.g. 1.1 mean the subsequent layer is 10% greater than the previous
+
+merge_max = 1                                       # How many mergings are allowed per time step
+density_threshold_merging = 5                       # If merging is true threshold for layer densities difference two layer try: 5-10 (kg m^-3)
+temperature_threshold_merging = 0.01                # If mering is true threshold for layer temperatures to merge  try: 0.05-0.1 (K)
+
+
+' PHYSICAL CONSTANTS '
+constant_density = 300.                             # constant density of freshly fallen snow [kg m-3], if densification_method is set to 'constant'
+
+albedo_fresh_snow = 0.85                            # albedo of fresh snow [-] (Moelg et al. 2012, TC)
+albedo_firn = 0.55                                  # albedo of firn [-] (Moelg et al. 2012, TC)
+albedo_ice = 0.3                                    # albedo of ice [-] (Moelg et al. 2012, TC)
+albedo_mod_snow_aging = 22                          # effect of ageing on snow albedo [days] (Oerlemans and Knap 1998, J. Glaciol.)
+albedo_mod_snow_depth = 3                           # effect of snow depth on albedo [cm] (Oerlemans and Knap 1998, J. Glaciol.)
+
+### For tropical glaciers or High Mountain Asia summer-accumulation glaciers (low latitude), the Moelg et al. 2012, TC should be tested for a possible better albedo fit 
+#albedo_mod_snow_aging = 6                          # effect of ageing on snow albedo [days] (Moelg et al. 2012, TC)
+#albedo_mod_snow_depth = 8                          # effect of snow depth on albedo [cm] (Moelg et al. 2012, TC)
+
+basal_heat_flux = 0.04                              # Basal/geothermal heat flux [W m^ (-2)]
+
+roughness_fresh_snow = 0.24                         # surface roughness length for fresh snow [mm] (Moelg et al. 2012, TC)
+roughness_ice = 1.7                                 # surface roughness length for ice [mm] (Moelg et al. 2012, TC)
+roughness_firn = 4.0                                # surface roughness length for aged snow [mm] (Moelg et al. 2012, TC)
+aging_factor_roughness = 0.0026                     # effect of ageing on roughness lenght (hours) 60 days from 0.24 to 4.0 => 0.0026
+
+snow_ice_threshold = 900.0                          # pore close of density [kg m^(-3)]
+
+lat_heat_melting = 3.34e5                           # latent heat for melting [J kg-1]
+lat_heat_vaporize = 2.5e6                           # latent heat for vaporization [J kg-1]
+lat_heat_sublimation = 2.834e6                      # latent heat for sublimation [J kg-1]
+
+spec_heat_air = 1004.67                             # specific heat of air [J kg-1 K-1]
+spec_heat_ice = 2050.00                             # specific heat of ice [J Kg-1 K-1]
+spec_heat_water = 4217.00                           # specific heat of water [J Kg-1 K-1]
+
+k_i = 2.22                                          # thermal conductivity ice [W m^-1 K^-1]
+k_w = 0.55                                          # thermal conductivity water [W m^-1 K^-1]
+k_a = 0.024                                         # thermal conductivity air [W m^-1 K^-1]
+
+water_density = 1000.0                              # density of water [kg m^(-3)]
+ice_density = 917.                                  # density of ice [kg m^(-3)]
+air_density = 1.1                                   # density of air [kg m^(-3)]
+
+sigma = 5.67e-8                                     # Stefan-Bolzmann constant [W m-2 K-4]
+zero_temperature = 273.16                           # Melting temperature [K]
+surface_emission_coeff = 0.99                       # surface emission coefficient [-]

cosipy/modules/heatEquation.py

@@ -1,51 +1,64 @@
+from ast import Or
 import numpy as np
+from constants import basal_heat_flux, heat_equation_lower_boundary
 from numba import njit
 
-@njit
 def solveHeatEquation(GRID, dt):
-    """ Solves the heat equation on a non-uniform grid
 
-    dt  ::  integration time
-
-    """
-    # number of layers
-    nl = GRID.get_number_layers()
+    heat_equation_lower_boundaries = ['basal_heat_flux','prescribed_temp']
+    if (heat_equation_lower_boundary == 'basal_heat_flux') or (heat_equation_lower_boundary == 'prescribed_temp'):
+        HeatEquation(GRID, dt)
+    else:
+        raise ValueError("heat_equation_lower_boundary = \"{:s}\" is not allowed, must be one of {:s}".format(heat_equation_lower_boundary, ", ".join(heat_equation_lower_boundaries)))
 
-    # Define index arrays 
-    k   = np.arange(1,nl-1) # center points
-    kl  = np.arange(2,nl)   # lower points
-    ku  = np.arange(0,nl-2) # upper points
-
-    # Get thermal diffusivity [m2 s-1]
-    K = np.asarray(GRID.get_thermal_diffusivity()) 
-
-    # Get snow layer heights    
-    hlayers = np.asarray(GRID.get_height())
+@njit
+def HeatEquation(GRID, dt):
+    """ Solves the heat equation on a non-uniform grid
 
-    # Get grid spacing
-    diff = ((hlayers[0:nl-1]/2.0)+(hlayers[1:nl]/2.0))
-    hk = diff[0:nl-2]  # between z-1 and z
-    hk1 = diff[1:nl-1] # between z and z+1
+    Constants:     
+                dt    ::  Integration time in a model time-step (hour) [s]
+                bhf   ::  Basal heat flux [W m-2]
+    Inputs:
+                h(z)  ::  Layer height [m]
+                K(z)  ::  Layer thermal diffusivity [m2 s-1]
+                k(z)  ::  Layer thermal conductivity [W m-1 K-1]
+                T(z)  ::  Layer temperature [K]    
+    Outputs:
+                T(z)  ::  Layer temperature (updated) [K]       
+    """
 
-    # Get temperature array from grid|
-    T = np.array(GRID.get_temperature())
+    # Get sub-surface layer properties:
+    h = np.asarray(GRID.get_height())
+    K = np.asarray(GRID.get_thermal_diffusivity())
+    T = np.asarray(GRID.get_temperature())
+    k = np.asarray(GRID.get_thermal_conductivity())
     Tnew = T.copy()
 
-    Kl = (K[1:nl-1]+K[2:nl])/2.0
-    Ku = (K[0:nl-2]+K[1:nl-1])/2.0
-
+    # Determine integration steps required in the solver to ensure numerical stability:
     stab_t = 0.0
     c_stab = 0.8
-    dt_stab  = c_stab * (min([min(diff[0:nl-2]**2/(2*Ku)),min(diff[1:nl-1]**2/(2*Kl))]))
+    dt_stab  = c_stab * min(((h[1:]+h[:-1])/2)**2/(K[1:]+K[:-1]))
 
+    # Numerically solve the Fourier heat equation using a finite central difference scheme in matrix form:
     while stab_t < dt:
-
         dt_use = np.minimum(dt_stab, dt-stab_t)
         stab_t = stab_t + dt_use
 
-        # Update the temperatures
-        Tnew[k] += ((Kl*dt_use*(T[kl]-T[k])/(hk1)) - (Ku*dt_use*(T[k]-T[ku])/(hk))) / (0.5*(hk+hk1))
+        # Update the temperatures of the intermediate sub-surface nodes:
+        Tnew[1:-1] = T[1:-1] + dt_use * (((0.5 * (K[2:]  + K[1:-1]))  * (T[2:]   - T[1:-1]) / (0.5 * (h[2:]  + h[1:-1])) - \
+                                        (0.5 * (K[:-2] + K[1:-1]))  * (T[1:-1] - T[:-2])  / (0.5 * (h[:-2] + h[1:-1]))) / \
+                                        (0.25 * h[:-2]  + 0.5 * h[1:-1] + 0.25 * h[2:]))
+
+        if heat_equation_lower_boundary == 'basal_heat_flux':
+
+            # Update the temperature of the base sub-surface node using the basal heat flux:
+            Tnew[-1]   = T[-1]   + dt_use * (((basal_heat_flux * K[-1] / k[-1]) - \
+                                            ((0.5 * (K[-2] + K[-1])) * (T[-1] - T[-2])  / (0.5 * (h[-2] + h[-1])))) / \
+                                            (0.25 * h[-2] + 0.75 * h[-1]))
+
         T = Tnew.copy()
 
     # Write results to GRID
     GRID.set_temperature(T)
+
+