earthsciencekit.glaciology

Domain hub · History · Examples · API reference

Inspectable glaciology models of glaciers, ice sheets, sea ice, rebound and ice cores; SI units unless stated otherwise.

earthsciencekit.glaciology.accumulation_area_ratio(accumulation_area, glacier_area)

Fraction of a glacier’s area above its equilibrium line.

A glacier in balance with its climate keeps a characteristic share of its area in the accumulation zone, about 0.6.

Parameters:

• accumulation_area (float or array_like) – m^2; constraint: between 0 and glacier_area.

• glacier_area (float or array_like) – m^2; constraint: positive.

Returns:

Model value [1]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Meier (1962); the steady-state ratio depends on the glacier’s hypsometry and is not universal. Equation: accumulation_area / glacier_area. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.dansgaard_delta(temperature, slope, offset)

Oxygen-isotope ratio of precipitation from the air temperature.

Air loses its heavy isotopes first as it cools and rains out, so the snow at colder sites is lighter.

Parameters:

• temperature (float or array_like) – K, mean annual surface air temperature; constraint: positive.

• slope (float or array_like) – permil/K, about 0.69; constraint: positive.

• offset (float or array_like) – permil, the delta at 0 degC, about -13.6; constraint: real.

Returns:

Model value [permil], delta-18O relative to SMOW. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Dansgaard (1964): a spatial relation for mid- and high-latitude coastal and polar stations; the temporal slope at one site can differ. Equation: slope * (temperature - 273.15) + offset. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.dansgaard_johnsen_age(depth, thickness, kink_height, accumulation)

Age of ice at depth in the Dansgaard-Johnsen flow model.

The first time scale for a deep ice core, Camp Century, integrated the thinning of annual layers down to the kink and below it.

Parameters:

• depth (float or array_like) – m of ice below the surface; constraint: between 0 and thickness, exclusive of thickness.

• thickness (float or array_like) – m of ice; constraint: positive.

• kink_height (float or array_like) – m above the bed; constraint: between 0 (exclusive) and thickness.

• accumulation (float or array_like) – m/s of ice; constraint: positive.

Returns:

Model value [s]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Dansgaard and Johnsen (1969), the integral of the reciprocal of dansgaard_johnsen_downward_velocity; it tends to nye_age as the kink height goes to zero. Below, z = thickness - depth, h = kink_height, a = accumulation and s = 2 * thickness - h. Equation: s / (2 * a) * np.log(s / (2 * np.maximum(z, h) - h)) + h * s / a * (1 / np.minimum(z, h) - 1 / h). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.dansgaard_johnsen_downward_velocity(depth, thickness, kink_height, accumulation)

Vertical ice velocity in the Dansgaard-Johnsen flow model.

Ice frozen to its bed shears only near the bottom, so the vertical velocity falls linearly down to a kink and then quadratically to zero at the bed.

Parameters:

• depth (float or array_like) – m of ice below the surface; constraint: between 0 and thickness.

• thickness (float or array_like) – m of ice; constraint: positive.

• kink_height (float or array_like) – m above the bed; constraint: between 0 (exclusive) and thickness.

• accumulation (float or array_like) – m/s of ice; constraint: positive.

Returns:

Model value [m/s], downward; a year’s layer is this speed times one year thick. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Dansgaard and Johnsen (1969): steady state, constant horizontal strain rate above the kink and a linear decrease below it, no basal melting. Below, z = thickness - depth and h = kink_height. Equation: np.where(z >= h, accumulation * (2 * z - h) / (2 * thickness - h), accumulation * z**2 / (h * (2 * thickness - h))). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.dansgaard_temperature(delta, slope, offset)

Air temperature recorded by the oxygen-isotope ratio of snow.

Inverting Dansgaard’s relation turns the isotope profile of an ice core into a temperature history.

Parameters:

• delta (float or array_like) – permil, delta-18O relative to SMOW; constraint: real.

• slope (float or array_like) – permil/K; constraint: positive.

• offset (float or array_like) – permil; constraint: real.

Returns:

Model value [K]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Dansgaard (1964), inverted; the spatial slope is an approximation of the temporal one. Equation: 273.15 + (delta - offset) / slope. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.driving_stress(thickness, surface_slope, density, gravity)

Gravitational driving stress of a glacier or ice sheet.

The weight of the ice resolved down the surface slope is what the bed and the ice must resist.

Parameters:

• thickness (float or array_like) – m; constraint: nonnegative.

• surface_slope (float or array_like) – 1, tangent of the surface slope; constraint: nonnegative.

• density (float or array_like) – kg/m^3; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [Pa], typically 50-150 kPa under glaciers. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Shallow-ice balance: horizontal scales much larger than the thickness, small slopes, longitudinal stresses neglected. Equation: density * gravity * thickness * surface_slope. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.effective_viscosity(stress, rate_factor, exponent)

Effective viscosity of ice under Glen’s law.

Ice under a large stress behaves as a much less viscous fluid than ice under a small one.

Parameters:

• stress (float or array_like) – Pa, effective stress; constraint: positive.

• rate_factor (float or array_like) – Pa^-n/s; constraint: positive.

• exponent (float or array_like) – 1, Glen’s n; constraint: at least 1.

Returns:

Model value [Pa s]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Glen (1955), written as a stress-dependent Newtonian viscosity; with n = 1 it is the viscosity 1 / (2 A). Equation: 1 / (2 * rate_factor * stress ** (exponent - 1)). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.ela_shift(temperature_change, lapse_rate)

Shift of the equilibrium-line altitude with a change of climate.

A colder climate brings the snow line down the mountains, and the glaciers grow to meet it.

Parameters:

• temperature_change (float or array_like) – K, change of summer air temperature at fixed altitude; constraint: real.

• lapse_rate (float or array_like) – K/m, decrease of temperature with height; constraint: positive.

Returns:

Model value [m], positive upward, so cooling lowers the line. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Temperature alone controls ablation and the lapse rate is uniform; changes of snowfall also move the line. Equation: temperature_change / lapse_rate. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.flotation_thickness(bed_depth, ice_density, water_density)

Thickness at which ice resting on a submarine bed begins to float.

The grounding line of a marine ice sheet sits where the ice is just thin enough to float, so its thickness is set by the water depth alone.

Parameters:

• bed_depth (float or array_like) – m, below sea level; constraint: nonnegative.

• ice_density (float or array_like) – kg/m^3; constraint: positive.

• water_density (float or array_like) – kg/m^3; constraint: positive.

Returns:

Model value [m]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Archimedes’ principle for an ice column of uniform density. Equation: water_density * bed_depth / ice_density. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.gas_age_difference(close_off_depth, mean_firn_density, ice_density, accumulation)

Age difference between the ice and the air trapped in it.

Air stays connected to the atmosphere through the porous firn until the bubbles close, so the gas in a core is younger than the ice around it.

Parameters:

• close_off_depth (float or array_like) – m, depth of bubble close-off; constraint: nonnegative.

• mean_firn_density (float or array_like) – kg/m^3, mean density of the firn above it; constraint: between 0 (exclusive) and ice_density.

• ice_density (float or array_like) – kg/m^3; constraint: positive.

• accumulation (float or array_like) – m/s of ice; constraint: positive.

Returns:

Model value [s]; about 4 kyr at present-day Vostok. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Schwander and Stauffer (1984), as used for Vostok by Barnola et al. (1987): steady firn densification with no thinning above close-off; diffusion through the open pores, which makes the gas slightly older, is neglected. Equation: close_off_depth * mean_firn_density / (ice_density * accumulation). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.glacier_response_time(thickness, terminus_ablation)

Volume response time of a glacier to a change of climate.

A glacier adjusts to a new climate in about the time its terminus ablation takes to remove its own thickness.

Parameters:

• thickness (float or array_like) – m, a characteristic thickness; constraint: positive.

• terminus_ablation (float or array_like) – m/s of ice, ablation rate at the terminus as a positive magnitude; constraint: positive.

Returns:

Model value [s], an e-folding time. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Johannesson, Raymond and Waddington (1989): a linear response of volume; decades for valley glaciers, millennia for ice sheets. Equation: thickness / terminus_ablation. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.glen_rate_factor(temperature, reference_factor, activation_energy, reference_temperature)

Arrhenius temperature dependence of Glen’s rate factor.

Ice near its melting point is softer by more than an order of magnitude than ice at -30 degC.

Parameters:

• temperature (float or array_like) – K, ice temperature (relative to the pressure-melting point); constraint: positive.

• reference_factor (float or array_like) – Pa^-n/s, rate factor at reference_temperature; constraint: positive.

• activation_energy (float or array_like) – J/mol, creep activation energy; constraint: nonnegative.

• reference_temperature (float or array_like) – K; constraint: positive.

Returns:

Model value [Pa^-n/s]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

An Arrhenius law with one activation energy, about 60 kJ/mol below -10 degC and 115 kJ/mol above (Cuffey and Paterson 2010), for temperate-ice water contents neglected. Equation: reference_factor * np.exp(-activation_energy / R_GAS * (1 / temperature - 1 / reference_temperature)). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.glen_strain_rate(stress, rate_factor, exponent)

Glen’s flow law: strain rate of polycrystalline ice.

Glen’s creep tests showed that ice deforms at a rate proportional to roughly the cube of the stress, so it is far from a Newtonian fluid.

Parameters:

• stress (float or array_like) – Pa, effective (or shear) stress; constraint: nonnegative.

• rate_factor (float or array_like) – Pa^-n/s, the softness A; constraint: positive.

• exponent (float or array_like) – 1, Glen’s n, about 3; constraint: at least 1.

Returns:

Model value [1/s], effective strain rate. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Glen (1955): secondary creep of isotropic ice at stresses of about 0.05-1 MPa. Equation: rate_factor * stress**exponent. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.glen_stress(strain_rate, rate_factor, exponent)

Stress needed to deform ice at a given strain rate.

Inverting Glen’s law turns a measured stretching rate into the stress that drives it.

Parameters:

• strain_rate (float or array_like) – 1/s, effective strain rate; constraint: nonnegative.

• rate_factor (float or array_like) – Pa^-n/s; constraint: positive.

• exponent (float or array_like) – 1, Glen’s n; constraint: at least 1.

Returns:

Model value [Pa]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Glen (1955), inverted; secondary creep of isotropic ice. Equation: (strain_rate / rate_factor) ** (1 / exponent). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.glen_surface_velocity(thickness, surface_slope, rate_factor, exponent, density, gravity)

Surface velocity of a glacier deforming by Glen’s law.

The surface speed rises with the fourth power of thickness and the cube of slope, which is why thick glaciers flow so much faster than thin ones.

Parameters:

• thickness (float or array_like) – m; constraint: positive.

• surface_slope (float or array_like) – 1, tangent of the surface slope; constraint: nonnegative.

• rate_factor (float or array_like) – Pa^-n/s; constraint: positive.

• exponent (float or array_like) – 1, Glen’s n; constraint: at least 1.

• density (float or array_like) – kg/m^3; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [m/s]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Nye (1952) laminar flow with Glen (1955) rheology, no sliding. Equation: 2 * rate_factor / (exponent + 1) * (density * gravity * surface_slope) ** exponent * thickness ** (exponent + 1). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.glen_velocity_profile(height, thickness, surface_slope, rate_factor, exponent, density, gravity)

Deformation velocity through a glacier obeying Glen’s law.

Most of the shearing in a frozen-bed glacier happens near the bed, so the velocity is nearly uniform in its upper half.

Parameters:

• height (float or array_like) – m, above the bed; constraint: between 0 and thickness.

• thickness (float or array_like) – m; constraint: positive.

• surface_slope (float or array_like) – 1, tangent of the surface slope; constraint: nonnegative.

• rate_factor (float or array_like) – Pa^-n/s; constraint: positive.

• exponent (float or array_like) – 1, Glen’s n; constraint: at least 1.

• density (float or array_like) – kg/m^3; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [m/s], downslope, zero at a frozen bed. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Nye (1952) laminar flow with Glen (1955) rheology: a wide slab, uniform rate factor and no sliding. Equation: 2 * rate_factor / (exponent + 1) * (density * gravity * surface_slope) ** exponent * (thickness ** (exponent + 1) - (thickness - height) ** (exponent + 1)). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.haskell_viscosity(relaxation_time, wavelength, density, gravity)

Mantle viscosity inferred from a rebound relaxation time.

Haskell turned the uplift of Fennoscandia into the first estimate of the viscosity of the mantle, about 1e21 Pa s.

Parameters:

• relaxation_time (float or array_like) – s; constraint: positive.

• wavelength (float or array_like) – m, of the depression; constraint: positive.

• density (float or array_like) – kg/m^3, of the mantle; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [Pa s]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Haskell (1935), the inverse of viscous_relaxation_time for a uniform half-space. Equation: density * gravity * wavelength * relaxation_time / (4 * np.pi). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.ice_shelf_strain_rate(thickness, rate_factor, exponent, ice_density, water_density, gravity)

Spreading rate of a freely floating ice shelf.

A floating shelf spreads under the difference between its weight and the water pressure, faster the thicker it is.

Parameters:

• thickness (float or array_like) – m; constraint: nonnegative.

• rate_factor (float or array_like) – Pa^-n/s; constraint: positive.

• exponent (float or array_like) – 1, Glen’s n; constraint: at least 1.

• ice_density (float or array_like) – kg/m^3; constraint: less than water_density.

• water_density (float or array_like) – kg/m^3; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [1/s], along-flow stretching. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Weertman (1957) for an unconfined shelf spreading in one direction, without buttressing or side drag. Equation: rate_factor * (ice_density * gravity * thickness * (1 - ice_density / water_density) / 4) ** exponent. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.isostatic_depression(ice_thickness, ice_density, mantle_density)

Equilibrium depression of the land under an ice sheet.

The mantle under a wide ice sheet gives way until the displaced rock weighs as much as the ice, about a third of its thickness.

Parameters:

• ice_thickness (float or array_like) – m; constraint: nonnegative.

• ice_density (float or array_like) – kg/m^3; constraint: positive.

• mantle_density (float or array_like) – kg/m^3; constraint: positive.

Returns:

Model value [m], positive downward. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Jamieson (1865): local isostasy for a load much wider than the flexural wavelength of the lithosphere. Equation: ice_thickness * ice_density / mantle_density. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.linear_mass_balance(elevation, equilibrium_line_altitude, balance_gradient)

Surface mass balance rising linearly with altitude.

A glacier gains mass above its equilibrium line and loses it below, at a rate that grows steadily with the distance from the line.

Parameters:

• elevation (float or array_like) – m; constraint: real.

• equilibrium_line_altitude (float or array_like) – m; constraint: real.

• balance_gradient (float or array_like) – 1/s, metres of ice per second per metre of altitude; constraint: nonnegative.

Returns:

Model value [m/s] of ice, positive for net accumulation. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Meier (1962): the balance gradient is about 0.003-0.01 m of water per year per metre on mid-latitude glaciers; real profiles flatten at high altitude. Equation: balance_gradient * (elevation - equilibrium_line_altitude). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.nye_age(depth, thickness, accumulation)

Age of ice at depth when layers thin uniformly.

In a steady ice sheet each layer thins in proportion to its height above the bed, so age grows logarithmically with depth.

Parameters:

• depth (float or array_like) – m, of ice equivalent below the surface; constraint: between 0 and thickness, exclusive of thickness.

• thickness (float or array_like) – m of ice; constraint: positive.

• accumulation (float or array_like) – m/s of ice; constraint: positive.

Returns:

Model value [s]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Nye (1963): a steady divide with uniform vertical strain rate, no basal melting or sliding; the age diverges at the bed. Equation: thickness / accumulation * np.log(thickness / (thickness - depth)). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.nye_crevasse_depth(strain_rate, rate_factor, exponent, density, gravity)

Depth of a dry crevasse in stretching ice.

A crevasse can only stay open where the tensile stress from stretching exceeds the closing pressure of the ice, which caps dry crevasses at a few tens of metres.

Parameters:

• strain_rate (float or array_like) – 1/s, longitudinal stretching rate; constraint: nonnegative.

• rate_factor (float or array_like) – Pa^-n/s; constraint: positive.

• exponent (float or array_like) – 1, Glen’s n; constraint: at least 1.

• density (float or array_like) – kg/m^3; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [m]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Nye (1955): a field of closely spaced crevasses in uniform ice of density independent of depth; fracture toughness is neglected. Equation: 2 * (strain_rate / rate_factor) ** (1 / exponent) / (density * gravity). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.nye_plastic_profile(distance, yield_stress, density, gravity)

Surface height of a perfectly plastic ice sheet on a flat bed.

If ice yields at a fixed stress, the ice sheet takes a parabolic profile whose height depends only on the distance from its margin.

Parameters:

• distance (float or array_like) – m, inland from the margin; constraint: nonnegative.

• yield_stress (float or array_like) – Pa, about 100 kPa; constraint: positive.

• density (float or array_like) – kg/m^3; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [m]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Nye (1951): a horizontal bed, basal shear stress everywhere equal to the yield stress and no isostatic depression. Equation: np.sqrt(2 * yield_stress * distance / (density * gravity)). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.rebound_uplift_rate(time, initial_depression, relaxation_time)

Rate of postglacial uplift.

The land rises fastest just after deglaciation; Fennoscandia still rises by up to 9 mm a year.

Parameters:

• time (float or array_like) – s, since the ice was removed; constraint: nonnegative.

• initial_depression (float or array_like) – m, at removal; constraint: nonnegative.

• relaxation_time (float or array_like) – s; constraint: positive.

Returns:

Model value [m/s], upward. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Haskell (1935), the time derivative of remaining_depression. Equation: initial_depression / relaxation_time * np.exp(-time / relaxation_time). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.remaining_depression(time, initial_depression, relaxation_time)

Depression still to be recovered after the ice has gone.

After the ice sheet melts the land rises exponentially, and the uplift still to come falls by e every relaxation time.

Parameters:

• time (float or array_like) – s, since the ice was removed; constraint: nonnegative.

• initial_depression (float or array_like) – m, at removal; constraint: nonnegative.

• relaxation_time (float or array_like) – s; constraint: positive.

Returns:

Model value [m], positive downward. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Haskell (1935): one wavelength on a viscous half-space, with the load removed at once. Equation: initial_depression * np.exp(-time / relaxation_time). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.schoof_grounding_line_flux(grounding_thickness, rate_factor, exponent, friction_coefficient, friction_exponent, ice_density, water_density, gravity)

Schoof’s ice flux across the grounding line of a marine ice sheet.

The flux out of a marine ice sheet is set by the ice thickness at its grounding line, and rises so steeply with it that a grounding line on a bed deepening inland cannot be steady.

Parameters:

• grounding_thickness (float or array_like) – m, the flotation thickness at the grounding line; constraint: nonnegative.

• rate_factor (float or array_like) – Pa^-n/s; constraint: positive.

• exponent (float or array_like) – 1, Glen’s n; constraint: at least 1.

• friction_coefficient (float or array_like) – Pa (s/m)^m, the C of a basal stress C u^m; constraint: positive.

• friction_exponent (float or array_like) – 1, the m of a basal stress C u^m, such as 1/3; constraint: positive.

• ice_density (float or array_like) – kg/m^3; constraint: less than water_density.

• water_density (float or array_like) – kg/m^3; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [m^2/s], per unit width. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Schoof (2007): a flowline with a power-law sliding bed, an unbuttressed shelf and a boundary layer at the grounding line; steady states are stable only where the bed rises inland. Below, A = rate_factor, n = exponent, C = friction_coefficient, m = friction_exponent, rho_i = ice_density, rho_w = water_density, h = grounding_thickness. Equation: (A * (rho_i * gravity)**(n + 1) * (1 - rho_i / rho_w)**n / (4**n * C))**(1 / (m + 1)) * h**((m + n + 3) / (m + 1)). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.sea_ice_freeboard(thickness, ice_density, water_density)

Height of floating snow-free ice above the waterline.

Satellite altimeters measure the freeboard, and Archimedes’ principle turns it into a thickness about ten times larger.

Parameters:

• thickness (float or array_like) – m; constraint: nonnegative.

• ice_density (float or array_like) – kg/m^3; constraint: less than water_density.

• water_density (float or array_like) – kg/m^3; constraint: positive.

Returns:

Model value [m]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Hydrostatic equilibrium of a uniform floe without snow load. Equation: thickness * (1 - ice_density / water_density). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.shallow_ice_diffusivity(thickness, surface_slope, rate_factor, exponent, density, gravity)

Nonlinear diffusivity of the shallow-ice equation.

The shallow-ice approximation makes the ice thickness obey a diffusion equation whose diffusivity grows steeply with thickness and slope.

Parameters:

• thickness (float or array_like) – m; constraint: nonnegative.

• surface_slope (float or array_like) – 1, magnitude of the surface gradient; constraint: nonnegative.

• rate_factor (float or array_like) – Pa^-n/s; constraint: positive.

• exponent (float or array_like) – 1, Glen’s n; constraint: at least 1.

• density (float or array_like) – kg/m^3; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [m^2/s], so that the flux is minus the diffusivity times the surface gradient. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Mahaffy (1976): the shallow-ice approximation for Glen ice frozen to its bed. Below, A = rate_factor, n = exponent. Equation: 2 * A / (n + 2) * (density * gravity) ** n * thickness ** (n + 2) * surface_slope ** (n - 1). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.shallow_ice_flux(thickness, surface_slope, rate_factor, exponent, density, gravity)

Ice flux per unit width in the shallow-ice approximation.

In a thin ice sheet the flux at each point depends only on the local thickness and surface slope.

Parameters:

• thickness (float or array_like) – m; constraint: nonnegative.

• surface_slope (float or array_like) – 1, magnitude of the surface gradient; constraint: nonnegative.

• rate_factor (float or array_like) – Pa^-n/s; constraint: positive.

• exponent (float or array_like) – 1, Glen’s n; constraint: at least 1.

• density (float or array_like) – kg/m^3; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [m^2/s], down the surface gradient. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Hutter (1983) and Mahaffy (1976): Glen ice frozen to its bed, aspect ratio small, uniform rate factor. Equation: 2 * rate_factor / (exponent + 2) * (density * gravity * surface_slope) ** exponent * thickness ** (exponent + 2). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.simulate_flowline_glacier(bed_elevation, spacing, equilibrium_line_altitude, balance_gradient, time, *, initial_thickness=None, rate_factor=2.4e-24, exponent=3.0, density=917.0, gravity=9.80665)

Shallow-ice flowline glacier under a linear mass balance.

Parameters:

• bed_elevation (array_like, shape (m,), m >= 3) – Bed elevation at equally spaced nodes along the flowline, head first [m].

• spacing (float) – Positive scalar node spacing [m].

• equilibrium_line_altitude (float) – Scalar equilibrium-line altitude held from the first sample on [m].

• balance_gradient (float) – Positive scalar mass-balance gradient [1/s].

• time (array_like) – Nonnegative strictly increasing output times [s]; the initial state belongs to the first sample.

• initial_thickness (array_like, shape (m,), optional) – Nonnegative ice thickness at the first sample [m]; ice-free by default. Pass a steady state for one climate with the equilibrium line of another to see the response to a step.

• rate_factor (float, optional) – Positive scalar Glen rate factor [Pa^-n/s], 2.4e-24 for ice at 0 degC (Cuffey and Paterson 2010).

• exponent (float, optional) – Scalar Glen exponent of at least 1 [1].

• density (float, optional) – Positive scalar ice density [kg/m^3] and gravity [m/s^2].

• gravity (float, optional) – Positive scalar ice density [kg/m^3] and gravity [m/s^2].

Returns:

Glacier length [m] (the spacing times the number of nodes with more than 1 m of ice), cross-sectional area per unit width [m^2], and the thickness at every node [m], labeled thickness_<node>.

Return type:

TimeSeries

Notes

Solves dH/dt = b(b_bed + H) - dq/dx with the shallow-ice flux q = -D ds/dx of shallow_ice_diffusivity() and the balance of linear_mass_balance(), by the method of lines on a staggered grid in Mahaffy’s (1976) form, with integrators.integrate_ode (BDF) in time. The ends carry no flux, the thickness is kept nonnegative (ablation tapers to zero over the last metre of ice, and a solver overshoot below zero relaxes back over about a year, so the right-hand side stays continuous), and there is no sliding, calving, width variation or isostasy.

earthsciencekit.glaciology.stefan_growth_rate(thickness, temperature_difference, conductivity, density, latent_heat)

Basal growth rate of ice in Stefan’s problem.

Thin ice grows fast and thick ice slowly, because the conductive heat flux falls with thickness.

Parameters:

• thickness (float or array_like) – m; constraint: positive.

• temperature_difference (float or array_like) – K, freezing temperature minus surface temperature; constraint: nonnegative.

• conductivity (float or array_like) – W/(m K); constraint: positive.

• density (float or array_like) – kg/m^3; constraint: positive.

• latent_heat (float or array_like) – J/kg; constraint: positive.

Returns:

Model value [m/s]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Stefan (1891), quasi-steady conduction with no ocean heat flux. Equation: conductivity * temperature_difference / (density * latent_heat * thickness). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.stefan_ice_thickness(freezing_degree_time, conductivity, density, latent_heat)

Stefan’s thickness of ice grown by conduction through itself.

The latent heat released at the base of the ice must be conducted up through ice already formed, so the thickness grows only as the square root of the accumulated cold.

Parameters:

• freezing_degree_time (float or array_like) – K s, time integral of the freezing temperature minus the surface temperature; constraint: nonnegative.

• conductivity (float or array_like) – W/(m K), about 2.2 for ice; constraint: positive.

• density (float or array_like) – kg/m^3; constraint: positive.

• latent_heat (float or array_like) – J/kg; constraint: positive.

Returns:

Model value [m]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Stefan (1891): ice free of snow, with its surface at the air temperature, a linear temperature profile, no ocean heat flux and no initial ice; it overestimates growth under snow. Equation: np.sqrt(2 * conductivity * freezing_degree_time / (density * latent_heat)). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.vialov_dome_height(half_width, accumulation, rate_factor, exponent, density, gravity)

Central height of a steady shallow-ice sheet.

An ice sheet in balance with a uniform snowfall has a dome height that grows only as the square root of its width.

Parameters:

• half_width (float or array_like) – m, divide to margin; constraint: positive.

• accumulation (float or array_like) – m/s of ice; constraint: positive.

• rate_factor (float or array_like) – Pa^-n/s; constraint: positive.

• exponent (float or array_like) – 1, Glen’s n; constraint: at least 1.

• density (float or array_like) – kg/m^3; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [m]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Vialov (1958): a flat bed, uniform accumulation up to a fixed margin, shallow-ice flow without sliding. Below, A = rate_factor, n = exponent. Equation: (2 * ((n + 2) * accumulation / (2 * A * (density * gravity) ** n)) ** (1 / n)) ** (n / (2 * n + 2)) * np.sqrt(half_width). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.vialov_profile(distance, half_width, dome_height, exponent)

Surface profile of a steady shallow-ice sheet.

Glen’s law flattens the dome and steepens the margin of a steady ice sheet compared with Nye’s parabola.

Parameters:

• distance (float or array_like) – m, from the divide; constraint: between 0 and half_width.

• half_width (float or array_like) – m; constraint: positive.

• dome_height (float or array_like) – m; constraint: nonnegative.

• exponent (float or array_like) – 1, Glen’s n; constraint: at least 1.

Returns:

Model value [m]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Vialov (1958), see vialov_dome_height; the shallow-ice approximation fails within a few thicknesses of the divide and the margin. Equation: dome_height * (1 - (distance / half_width) ** ((exponent + 1) / exponent)) ** (exponent / (2 * exponent + 2)). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.viscous_relaxation_time(wavelength, viscosity, density, gravity)

Relaxation time of a surface depression on a viscous half-space.

A broad depression in a viscous mantle recovers faster than a narrow one, so the timing of rebound measures the viscosity.

Parameters:

• wavelength (float or array_like) – m, of the depression; constraint: positive.

• viscosity (float or array_like) – Pa s, of the mantle; constraint: positive.

• density (float or array_like) – kg/m^3, of the mantle; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [s], an e-folding time. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Haskell (1935): a uniform Newtonian half-space without an elastic lithosphere, which slows the shortest wavelengths. Equation: 4 * np.pi * viscosity / (density * gravity * wavelength). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.viscous_slab_velocity(height, thickness, surface_slope, viscosity, density, gravity)

Velocity in a Newtonian viscous slab flowing down a slope.

Forbes argued that a glacier flows like a very viscous fluid, fastest at the surface and slowest against its bed.

Parameters:

• height (float or array_like) – m, above the bed; constraint: between 0 and thickness.

• thickness (float or array_like) – m; constraint: positive.

• surface_slope (float or array_like) – 1, tangent of the surface slope; constraint: nonnegative.

• viscosity (float or array_like) – Pa s; constraint: positive.

• density (float or array_like) – kg/m^3; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [m/s], downslope, zero at a frozen bed. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Forbes (1843): a wide parallel-sided slab with a linear viscosity, no sliding and small slope. Equation: density * gravity * surface_slope * (2 * thickness * height - height**2) / (2 * viscosity). Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.water_filled_crevasse_depth(strain_rate, rate_factor, exponent, water_depth, ice_density, water_density, gravity)

Depth of a crevasse partly filled with meltwater.

Water is denser than ice, so a water-filled crevasse is pushed deeper than a dry one and can cut through a whole ice shelf.

Parameters:

• strain_rate (float or array_like) – 1/s, longitudinal stretching rate; constraint: nonnegative.

• rate_factor (float or array_like) – Pa^-n/s; constraint: positive.

• exponent (float or array_like) – 1, Glen’s n; constraint: at least 1.

• water_depth (float or array_like) – m, height of water standing in the crevasse; constraint: nonnegative.

• ice_density (float or array_like) – kg/m^3; constraint: positive.

• water_density (float or array_like) – kg/m^3; constraint: positive.

• gravity (float or array_like) – m/s^2; constraint: positive.

Returns:

Model value [m]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Nye (1955) with the water pressure added (Benn et al. 2007); the result is meaningful only while it is less than the ice thickness and at least the water depth. Equation: 2 * (strain_rate / rate_factor) ** (1 / exponent) / (ice_density * gravity) + water_density / ice_density * water_depth. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.

earthsciencekit.glaciology.weertman_sliding_velocity(basal_stress, sliding_coefficient, exponent)

Weertman’s sliding velocity of temperate ice over a rough bed.

Ice gets past bed obstacles by regelation around small ones and enhanced creep around large ones, which makes the sliding speed a power of the basal stress.

Parameters:

• basal_stress (float or array_like) – Pa; constraint: nonnegative.

• sliding_coefficient (float or array_like) – m/(s Pa^m); constraint: positive.

• exponent (float or array_like) – 1, m = (n + 1) / 2, about 2 for Glen’s n = 3; constraint: at least 1.

Returns:

Model value [m/s]. Scalar input returns a NumPy scalar; arrays broadcast.

Return type:

numpy scalar or ndarray

Notes

Weertman (1957): a bed at the melting point, free of cavities and sediment; the coefficient lumps the bed roughness. Equation: sliding_coefficient * basal_stress**exponent. Invalid, empty, nonreal or nonfinite inputs raise ValueError. NumPy broadcasting applies; singular and overflowing outputs raise ValueError.