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) 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 0a2481ff..57937cc9 100644 --- a/cosipy/modules/heatEquation.py +++ b/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) + +