earths 0.0.1

High-fidelity Earth simulation engine — orbit, atmosphere, geology, hydrology, biosphere, terrain, lighting, rendering, satellites, and temporal systems with full scientific coupling
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102

use sciforge::hub::domain::common::constants::{EARTH_MASS, EARTH_RADIUS, G, R_GAS};
use sciforge::hub::prelude::constants::elements::atomic_mass;
pub fn ice_density() -> f64 {
    *crate::ICE_DENSITY
}
pub const ICE_CREEP_EXPONENT: f64 = 3.0;
pub const GLEN_FLOW_LAW_A: f64 = 2.4e-24;
pub fn water_molar_mass_from_elements() -> f64 {
    2.0 * atomic_mass(1) + atomic_mass(8)
}
pub struct Glacier {
    pub name: &'static str,
    pub length_km: f64,
    pub thickness_m: f64,
    pub surface_slope_deg: f64,
    pub accumulation_m_yr: f64,
    pub ablation_m_yr: f64,
}
impl Glacier {
    pub fn mass_balance_m_yr(&self) -> f64 {
        self.accumulation_m_yr - self.ablation_m_yr
    }
    pub fn basal_shear_stress_pa(&self) -> f64 {
        let g = G * EARTH_MASS / (EARTH_RADIUS * EARTH_RADIUS);
        let slope_rad = self.surface_slope_deg.to_radians();
        ice_density() * g * self.thickness_m * slope_rad.sin()
    }
    pub fn deformation_velocity_m_yr(&self) -> f64 {
        let tau = self.basal_shear_stress_pa();
        let n = ICE_CREEP_EXPONENT;
        let velocity_m_s = 2.0 * GLEN_FLOW_LAW_A * tau.powf(n) * self.thickness_m / (n + 1.0);
        velocity_m_s * crate::SECONDS_PER_YEAR
    }
    pub fn equilibrium_line_altitude(&self, base_elevation_m: f64) -> f64 {
        base_elevation_m
            + self.length_km * 1000.0 * 0.5 * (self.surface_slope_deg.to_radians()).tan()
    }
    pub fn volume_m3(&self, width_m: f64) -> f64 {
        self.length_km * 1000.0 * width_m * self.thickness_m * 0.7
    }
    pub fn sea_level_equivalent_mm(&self, width_m: f64) -> f64 {
        let vol = self.volume_m3(width_m);
        let ocean_area = crate::OCEAN_SURFACE_AREA;
        vol * ice_density() / 1000.0 / ocean_area * 1000.0
    }
}
pub fn antarctic_ice_sheet() -> Glacier {
    Glacier {
        name: "Antarctic Ice Sheet",
        length_km: 4500.0,
        thickness_m: 2160.0,
        surface_slope_deg: 0.1,
        accumulation_m_yr: 0.166,
        ablation_m_yr: 0.150,
    }
}
pub fn greenland_ice_sheet() -> Glacier {
    Glacier {
        name: "Greenland Ice Sheet",
        length_km: 2400.0,
        thickness_m: 1500.0,
        surface_slope_deg: 0.3,
        accumulation_m_yr: 0.34,
        ablation_m_yr: 0.32,
    }
}
pub fn ice_rheology_viscosity(temperature_c: f64, stress_pa: f64) -> f64 {
    let temp_k = temperature_c + crate::CELSIUS_TO_KELVIN;
    let q = crate::Q_ICE_CREEP;
    let a = GLEN_FLOW_LAW_A * (-q / (R_GAS * temp_k)).exp();
    1.0 / (2.0 * a * stress_pa.powf(ICE_CREEP_EXPONENT - 1.0))
}
pub fn glen_strain_rate(temperature_c: f64, stress_pa: f64) -> f64 {
    let temp_k = temperature_c + crate::CELSIUS_TO_KELVIN;
    let q = crate::Q_ICE_CREEP;
    let a_temp = GLEN_FLOW_LAW_A * (-q / (R_GAS * temp_k)).exp();
    a_temp * stress_pa.powf(ICE_CREEP_EXPONENT)
}
pub fn shallow_ice_velocity(thickness_m: f64, slope_rad: f64, temperature_c: f64) -> f64 {
    let temp_k = temperature_c + crate::CELSIUS_TO_KELVIN;
    let q = crate::Q_ICE_CREEP;
    let a_temp = GLEN_FLOW_LAW_A * (-q / (R_GAS * temp_k)).exp();
    let g = G * EARTH_MASS / (EARTH_RADIUS * EARTH_RADIUS);
    let n = ICE_CREEP_EXPONENT;
    let driving_stress = ice_density() * g * slope_rad.sin().abs();
    2.0 * a_temp / (n + 1.0) * driving_stress.powf(n) * thickness_m.powf(n + 1.0)
}
pub fn glacial_bed_erosion_rate(
    basal_velocity_m_s: f64,
    ice_thickness_m: f64,
    water_pressure_pa: f64,
) -> f64 {
    const K_G: f64 = 1e-4;
    const L_HALLET: f64 = 2.0;
    const P_STRESS_EXP: f64 = 1.0;
    const N_REF_PA: f64 = 1.0e6;
    let g = G * EARTH_MASS / (EARTH_RADIUS * EARTH_RADIUS);
    let overburden = ice_density() * g * ice_thickness_m;
    let n_eff = (overburden - water_pressure_pa).max(0.0);
    let stress_ratio = (n_eff / N_REF_PA).powf(P_STRESS_EXP);
    K_G * basal_velocity_m_s.abs().powf(L_HALLET) * stress_ratio
}