From 1a205829e5d92313350f7585a6b4499996181143 Mon Sep 17 00:00:00 2001 From: Marcus Gastaldello <116652649+MarcusGastaldello@users.noreply.github.com> Date: Thu, 27 Apr 2023 15:54:40 +0200 Subject: [PATCH 1/4] Update constants.py --- constants.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/constants.py b/constants.py index ee1f0fcc..09fbcc85 100644 --- a/constants.py +++ b/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) From bfa35ed076aeeea48aa046f6cabcc98d3a7514d7 Mon Sep 17 00:00:00 2001 From: Marcus Gastaldello <116652649+MarcusGastaldello@users.noreply.github.com> Date: Thu, 27 Apr 2023 16:09:42 +0200 Subject: [PATCH 2/4] Update heatEquation.py --- cosipy/modules/heatEquation.py | 56 ++++++++++++++++------------------ 1 file changed, 27 insertions(+), 29 deletions(-) diff --git a/cosipy/modules/heatEquation.py b/cosipy/modules/heatEquation.py index 0a2481ff..8297fb62 100644 --- a/cosipy/modules/heatEquation.py +++ b/cosipy/modules/heatEquation.py @@ -1,50 +1,48 @@ import numpy as np +from constants import basal_heat_flux from numba import njit @njit def solveHeatEquation(GRID, dt): """ Solves the heat equation on a non-uniform grid - dt :: integration time - + 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] """ - # number of layers - nl = GRID.get_number_layers() - - # 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()) - - # 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 - - # 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(((z[1:]+z[:-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 * (z[2:] + z[1:-1])) - \ + (0.5 * (K[:-2] + K[1:-1])) * (T[1:-1] - T[:-2]) / (0.5 * (z[:-2] + z[1:-1]))) / \ + (0.25 * z[:-2] + 0.5 * z[1:-1] + 0.25 * z[2:])) + # Update the temperature of the base sub-surface node: + 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 * (z[-2] + z[-1])))) / \ + (0.25 * z[-2] + 0.75 * z[-1])) T = Tnew.copy() # Write results to GRID From 3525021160c636b1ec0203fa7f34991df8cf6c13 Mon Sep 17 00:00:00 2001 From: Marcus Gastaldello <116652649+MarcusGastaldello@users.noreply.github.com> Date: Thu, 27 Apr 2023 16:10:19 +0200 Subject: [PATCH 3/4] Update heatEquation.py --- cosipy/modules/heatEquation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cosipy/modules/heatEquation.py b/cosipy/modules/heatEquation.py index 8297fb62..03c81fc2 100644 --- a/cosipy/modules/heatEquation.py +++ b/cosipy/modules/heatEquation.py @@ -30,7 +30,7 @@ def solveHeatEquation(GRID, dt): c_stab = 0.8 dt_stab = c_stab * min(((z[1:]+z[:-1])/2)**2/(K[1:]+K[:-1])) - # Numerically Solve the Fourier heat equation using a finite central difference scheme in matrix form: + # 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 From 38781883405d0dc2784a5afa660aa47f9f0c70be Mon Sep 17 00:00:00 2001 From: Marcus Gastaldello <116652649+MarcusGastaldello@users.noreply.github.com> Date: Fri, 2 Feb 2024 17:17:45 +0100 Subject: [PATCH 4/4] Add files via upload --- cosipy/modules/constants.py | 102 +++++++++++++++++++++++++++++++++ cosipy/modules/heatEquation.py | 49 ++++++++++------ 2 files changed, 134 insertions(+), 17 deletions(-) create mode 100644 cosipy/modules/constants.py diff --git a/cosipy/modules/constants.py b/cosipy/modules/constants.py new file mode 100644 index 00000000..5c9ac6c5 --- /dev/null +++ b/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 [-] diff --git a/cosipy/modules/heatEquation.py b/cosipy/modules/heatEquation.py index 03c81fc2..57937cc9 100644 --- a/cosipy/modules/heatEquation.py +++ b/cosipy/modules/heatEquation.py @@ -1,21 +1,30 @@ +from ast import Or import numpy as np -from constants import basal_heat_flux +from constants import basal_heat_flux, heat_equation_lower_boundary from numba import njit -@njit def solveHeatEquation(GRID, dt): + + 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))) + +@njit +def HeatEquation(GRID, dt): """ Solves the heat equation on a non-uniform grid Constants: - dt :: Integration time in a model time-step (hour) [s] - bhf :: Basal heat flux [W m-2] + 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] + 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] + T(z) :: Layer temperature (updated) [K] """ # Get sub-surface layer properties: @@ -28,7 +37,7 @@ def solveHeatEquation(GRID, dt): # Determine integration steps required in the solver to ensure numerical stability: stab_t = 0.0 c_stab = 0.8 - dt_stab = c_stab * min(((z[1:]+z[:-1])/2)**2/(K[1:]+K[:-1])) + 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: @@ -36,14 +45,20 @@ def solveHeatEquation(GRID, dt): stab_t = stab_t + dt_use # 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 * (z[2:] + z[1:-1])) - \ - (0.5 * (K[:-2] + K[1:-1])) * (T[1:-1] - T[:-2]) / (0.5 * (z[:-2] + z[1:-1]))) / \ - (0.25 * z[:-2] + 0.5 * z[1:-1] + 0.25 * z[2:])) - # Update the temperature of the base sub-surface node: - 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 * (z[-2] + z[-1])))) / \ - (0.25 * z[-2] + 0.75 * z[-1])) + 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) + +