critical_crack_length module#
Analytical parametrizations for the critical crack length
This module implements the parametrization by Richter et al (2019) for the critical crack length (flat terrain). It allows computation on CPU or GPU depending on user environment and install.
Other parametrizations are welcome to be added.
- xsnow.extensions.critical_crack_length.compute_rho_slab(density, thickness, xp)#
Compute the average slab density above each physical layer.
For each layer i, we define the “slab” as all layers strictly above i (shallower, closer to the snow surface). We return the mean density of that overlying slab. The topmost physical layer has no slab above, so it receives NaN. Padded layers above the real surface are also NaN.
Assumptions#
Axis -1 is the ‘layer’ dimension.
layer=0 is the deepest layer (near ground).
layer increases upward toward the snow surface.
After the true surface layer, remaining entries are padding with NaN in both density and thickness.
thickness is the physical thickness [m] of each layer. Zero/NaN thickness means “not a real layer” (padding).
- type density:
- param density:
Shape (…, L). Layer density [kg m^-3], ordered bottom (0) -> surface.
- type density:
xp.ndarray
- type thickness:
- param thickness:
Shape (…, L). Layer thickness [m], same ordering and shape.
- type thickness:
xp.ndarray
- type xp:
- param xp:
Array backend (e.g. numpy or cupy).
- type xp:
module
- returns:
rho_slab – Shape (…, L). For each layer i, rho_slab[i] is the mean density of all VALID layers above i. The surface layer itself (no overlying slab) is NaN. Padded layers beyond the surface remain NaN.
- rtype:
xp.ndarray
Notes
We mask out any padded layers (thickness <= 0 or NaN) so they don’t contaminate the cumulative sums.
We reverse along the layer axis, do cumsums from the top down, then reverse back. Finally, we shift so that layer i gets only the slab strictly above i.
- xsnow.extensions.critical_crack_length.compute_rc_flat_richter2019(rho_layer, grain_size_mm, shear_strength_kpa, rho_slab, xp)#
Compute critical crack length rc_flat [m] following Richter et al. (2019).
This returns the critical self-propagating crack length for a weak layer under a slab on flat terrain (i.e. slope effects are NOT included). Layers with no overlying slab (surface layer, padding) get NaN.
- Parameters:
rho_layer (xp.ndarray) – Layer density [kg m^-3], shape (…, L). This is the candidate failure layer density.
grain_size_mm (xp.ndarray) – Grain size [mm] for each layer, shape (…, L).
shear_strength_kpa (xp.ndarray) – Shear strength [kPa] for each layer, shape (…, L). (kPa = kN/m^2; will be converted to Pa.)
rho_slab (xp.ndarray) – Mean density [kg m^-3] of the slab ABOVE each layer (output of compute_rho_slab), shape (…, L). Should be NaN for the surface layer and padded layers.
xp (module) – Array backend (numpy or cupy).
- Returns:
rc_flat – Shape (…, L). Critical crack length [m] for each layer. NaN where no slab overlies the layer.
- Return type:
xp.ndarray
Notes
- This uses:
E’ = E0 * (rho_slab / RHO_ICE)**E_EXP / (1 - POISSON^2) D_sl/sigma_n = 1 / (g * rho_slab) term1 ~ sqrt( A * (…)**B ) term2 ~ sqrt( 2 * tau_p * E’ * D_sl/sigma_n ) rc_flat = term1 * term2
where tau_p is shear strength converted from kPa to Pa, and grain size is converted from mm to m.
No slope angle correction is applied here; this is for flat terrain (normal loading only), like the original xsnow rc_flat calculation.
- xsnow.extensions.critical_crack_length.compute_critical_crack_length(ds)#
Compute and attach the variable rc_flat (critical crack length on flat terrain, in meters) to this xsnowDataset.
This wraps the physics in compute_rho_slab() and compute_rc_flat_richter2019() and handles: - pulling data out of the xarray-backed dataset - computing layer thickness from cumulative height - dispatching to NumPy or CuPy depending on backend - wrapping the result back into xarray with metadata
Requirements#
- The dataset must contain:
density [kg m^-3]
grain_size [mm]
shear_strength [kPa]
height [m], cumulative height from ground upward
- We assume:
Dimension layer exists and indexes from ground (layer=0) upward.
height is monotonically increasing along layer.
Layers above the true snow surface are padded with NaN thickness.
- Slope angle / aspect:
Slope angle (inclination) and aspect (azimuth) are NOT used here. This is a flat-terrain formulation (same as in earlier xsnow versions).
- returns:
new_ds – A new xsnowDataset with a new data variable rc_flat [m].
- rtype:
xsnowDataset