diff --git a/chained_eclipse/models.py b/chained_eclipse/models.py
index bd07287..0bf30d0 100644
--- a/chained_eclipse/models.py
+++ b/chained_eclipse/models.py
@@ -2,6 +2,7 @@
 
 from __future__ import annotations
 
+import math
 from dataclasses import asdict, dataclass, field
 from typing import Any, Callable
 
@@ -14,11 +15,25 @@
 )
 
 
+def spherical_mass_kg(radius_km: float, density_kg_m3: float) -> float:
+    """Return the mass of a spherical body from its radius and bulk density."""
+
+    if not math.isfinite(radius_km) or radius_km <= 0.0:
+        raise ValueError("radius_km must be a finite positive value")
+    if not math.isfinite(density_kg_m3) or density_kg_m3 <= 0.0:
+        raise ValueError("density_kg_m3 must be a finite positive value")
+    return 4.0 / 3.0 * math.pi * (radius_km * 1_000.0) ** 3 * density_kg_m3
+
+
 @dataclass(slots=True)
 class OrbitalElements:
     """Osculating Earth-centred elements at the configured epoch.
 
-    Angular fields are degrees in the mean J2000 ecliptic frame.
+    Angular fields are degrees in the mean J2000 ecliptic frame. When callers
+    customize the radius or density but leave ``mass_kg`` at the package
+    default, the mass is recalculated to keep the physical parameters
+    internally consistent. An explicitly supplied non-default mass is retained
+    for independently constrained bodies and massless control integrations.
     """
 
     semimajor_axis_km: float = 180_000.0
@@ -33,6 +48,17 @@ class OrbitalElements:
     mass_kg: float = SECOND_MOON_MASS_KG
     frame: str = "mean_ecliptic_J2000"
 
+    def __post_init__(self) -> None:
+        derived_mass = spherical_mass_kg(self.radius_km, self.density_kg_m3)
+        uses_custom_bulk_properties = (
+            self.radius_km != SECOND_MOON_RADIUS_KM
+            or self.density_kg_m3 != SECOND_MOON_DENSITY_KG_M3
+        )
+        if uses_custom_bulk_properties and self.mass_kg == SECOND_MOON_MASS_KG:
+            self.mass_kg = derived_mass
+        elif not math.isfinite(self.mass_kg) or self.mass_kg < 0.0:
+            raise ValueError("mass_kg must be a finite non-negative value")
+
     def to_dict(self) -> dict[str, Any]:
         return asdict(self)
 
@@ -136,4 +162,3 @@ def position(self, jd: np.ndarray | float) -> np.ndarray:
         if state.ndim == 1:
             return state[:3]
         return state[..., :3]
-
diff --git a/tests/test_orbital_dynamics.py b/tests/test_orbital_dynamics.py
index 3db0831..a34ddc3 100644
--- a/tests/test_orbital_dynamics.py
+++ b/tests/test_orbital_dynamics.py
@@ -1,10 +1,12 @@
 """Checks for mass calculation, element conversion, and two-body propagation."""
 
+import math
+
 import numpy as np
 import pytest
 
 from chained_eclipse.constants import SECOND_MOON_MASS_KG
-from chained_eclipse.models import OrbitalElements
+from chained_eclipse.models import OrbitalElements, spherical_mass_kg
 from chained_eclipse.orbital_dynamics import (
     elements_to_state,
     mean_motion_rad_s,
@@ -17,6 +19,52 @@ def test_second_moon_mass_matches_radius_density_calculation() -> None:
     assert SECOND_MOON_MASS_KG == pytest.approx(8.2430310159e21, rel=1e-10)
 
 
+def test_custom_bulk_properties_recalculate_default_mass() -> None:
+    elements = OrbitalElements(radius_km=1_000.0, density_kg_m3=2_500.0)
+
+    assert elements.mass_kg == pytest.approx(spherical_mass_kg(1_000.0, 2_500.0))
+    assert elements.mass_kg != SECOND_MOON_MASS_KG
+
+
+def test_explicit_independently_constrained_mass_is_retained() -> None:
+    elements = OrbitalElements(
+        radius_km=1_000.0,
+        density_kg_m3=2_500.0,
+        mass_kg=9.0e21,
+    )
+
+    assert elements.mass_kg == 9.0e21
+
+
+def test_massless_control_is_retained() -> None:
+    elements = OrbitalElements(mass_kg=0.0)
+
+    assert elements.mass_kg == 0.0
+
+
+@pytest.mark.parametrize(
+    ("radius_km", "density_kg_m3"),
+    [
+        (0.0, 2_500.0),
+        (-1.0, 2_500.0),
+        (math.inf, 2_500.0),
+        (1_000.0, 0.0),
+        (1_000.0, -1.0),
+        (1_000.0, math.nan),
+    ],
+)
+def test_invalid_bulk_properties_are_rejected(
+    radius_km: float, density_kg_m3: float
+) -> None:
+    with pytest.raises(ValueError):
+        OrbitalElements(radius_km=radius_km, density_kg_m3=density_kg_m3)
+
+
+def test_invalid_explicit_mass_is_rejected() -> None:
+    with pytest.raises(ValueError, match="mass_kg"):
+        OrbitalElements(mass_kg=-1.0)
+
+
 def test_element_state_round_trip() -> None:
     elements = OrbitalElements(
         longitude_ascending_node_deg=119.7,
@@ -43,4 +91,3 @@ def test_two_body_state_repeats_after_one_orbital_period() -> None:
     final = propagate_two_body(elements, period)
     assert np.linalg.norm(final[:3] - initial[:3]) < 1e-6
     assert np.linalg.norm(final[3:] - initial[3:]) < 1e-10
-