Use BoldM type for fault magnitudes in Magnitudes config

pyproject.toml

@@ -16,7 +16,7 @@ dependencies = [
   "nshmdb>=2025.12.1",
   "oq_wrapper>=2025.12.3",
   "qcore-utils>=2025.12.2",
-  "source_modelling>=2025.12.1",
+  "source_modelling>=2026.6.2",
   # Data Formats
   "geopandas",
   "pandas[parquet, hdf5]",

tests/test_generate_domain.py

@@ -6,7 +6,7 @@
 from hypothesis import given
 from hypothesis import strategies as st
 
-from source_modelling import sources
+from source_modelling import magnitude_scaling, sources
 from workflow.realisations import (
     Magnitudes,
     Rakes,
@@ -90,7 +90,7 @@ def test_generate_domain() -> None:
         dip_dir=180.0,
     )
     source_config = SourceConfig(dict(source=source))
-    magnitudes = Magnitudes(dict(source=6.0))
+    magnitudes = Magnitudes(dict(source=magnitude_scaling.BoldM(6.0)))
     rakes = Rakes(dict(source=180.0))
 
     velocity_model_parameters = VelocityModelParameters(

tests/test_realisation.py

@@ -13,7 +13,7 @@
 import schema
 
 from IM import im_calculation
-from source_modelling import rupture_propagation
+from source_modelling import magnitude_scaling, rupture_propagation
 from velocity_modelling import bounding_box
 from workflow import defaults, realisations, schemas
 from workflow.realisations import SourceConfig
@@ -493,7 +493,11 @@ def test_rupture_prop_config(tmp_path: Path) -> None:
 
 def test_magnitudes(tmp_path: Path) -> None:
     magnitudes = realisations.Magnitudes(
-        magnitudes={"A": 6.5, "B": 6.7, "C": 6.9},
+        magnitudes={
+            "A": magnitude_scaling.BoldM(6.5),
+            "B": magnitude_scaling.BoldM(6.7),
+            "C": magnitude_scaling.BoldM(6.9),
+        },
     )
 
     realisation_ffp = tmp_path / "realisation.json"

tests/test_realisation_to_srf.py

@@ -158,7 +158,7 @@ def test_build_genslip_command_static_args() -> None:
         "nh=1",
         "read_erf=0",
         "plane_header=1",
-        "srf_version=1.0",
+        "srf_version=2.0",
         "read_gsf=1",
         "nstk=50",
         "ndip=25",

workflow/realisations.py

@@ -26,6 +26,7 @@
 
 from IM import im_calculation
 from source_modelling import sources
+from source_modelling.magnitude_scaling import BoldM
 from source_modelling.rupture_propagation import JumpPair
 from source_modelling.sources import IsSource
 from velocity_modelling.bounding_box import BoundingBox
@@ -649,10 +650,10 @@ class Magnitudes(RealisationConfiguration):
     _config_key: ClassVar[str] = "magnitudes"
     _schema: ClassVar[Schema] = schemas.MAGNITUDE_SCHEMA
 
-    magnitudes: dict[str, float]
+    magnitudes: dict[str, BoldM]
     """A map from faults to their magnitudes."""
 
-    def __getitem__(self, key: str) -> float:
+    def __getitem__(self, key: str) -> BoldM:
         """Get the magnitude for a fault name.
 
         Parameters
@@ -662,7 +663,7 @@ def __getitem__(self, key: str) -> float:
 
         Returns
         -------
-        float
+        BoldM
             The magnitude.
         """
         return self.magnitudes[key]

workflow/scripts/gcmt_to_realisation.py

@@ -192,8 +192,7 @@ def gcmt_to_realisation(
             param_hint="GCMT_EVENT_ID",
         )
 
-    magnitude = moment.moment_to_magnitude(solution_moment)
-
+    magnitude = moment.moment_to_magnitude(solution_moment, bold_m=True)
     model = community_fault_model.get_community_fault_model()
 
     match nodal_plane:
@@ -206,43 +205,34 @@ def gcmt_to_realisation(
                 model, np.array([latitude, longitude]), nodal_plane_1, nodal_plane_2
             )
 
-    rake = selected_nodal_plane.rake
-
     # Calculate dip direction from strike (strike + 90 degrees for right-hand rule)
     dip_direction = (selected_nodal_plane.strike + 90) % 360
 
     length, width = magnitude_scaling.magnitude_to_length_width(
-        scaling_relation, magnitude, rake
+        scaling_relation, magnitude, selected_nodal_plane.rake
     )
 
     centroid = np.array([latitude, longitude, centroid_depth])
 
     # Create source based on source_type parameter
     if source_type == SourceType.POINT_SOURCE:
-        length_km, width_km = magnitude_scaling.magnitude_to_length_width(
-            scaling_relation, magnitude, rake
-        )
-        length_m = length_km * 1000  # Convert km to meters
-        width_m = width_km * 1000  # Convert km to meters
-
         source_geometry = sources.Point.from_lat_lon_depth(
             point_coordinates=np.array(
                 [latitude, longitude, centroid_depth * 1000]
             ),  # Convert km to meters
-            length_m=length_m,  # Use calculated length from area
-            width_m=width_m,  # Use calculated width from area
+            length_m=length*1000,  # convert km to metres
+            width_m=width*1000,  # convert km to metres
             strike=selected_nodal_plane.strike,
             dip=selected_nodal_plane.dip,
             dip_dir=dip_direction,
         )
 
     else:
-        # Create plane source (default behavior)
         plane = sources.Plane.from_centroid_strike_dip(
             centroid,
             selected_nodal_plane.dip,
-            length,
-            width,
+            length*1000,  # convert km to metres
+            width*1000,  # convert km to metres
             strike=selected_nodal_plane.strike,
         )
 
@@ -284,8 +274,10 @@ def gcmt_to_realisation(
         hypocentre = np.array([1 / 2, 1 / 2])
 
     source_config = SourceConfig(source_geometries={gcmt_event_id: source_geometry})
-    magnitudes = Magnitudes(magnitudes={gcmt_event_id: float(magnitude)})
-    rakes = Rakes(rakes={gcmt_event_id: float(rake)})
+    magnitudes = Magnitudes(
+        magnitudes={gcmt_event_id: magnitude}
+    )
+    rakes = Rakes(rakes={gcmt_event_id: float(selected_nodal_plane.rake)})
     rupture_config = RupturePropagationConfig(
         rupture_causality_tree={gcmt_event_id: None},
         jump_points={},

workflow/scripts/generate_domain.py

@@ -38,7 +38,7 @@
 import typer
 
 from qcore import cli, geo
-from source_modelling import moment, sources
+from source_modelling import magnitude_scaling, moment, sources
 from velocity_modelling import bounding_box
 from velocity_modelling.bounding_box import BoundingBox
 from workflow import log_utils, realisations, utils
@@ -129,13 +129,13 @@ def boundary_distance(
     return float(shapely.hausdorff_distance(source_geometry, bounding_box_geometry))
 
 
-def total_magnitude(magnitudes: Iterable[float]) -> float:
+def total_magnitude(magnitudes: Iterable[magnitude_scaling.BoldM]) -> float:
     """
     Compute the total magnitude from an array of individual magnitudes.
 
     Parameters
     ----------
-    magnitudes : Iterable[float]
+    magnitudes : Iterable[magnitude_scaling.BoldM]
         An array of magnitudes.
 
     Returns
@@ -144,7 +144,7 @@ def total_magnitude(magnitudes: Iterable[float]) -> float:
         The total magnitude, computed from the summed moment of the input magnitudes.
     """
     total_moment = sum(
-        moment.magnitude_to_moment(magnitude) for magnitude in magnitudes
+        moment.magnitude_to_moment(magnitude, bold_m=True) for magnitude in magnitudes
     )
     return moment.moment_to_magnitude(total_moment)
 
@@ -228,7 +228,7 @@ def average_rake(rakes: Rakes, magnitudes: Magnitudes) -> float:
         The moment-weighted average rake.
     """
     moments = {
-        k: moment.magnitude_to_moment(magnitude)
+        k: moment.magnitude_to_moment(magnitude, bold_m=True)
         for k, magnitude in magnitudes.magnitudes.items()
     }
     max_moment = max(moments.values())

workflow/scripts/nshm2022_to_realisation.py

@@ -45,8 +45,8 @@
 
 from nshmdb import nshmdb
 from qcore import cli
-from qcore.uncertainties import distributions, mag_scaling
-from source_modelling import moment, rupture_propagation, sources
+from qcore.uncertainties import distributions
+from source_modelling import magnitude_scaling, moment, rupture_propagation, sources
 from source_modelling.sources import Fault
 from workflow import realisations
 from workflow.defaults import DefaultsVersion
@@ -64,39 +64,13 @@
 app = typer.Typer()
 
 
-def a_to_mw_leonard(area: float, rake: float) -> float:
-    """
-    Convert fault area and rake to moment magnitude using the Leonard scaling relation.
-
-    Parameters
-    ----------
-    area : float
-        The area of the fault in square kilometres.
-    rake : float
-        The rake angle of the fault in degrees.
-
-    Returns
-    -------
-    float
-        The estimated moment magnitude of the fault.
-
-    References
-    ----------
-    Leonard, M. (2010). Earthquake fault scaling: Self-consistent
-    relating of rupture length, width, average displacement, and
-    moment release. Bulletin of the Seismological Society of America,
-    100(5A), 1971-1988.
-    """
-    return mag_scaling.a_to_mw_leonard(area, 4, 3.99, rake)
-
-
 def default_magnitude_estimation(
     faults: dict[str, Fault],
     # NOTE: this must be in quotes because the runtime class DisjointSet is
     # not generic, just the stub implementation.
     components: "DisjointSet[str]",
     avg_rake: float,
-) -> dict[str, float]:
+) -> dict[str, magnitude_scaling.BoldM]:
     """Estimate the magnitudes for a set of faults based on their areas and average rake.
 
     Parameters
@@ -113,11 +87,14 @@ def default_magnitude_estimation(
     Returns
     -------
     dict
-        A dictionary where the keys are fault names and the values are the estimated magnitudes for each fault.
+        A dictionary where the keys are fault names and the values are the
+        estimated magnitudes (in the `BoldM` convention) for each fault.
     """
     total_area = sum(fault.area() for fault in faults.values())
-    estimated_mw = a_to_mw_leonard(total_area, avg_rake)
-    estimated_moment = mag_scaling.mag2mom(estimated_mw)
+    estimated_magnitude = magnitude_scaling.area_to_magnitude(
+        magnitude_scaling.ScalingRelation.LEONARD2014, total_area, avg_rake
+    )
+    estimated_moment = moment.magnitude_to_moment(estimated_magnitude, bold_m=True)
     roots = {components[fault_name] for fault_name in faults}
     component_areas = {
         root: sum(faults[name].area() for name in components.subset(root))
@@ -139,7 +116,7 @@ def default_magnitude_estimation(
         segment_moments[fault_name] = (fault.area() / component_area) * component_moment
 
     return {
-        fault_name: mag_scaling.mom2mag(fault_moment)
+        fault_name: moment.moment_to_magnitude(fault_moment, bold_m=True)
         for fault_name, fault_moment in segment_moments.items()
     }
 
@@ -355,7 +332,9 @@ def generate_realisation(
             ]
         )
     else:
-        mfds_rates = db.most_likely_fault(rupture_id, magnitudes)
+        # The ty ignore below can be removed once NSHM2022DB merges the change to
+        # accept dict[str, BoldM] in most_likely_fault (branch support-BoldM-in-workflow).
+        mfds_rates = db.most_likely_fault(rupture_id, magnitudes)  # ty: ignore[invalid-argument-type]
         mfds_probabilities = np.array(list(mfds_rates.values()))
         if np.allclose(mfds_probabilities, 0):
             mfds_probabilities = np.ones_like(mfds_probabilities)

workflow/scripts/realisation_to_srf.py

@@ -456,7 +456,7 @@ def _build_genslip_command(
     cmd = [
         str(genslip_path),
         "plane_header=1",
-        "srf_version=1.0",
+        "srf_version=2.0",
         "read_erf=0",
         "write_srf=1",
         "read_gsf=1",
@@ -657,7 +657,7 @@ def generate_point_source_srf(
 
     # Get magnitude and convert to seismic moment
     magnitude = params.magnitudes.magnitudes[name]
-    moment_newton_metre = moment.magnitude_to_moment(magnitude)
+    moment_newton_metre = moment.magnitude_to_moment(magnitude, bold_m=True)
 
     velocity_model_df = params.velocity_model_1d.model.copy()
     velocity_model_df["depth_km"] = (