From 5297ab4ae7f2ee94d463a6722dd85a8c06e45a7e Mon Sep 17 00:00:00 2001 From: Marcus Gastaldello <116652649+MarcusGastaldello@users.noreply.github.com> Date: Sat, 27 Jan 2024 20:39:57 +0100 Subject: [PATCH 1/2] Add files via upload Added Ligtenberg densification --- cosipy/modules/densification.py | 64 +++++++++++++++++++++++++++++++-- 1 file changed, 61 insertions(+), 3 deletions(-) diff --git a/cosipy/modules/densification.py b/cosipy/modules/densification.py index fa7e16a1..04109655 100644 --- a/cosipy/modules/densification.py +++ b/cosipy/modules/densification.py @@ -1,6 +1,5 @@ import numpy as np -from constants import densification_method, snow_ice_threshold, minimum_snow_layer_height, \ - zero_temperature, ice_density +from constants import * from numba import njit def densification(GRID,SLOPE,dt): @@ -10,11 +9,13 @@ def densification(GRID,SLOPE,dt): dt :: integration time """ - densification_allowed = ['Boone', 'Vionnet', 'empirical', 'constant'] + densification_allowed = ['Boone', 'Vionnet', 'empirical', 'constant','Ligtenberg11'] if densification_method == 'Boone': method_Boone(GRID,SLOPE,dt) elif densification_method == 'Vionnet': method_Vionnet(GRID,SLOPE,dt) + if densification_method == 'Ligtenberg11': + method_Ligtenberg(GRID,SLOPE,dt) elif densification_method == 'empirical': method_empirical(GRID,SLOPE,dt) elif densification_method == 'constant': @@ -162,6 +163,63 @@ def method_Vionnet(GRID,SLOPE,dt): GRID.get_node_porosity(idxNode)) print('Fraction > 1: %.5f' % (GRID.get_node_ice_fraction(idxNode)+GRID.get_node_liquid_water_content(idxNode)+GRID.get_node_porosity(idxNode))) +@njit +def method_Ligtenberg(GRID,SLOPE,dt): + """ Description: Densification based on in situ measurements of Antarctic snow compaction + after Arthern et al. 2010 (modified by Ligtenberg et al. 2011) + """ + + # Constants + g = 9.81 # gravitational acceleration [m s-2] + R = 8.314 # universal gas constant [J mol-1] + Ec = 60e3 # creep by lattice diffusion activation energy [J mol-1] + Eg = 42.4e3 # grain growth activation energy [J mol-1] + + # Extract variables: + rho = np.array(GRID.get_density()) + h = np.array(GRID.get_height()) + z = np.array(GRID.get_depth()) + T = np.array(GRID.get_temperature()) + theta_ice = np.array(GRID.get_ice_fraction()) + theta_water = np.array(GRID.get_liquid_water_content()) + + # Convert units: + b = accumulation # accumulation [mm a-1] + dt = dt / 31536000 # timestep as a fraction of a calendar year + + # Binary mask for snow/ice determination: + mask = np.where(rho < snow_ice_threshold,1,0) + + # Gravitational Constant: + C = np.where(rho < 550 , 0.07 * np.maximum(1.435 - 0.151 * np.log(b) , 0.25) , 0.03 * np.maximum(2.366 - 0.293 * np.log(b) , 0.25)) + + """ Note: Ligtenberg's method required average annual temperature of layers - a value difficult to calculate given layer movement and remeshing. + An approximation is made by assuming a thermal regime ressembling a trumpet curve and using linear interpolation""" + + # Obtain approximate temperatures from interpolation depths (domain base and depth of zero annual amplitude) + Tz1 = T[np.abs(z - z_zaa).argmin()] + Tz2 = T[-1] + + # Calculate linear temperature profile model constants + m = (Tz2 - Tz1) / (z[-1] - z_zaa) + q = Tz1 - (m * z_zaa) + + # Calculate approximate annual average temperature from linear gradient interpolation + T_AVG = z * m + q + + # Layer density change: + drho = mask * dt * C * b * g * (ice_density - rho) * np.exp((-Ec / (R * T)) + (Eg / (R * T_AVG))) + + # Calculate change in volumetric ice fraction & liquid water content: + dtheta_i = np.where(theta_water == 0, drho / ice_density , drho / 2 * ice_density) + dtheta_w = np.where(theta_water == 0, 0 , drho / 2 * water_density) + + # Set updated volumetric ice fraction & liquid water content: + GRID.set_ice_fraction(dtheta_i + theta_ice) + GRID.set_liquid_water_content(dtheta_w + theta_water) + + # Set updated layer height due to compaction: + GRID.set_height(np.maximum(minimum_snow_layer_height, h * (rho / (rho + drho)))) def method_empirical(GRID,SLOPE,dt): From 5de97b60984a7a0825024ca5a41948cebdb9fecc Mon Sep 17 00:00:00 2001 From: Marcus Gastaldello <116652649+MarcusGastaldello@users.noreply.github.com> Date: Sat, 27 Jan 2024 20:58:17 +0100 Subject: [PATCH 2/2] Update constants.py --- constants.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/constants.py b/constants.py index 45340ec7..b2bbb284 100644 --- a/constants.py +++ b/constants.py @@ -12,7 +12,7 @@ ' PARAMETERIZATIONS ' stability_correction = 'Ri' # possibilities: 'Ri','MO' albedo_method = 'Oerlemans98' # possibilities: 'Oerlemans98','Bougamont05' -densification_method = 'Boone' # possibilities: 'Boone','empirical','constant' TODO: solve error Vionnet +densification_method = 'Boone' # possibilities: 'Boone','Ligtenberg11','empirical','constant' TODO: solve error Vionnet penetrating_method = 'Bintanja95' # possibilities: 'Bintanja95' roughness_method = 'Moelg12' # possibilities: 'Moelg12' saturation_water_vapour_method = 'Sonntag90' # possibilities: 'Sonntag90' @@ -48,6 +48,8 @@ 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 +accumulation = 1000 # accumulation rate (mm year^-1) (for Ligtenberg 2011 densification) +z_zaa = 5 # estimated depth of zero annual accumulation (for Ligtenberg 2011 densification) ' REMESHING OPTIONS' remesh_method = 'log_profile' # Remeshing (log_profile or adaptive_profile)