ground_motion_lib/
mf2013.rs

1//! Implementation of Morikawa & Fujiwara (2013) Ground Motion Prediction Equations (GMPE).
2//!
3//! This module defines the parameters and calculation logic for predicting
4//! ground motion values (PGA, PGV, PSA) based on earthquake and site characteristics.
5
6use crate::auxilary::{DL, G_GLOBAL};
7use crate::gmm::{Earthquake, GmpePoint, GmpePointKind, GroundMotionModeling, Vs30Point};
8use geo::{Distance, Haversine, Point};
9
10/// Morikawa & Fujiwara (2013) Ground Motion Prediction Equation parameters.
11#[derive(Debug)]
12pub struct MF2013 {
13    /// Magnitude upper limit (Mw0)
14    pub mw0: f64,
15    /// Coefficient for magnitude scaling
16    pub a: f64,
17    /// Coefficient for distance scaling
18    pub b: f64,
19    /// Constant term
20    pub c: f64,
21    /// Distance damping parameter
22    pub d: f64,
23    /// Exponent scaling factor for distance damping
24    pub e: f64,
25    /// Standard deviation of the log ground motion (not currently used)
26    pub sigma: f64,
27    /// Coefficient for deep sedimentary layer correction
28    pub pd: f64,
29    /// Minimum depth for deep sedimentary layer correction
30    pub dl_min: f64,
31    /// Reference depth for deep layer correction
32    pub d0: f64,
33    /// Coefficient for Vs30 amplification term
34    pub ps: f64,
35    /// Maximum Vs30 considered for amplification (Vs_max)
36    pub vs_max: f64,
37    /// Reference Vs30 value (V0)
38    pub v0: f64,
39    /// Coefficient for anomalous seismic intensity distribution (ASID)
40    pub gamma: f64,
41    /// Whether ASID correction is enabled
42    pub asid: bool,
43    /// Type of motion (PGA, PGV, PSA etc.)
44    pub motion_kind: GmpePointKind,
45}
46
47impl MF2013 {
48    /// Calculate predicted ground motion value (in physical units) for a site and earthquake.
49    ///
50    /// Note: Currently assumes a point source (no finite fault modeling).
51    ///
52    /// # Arguments
53    ///
54    /// * `epicentral_distance` - Horizontal distance from the site to the earthquake epicenter (km).
55    /// * `eq_mag` - Earthquake moment magnitude (Mw).
56    /// * `eq_depth` - Hypocentral depth (km).
57    /// * `vs_30` - Average shear-wave velocity in the top 30 meters at the site (m/s).
58    /// * `dl` - Depth to the 1400 m/s shear-wave velocity layer (m).
59    /// * `xvf` - Binary flag for volcanic front effect (1.0 if oceanward of front, 0.0 otherwise).
60    ///
61    /// # Returns
62    ///
63    /// Predicted ground motion value in cm/s² (PGA, PSA) or cm/s (PGV).
64    fn get_gmpe_by_distnace(
65        &self,
66        epicentral_distance: f64,
67        eq_mag: f64,
68        eq_depth: f64,
69        vs_30: f64,
70        dl: f64,
71        xvf: f64,
72    ) -> f64 {
73        // Rupture distance assuming point source
74        let r_rup = (epicentral_distance.powi(2) + eq_depth.powi(2)).sqrt();
75
76        let magnitude = eq_mag.min(self.mw0);
77        let a_m_w = self.a * magnitude;
78
79        // Deep sedimentary layer correction
80        let g_d = self.pd * (dl.max(self.dl_min) / self.d0).log10();
81
82        // Main GMPE equation (log10 of predicted motion)
83        // logA where A in cm/s^2 (pga,psa) or cm/s (pgv)
84        let log_a = (a_m_w + self.b * r_rup + self.c)
85            - (r_rup + self.d * 10.0_f64.powf(self.e * magnitude)).log10();
86
87        // Amplification by Deep Sedimentary Layers
88        // Apply deep layer correction
89        let log_agd = log_a + g_d;
90
91        // Vs30 site amplification
92        let gs = self.ps * (vs_30.min(self.vs_max) / self.v0).log10();
93        let log_ags = log_agd + gs;
94
95        // Optional anomalous seismic intensity distribution correction
96        if self.asid {
97            let ai = self.gamma + xvf * (eq_depth - 30.);
98            10.0_f64.powf(log_ags + ai)
99        } else {
100            10.0_f64.powf(log_ags)
101        }
102    }
103}
104
105impl GroundMotionModeling for MF2013 {
106    /// Compute ground motion prediction at a given site point for a specified earthquake event.
107    ///
108    /// # Arguments
109    ///
110    /// * `point` - The site location and properties (longitude, latitude, Vs30, depth to 1400 m/s
111    ///   layer, etc.).
112    /// * `eq` - The earthquake event (magnitude, depth, hypocenter location).
113    ///
114    /// # Returns
115    ///
116    /// A `GmpePoint` containing the predicted ground motion value and associated metadata.
117    fn calc_from_point(&self, point: &Vs30Point, eq: &Earthquake) -> GmpePoint {
118        let epicentral_distance = Haversine
119            .distance(Point::new(eq.lon, eq.lat), Point::new(point.lon, point.lat))
120            / 1000.;
121        let vs_30 = point.vs30;
122        let dl = match point.dl {
123            None => DL as f64,
124            Some(dl) => dl,
125        };
126        let xvf = match point.xvf {
127            None => 0.,
128            Some(_) => 1.,
129        };
130        let mut ground_motion =
131            self.get_gmpe_by_distnace(epicentral_distance, eq.magnitude, eq.depth, vs_30, dl, xvf);
132        // convert cm/c^2 to %g
133        if matches!(self.motion_kind, GmpePointKind::Pga | GmpePointKind::Psa) {
134            ground_motion = ((ground_motion / 100.) / G_GLOBAL) * 100.;
135        };
136        GmpePoint {
137            lon: point.lon,
138            lat: point.lat,
139            value: ground_motion,
140            kind: self.motion_kind,
141        }
142    }
143}
ground_motion_lib/mf2013.rs

ground_motion_lib/
mf2013.rs
