earths 0.0.4

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

use sciforge::hub::domain::common::constants::{EARTH_MASS, EARTH_RADIUS, G, R_GAS};
use sciforge::hub::prelude::constants::elements::atomic_mass;
pub fn icedensity() -> f64 {
    *crate::ICEDENSITY
}
pub const ICECREEPEXPONENT: f64 = 3.0;
pub const GLENFLOWLAWA: f64 = 2.4e-24;
pub fn watermolarmassfromelements() -> f64 {
    2.0 * atomic_mass(1) + atomic_mass(8)
}
pub struct Glacier {
    pub name: &'static str,
    pub lengthkm: f64,
    pub thicknessm: f64,
    pub surfaceslopedeg: f64,
    pub accumulationmyr: f64,
    pub ablationmyr: f64,
}
impl Glacier {
    pub fn massbalancemyr(&self) -> f64 {
        self.accumulationmyr - self.ablationmyr
    }
    pub fn basalshearstresspa(&self) -> f64 {
        let g = G * EARTH_MASS / (EARTH_RADIUS * EARTH_RADIUS);
        let sloperad = self.surfaceslopedeg.to_radians();
        icedensity() * g * self.thicknessm * sloperad.sin()
    }
    pub fn deformationvelocitymyr(&self) -> f64 {
        let tau = self.basalshearstresspa();
        let n = ICECREEPEXPONENT;
        let velocityms = 2.0 * GLENFLOWLAWA * tau.powf(n) * self.thicknessm / (n + 1.0);
        velocityms * crate::SECONDSPERYEAR
    }
    pub fn equilibriumlinealtitude(&self, baseelevationm: f64) -> f64 {
        baseelevationm + self.lengthkm * 1000.0 * 0.5 * (self.surfaceslopedeg.to_radians()).tan()
    }
    pub fn volumem3(&self, widthm: f64) -> f64 {
        self.lengthkm * 1000.0 * widthm * self.thicknessm * 0.7
    }
    pub fn sealevelequivalentmm(&self, widthm: f64) -> f64 {
        let vol = self.volumem3(widthm);
        let oceanarea = crate::OCEANSURFACEAREA;
        vol * icedensity() / 1000.0 / oceanarea * 1000.0
    }
}
pub fn antarcticicesheet() -> Glacier {
    Glacier {
        name: "Antarctic Ice Sheet",
        lengthkm: 4500.0,
        thicknessm: 2160.0,
        surfaceslopedeg: 0.1,
        accumulationmyr: 0.166,
        ablationmyr: 0.150,
    }
}
pub fn greenlandicesheet() -> Glacier {
    Glacier {
        name: "Greenland Ice Sheet",
        lengthkm: 2400.0,
        thicknessm: 1500.0,
        surfaceslopedeg: 0.3,
        accumulationmyr: 0.34,
        ablationmyr: 0.32,
    }
}
pub fn icerheologyviscosity(temperaturec: f64, stresspa: f64) -> f64 {
    let tempk = temperaturec + crate::CELSIUSTOKELVIN;
    let q = crate::QICECREEP;
    let a = GLENFLOWLAWA * (-q / (R_GAS * tempk)).exp();
    1.0 / (2.0 * a * stresspa.powf(ICECREEPEXPONENT - 1.0))
}
pub fn glenstrainrate(temperaturec: f64, stresspa: f64) -> f64 {
    let tempk = temperaturec + crate::CELSIUSTOKELVIN;
    let q = crate::QICECREEP;
    let atemp = GLENFLOWLAWA * (-q / (R_GAS * tempk)).exp();
    atemp * stresspa.powf(ICECREEPEXPONENT)
}
pub fn shallowicevelocity(thicknessm: f64, sloperad: f64, temperaturec: f64) -> f64 {
    let tempk = temperaturec + crate::CELSIUSTOKELVIN;
    let q = crate::QICECREEP;
    let atemp = GLENFLOWLAWA * (-q / (R_GAS * tempk)).exp();
    let g = G * EARTH_MASS / (EARTH_RADIUS * EARTH_RADIUS);
    let n = ICECREEPEXPONENT;
    let drivingstress = icedensity() * g * sloperad.sin().abs();
    2.0 * atemp / (n + 1.0) * drivingstress.powf(n) * thicknessm.powf(n + 1.0)
}
pub fn glacialbederosionrate(
    basalvelocityms: f64,
    icethicknessm: f64,
    waterpressurepa: f64,
) -> f64 {
    const KG: f64 = 1e-4;
    const LHALLET: f64 = 2.0;
    const PSTRESSEXP: f64 = 1.0;
    const NREFPA: f64 = 1.0e6;
    let g = G * EARTH_MASS / (EARTH_RADIUS * EARTH_RADIUS);
    let overburden = icedensity() * g * icethicknessm;
    let neff = (overburden - waterpressurepa).max(0.0);
    let stressratio = (neff / NREFPA).powf(PSTRESSEXP);
    KG * basalvelocityms.abs().powf(LHALLET) * stressratio
}