diff --git a/flow360/__init__.py b/flow360/__init__.py
index 6961851e8..124ef8ddb 100644
--- a/flow360/__init__.py
+++ b/flow360/__init__.py
@@ -118,6 +118,7 @@
     Fluid,
     ForcePerArea,
     FromUserDefinedDynamics,
+    Gravity,
     HeatEquationInitialCondition,
     NavierStokesInitialCondition,
     NavierStokesModifiedRestartSolution,
diff --git a/flow360/component/simulation/exposed_units.py b/flow360/component/simulation/exposed_units.py
index 8461f77ed..13b552213 100644
--- a/flow360/component/simulation/exposed_units.py
+++ b/flow360/component/simulation/exposed_units.py
@@ -15,6 +15,7 @@
     "temperature": [],
     "delta_temperature": [],
     "velocity": [],
+    "acceleration": [],
     "area": [],
     "force": [],
     "pressure": [],
@@ -71,6 +72,7 @@
     "(length)**2": {"SI": "m**2", "CGS": "cm**2", "Imperial": "ft**2"},
     "(length)**(-2)": {"SI": "1/m**2", "CGS": "1/cm**2", "Imperial": "1/ft**2"},
     "(length)/(time)": {"SI": "m/s", "CGS": "cm/s", "Imperial": "ft/s"},
+    "(length)/(time)**2": {"SI": "m/s**2", "CGS": "cm/s**2", "Imperial": "ft/s**2"},  # Acceleration
     "(angle)/(time)": ["rad/s", "degree/s", "rpm"],  # list --> Unit system agnostic dimensions
     "(angle)": ["degree", "rad"],
     "(temperature)": {"SI": "K", "CGS": "K", "Imperial": "degF"},
diff --git a/flow360/component/simulation/meshing_param/volume_params.py b/flow360/component/simulation/meshing_param/volume_params.py
index 50bbe9ca8..eb53ca4f2 100644
--- a/flow360/component/simulation/meshing_param/volume_params.py
+++ b/flow360/component/simulation/meshing_param/volume_params.py
@@ -661,7 +661,7 @@ class AutomatedFarfield(_FarfieldBase):
     def farfield(self):
         """Returns the farfield boundary surface."""
         # Make sure the naming is the same here and what the geometry/surface mesh pipeline generates.
-        return GhostSurface(name="farfield", private_attribute_id="farfield")
+        return GhostSurface(name="farField", private_attribute_id="farField")
 
     @property
     def symmetry_plane(self) -> GhostSurface:
diff --git a/flow360/component/simulation/models/volume_models.py b/flow360/component/simulation/models/volume_models.py
index a46a390e6..9ce359029 100644
--- a/flow360/component/simulation/models/volume_models.py
+++ b/flow360/component/simulation/models/volume_models.py
@@ -63,6 +63,7 @@
     SeedpointVolume,
 )
 from flow360.component.simulation.unit_system import (
+    AccelerationType,
     AngleType,
     AngularVelocityType,
     HeatSourceType,
@@ -1485,6 +1486,53 @@ def _validate_volumetric_heat_source_for_liquid(
         return value
 
 
+class Gravity(Flow360BaseModel):
+    """
+    :class:`Gravity` class for specifying gravitational body force.
+
+    The gravity model applies a body force ρg to the momentum equations and ρ(g·u) to the
+    energy equation, enabling simulation of buoyancy-driven flows.
+
+    Example
+    -------
+
+    Define gravity acting in the negative z-direction with Earth's gravitational acceleration:
+
+    >>> fl.Gravity(
+    ...     direction=(0, 0, -1),
+    ...     magnitude=9.81 * fl.u.m / fl.u.s**2,
+    ... )
+
+    Define gravity applied only to specific volume zones:
+
+    >>> fl.Gravity(
+    ...     entities=[volume_mesh["fluid_zone"]],
+    ...     direction=(0, 0, -1),
+    ...     magnitude=9.81 * fl.u.m / fl.u.s**2,
+    ... )
+
+    ====
+    """
+
+    name: Optional[str] = pd.Field("Gravity", description="Name of the `Gravity` model.")
+    type: Literal["Gravity"] = pd.Field("Gravity", frozen=True)
+    entities: Optional[EntityList[GenericVolume]] = pd.Field(
+        default=None,
+        alias="volumes",
+        description="The entity list for the `Gravity` model. "
+        + "If not specified, gravity is applied to all fluid zones.",
+    )
+    direction: Axis = pd.Field(
+        description="The direction of the gravitational acceleration vector. "
+        + "This vector will be normalized automatically."
+    )
+    magnitude: AccelerationType = pd.Field(
+        description="The magnitude of the gravitational acceleration. "
+        + "For Earth's surface gravity, use 9.81 m/s²."
+    )
+    private_attribute_id: str = pd.Field(default_factory=generate_uuid, frozen=True)
+
+
 VolumeModelTypes = Union[
     Fluid,
     Solid,
@@ -1492,4 +1540,5 @@ def _validate_volumetric_heat_source_for_liquid(
     BETDisk,
     Rotation,
     PorousMedium,
+    Gravity,
 ]
diff --git a/flow360/component/simulation/simulation_params.py b/flow360/component/simulation/simulation_params.py
index d07ba044e..a6ba0b15a 100644
--- a/flow360/component/simulation/simulation_params.py
+++ b/flow360/component/simulation/simulation_params.py
@@ -117,6 +117,7 @@
     _check_duplicate_entities_in_models,
     _check_duplicate_isosurface_names,
     _check_duplicate_surface_usage,
+    _check_gravity_model_conflicts,
     _check_hybrid_model_to_use_zonal_enforcement,
     _check_low_mach_preconditioner_output,
     _check_numerical_dissipation_factor_output,
@@ -587,6 +588,11 @@ def check_duplicate_entities_in_models(self, param_info: ParamsValidationInfo):
         """Only allow each Surface/Volume entity to appear once in the Surface/Volume model"""
         return _check_duplicate_entities_in_models(self, param_info)
 
+    @pd.model_validator(mode="after")
+    def check_gravity_model_conflicts(self):
+        """Disallow multiple Gravity models when one applies to all zones"""
+        return _check_gravity_model_conflicts(self)
+
     @contextual_model_validator(mode="after")
     def check_unique_selector_names(self):
         """Ensure all EntitySelector names are unique"""
diff --git a/flow360/component/simulation/translator/solver_translator.py b/flow360/component/simulation/translator/solver_translator.py
index 132c285bd..6d947986c 100644
--- a/flow360/component/simulation/translator/solver_translator.py
+++ b/flow360/component/simulation/translator/solver_translator.py
@@ -49,6 +49,7 @@
     AngularVelocity,
     BETDisk,
     Fluid,
+    Gravity,
     NavierStokesInitialCondition,
     NavierStokesModifiedRestartSolution,
     PorousMedium,
@@ -1335,6 +1336,43 @@ def porous_media_translator(model: PorousMedium):
     return porous_medium
 
 
+def gravity_entity_info_serializer(volume):
+    """Gravity entity serializer"""
+    if volume is None:
+        return {"zoneType": "global"}
+    return {
+        "zoneType": "mesh",
+        "zoneName": volume.full_name,
+    }
+
+
+def gravity_translator(model: Gravity, params):
+    """Gravity translator - converts dimensional gravity to non-dimensional form.
+
+    Non-dimensionalization: g* = g * L_ref / a_∞²
+    where L_ref is the reference length (mesh unit) and a_∞ is the speed of sound.
+    """
+    # Get the magnitude value in m/s²
+    magnitude = model.magnitude.to("m/s**2").value.item()
+
+    # Get direction (already normalized by validator)
+    direction = model.direction
+
+    # Compute the dimensional gravity vector
+    gravity_vector_dimensional = [d * magnitude for d in direction]
+
+    # Compute non-dimensionalization factor: L_ref / a_∞²
+    base_length = params.base_length.to("m").value
+    base_velocity = params.base_velocity.to("m/s").value  # speed of sound
+
+    nondim_factor = base_length / (base_velocity**2)
+
+    # Non-dimensionalize the gravity vector
+    gravity_vector_nondim = [g * nondim_factor for g in gravity_vector_dimensional]
+
+    return {"gravityVector": gravity_vector_nondim}
+
+
 def bet_disk_entity_info_serializer(volume):
     """BET disk entity serializer"""
     v = convert_tuples_to_lists(remove_units_in_dict(dump_dict(volume)))
@@ -2351,6 +2389,23 @@ def get_solver_json(
             entity_injection_func=porous_media_entity_info_serializer,
         )
 
+    ##:: Step 9b: Get gravity
+    if has_instance_in_list(input_params.models, Gravity):
+        gravity_list = []
+        for gravity_model in get_all_entries_of_type(input_params.models, Gravity):
+            gravity_dict = gravity_translator(gravity_model, input_params)
+            if gravity_model.entities is None:
+                # Apply to all zones (global)
+                gravity_dict.update(gravity_entity_info_serializer(None))
+                gravity_list.append(gravity_dict)
+            else:
+                # Apply to specified zones
+                for entity in gravity_model.entities.stored_entities:
+                    entity_gravity_dict = gravity_dict.copy()
+                    entity_gravity_dict.update(gravity_entity_info_serializer(entity))
+                    gravity_list.append(entity_gravity_dict)
+        translated["gravity"] = gravity_list
+
     ##:: Step 10: Get heat transfer zones
     solid_zone_boundaries = set()
     if has_instance_in_list(input_params.models, Solid):
diff --git a/flow360/component/simulation/unit_system.py b/flow360/component/simulation/unit_system.py
index 5b6c9acc8..9ebf1983d 100644
--- a/flow360/component/simulation/unit_system.py
+++ b/flow360/component/simulation/unit_system.py
@@ -31,6 +31,7 @@
 udim.viscosity = udim.pressure * udim.time
 udim.kinematic_viscosity = udim.length * udim.length / udim.time
 udim.angular_velocity = udim.angle / udim.time
+udim.acceleration = udim.length / udim.time**2
 udim.heat_flux = udim.mass / udim.time**3
 udim.moment = udim.force * udim.length
 udim.heat_source = udim.mass / udim.time**3 / udim.length
@@ -918,6 +919,16 @@ class _VelocityType(_DimensionedType):
 VelocityType = Annotated[_VelocityType, PlainSerializer(_dimensioned_type_serializer)]
 
 
+class _AccelerationType(_DimensionedType):
+    """:class: AccelerationType"""
+
+    dim = udim.acceleration
+    dim_name = "acceleration"
+
+
+AccelerationType = Annotated[_AccelerationType, PlainSerializer(_dimensioned_type_serializer)]
+
+
 class _AreaType(_DimensionedType):
     """:class: AreaType"""
 
@@ -1379,6 +1390,13 @@ class Flow360VelocityUnit(_Flow360BaseUnit):
     unit_name = "flow360_velocity_unit"
 
 
+class Flow360AccelerationUnit(_Flow360BaseUnit):
+    """:class: Flow360AccelerationUnit"""
+
+    dimension_type = AccelerationType
+    unit_name = "flow360_acceleration_unit"
+
+
 class Flow360AreaUnit(_Flow360BaseUnit):
     """:class: Flow360AreaUnit"""
 
@@ -1548,6 +1566,7 @@ class BaseSystemType(Enum):
     "time",
     "temperature",
     "velocity",
+    "acceleration",
     "area",
     "force",
     "pressure",
@@ -1581,6 +1600,7 @@ class UnitSystem(pd.BaseModel):
     time: TimeType = pd.Field()
     temperature: AbsoluteTemperatureType = pd.Field()
     velocity: VelocityType = pd.Field()
+    acceleration: AccelerationType = pd.Field()
     area: AreaType = pd.Field()
     force: ForceType = pd.Field()
     pressure: PressureType = pd.Field()
@@ -1739,6 +1759,7 @@ def __exit__(self, exc_type, exc_val, exc_tb):
 flow360_time_unit = Flow360TimeUnit()
 flow360_temperature_unit = Flow360TemperatureUnit()
 flow360_velocity_unit = Flow360VelocityUnit()
+flow360_acceleration_unit = Flow360AccelerationUnit()
 flow360_area_unit = Flow360AreaUnit()
 flow360_force_unit = Flow360ForceUnit()
 flow360_pressure_unit = Flow360PressureUnit()
@@ -1766,6 +1787,7 @@ def __exit__(self, exc_type, exc_val, exc_tb):
     flow360_time_unit,
     flow360_temperature_unit,
     flow360_velocity_unit,
+    flow360_acceleration_unit,
     flow360_area_unit,
     flow360_force_unit,
     flow360_pressure_unit,
@@ -1807,6 +1829,9 @@ class Flow360ConversionUnitSystem(pd.BaseModel):
     base_velocity: float = pd.Field(
         np.inf, json_schema_extra={"target_dimension": Flow360VelocityUnit}
     )
+    base_acceleration: float = pd.Field(
+        np.inf, json_schema_extra={"target_dimension": Flow360AccelerationUnit}
+    )
     base_area: float = pd.Field(np.inf, json_schema_extra={"target_dimension": Flow360AreaUnit})
     base_force: float = pd.Field(np.inf, json_schema_extra={"target_dimension": Flow360ForceUnit})
     base_density: float = pd.Field(
@@ -1886,6 +1911,7 @@ def __init__(self):
         )
 
         conversion_system["velocity"] = "flow360_velocity_unit"
+        conversion_system["acceleration"] = "flow360_acceleration_unit"
         conversion_system["area"] = "flow360_area_unit"
         conversion_system["force"] = "flow360_force_unit"
         conversion_system["density"] = "flow360_density_unit"
@@ -1936,6 +1962,7 @@ class _PredefinedUnitSystem(UnitSystem):
     time: TimeType = pd.Field(exclude=True)
     temperature: AbsoluteTemperatureType = pd.Field(exclude=True)
     velocity: VelocityType = pd.Field(exclude=True)
+    acceleration: AccelerationType = pd.Field(exclude=True)
     area: AreaType = pd.Field(exclude=True)
     force: ForceType = pd.Field(exclude=True)
     pressure: PressureType = pd.Field(exclude=True)
diff --git a/flow360/component/simulation/validation/validation_simulation_params.py b/flow360/component/simulation/validation/validation_simulation_params.py
index 2d422f37e..8c20b5e4b 100644
--- a/flow360/component/simulation/validation/validation_simulation_params.py
+++ b/flow360/component/simulation/validation/validation_simulation_params.py
@@ -31,6 +31,7 @@
 from flow360.component.simulation.models.volume_models import (
     ActuatorDisk,
     Fluid,
+    Gravity,
     Rotation,
     Solid,
 )
@@ -115,7 +116,7 @@ def _check_duplicate_entities_in_models(params, param_info: ParamsValidationInfo
     usage = EntityUsageMap()
 
     for model in models:
-        if hasattr(model, "entities"):
+        if hasattr(model, "entities") and model.entities is not None:
             expanded_entities = param_info.expand_entity_list(model.entities)
             # seen_entity_hashes: set[str] = set()
             for entity in expanded_entities:
@@ -145,6 +146,39 @@ def _check_duplicate_entities_in_models(params, param_info: ParamsValidationInfo
     return params
 
 
+def _check_gravity_model_conflicts(params):
+    """Check that at most one Gravity model applies to all zones (entities=None).
+
+    When a Gravity model has entities=None it applies to every fluid zone,
+    so having more than one such model or mixing a global one with
+    zone-specific ones would create conflicting gravity definitions.
+    """
+    if not params.models:
+        return params
+
+    gravity_models = [m for m in params.models if isinstance(m, Gravity)]
+    if len(gravity_models) <= 1:
+        return params
+
+    global_gravity = [m for m in gravity_models if m.entities is None]
+
+    if len(global_gravity) > 1:
+        raise ValueError(
+            "Multiple Gravity models with unspecified entities (applying to all zones) "
+            "are not allowed. Please specify explicit entities for each Gravity model "
+            "or use a single Gravity model for all zones."
+        )
+
+    if len(global_gravity) == 1:
+        raise ValueError(
+            "A Gravity model that applies to all zones (entities not specified) "
+            "cannot coexist with other Gravity models. Please specify explicit "
+            "entities for each Gravity model or use a single Gravity model for all zones."
+        )
+
+    return params
+
+
 def _check_low_mach_preconditioner_output(v):
     models = v.models
 
diff --git a/tests/simulation/params/test_gravity.py b/tests/simulation/params/test_gravity.py
new file mode 100644
index 000000000..881f99b66
--- /dev/null
+++ b/tests/simulation/params/test_gravity.py
@@ -0,0 +1,300 @@
+import math
+import re
+
+import pytest
+
+import flow360.component.simulation.units as u
+from flow360.component.simulation.models.volume_models import Gravity
+from flow360.component.simulation.primitives import GenericVolume
+from flow360.component.simulation.simulation_params import SimulationParams
+from flow360.component.simulation.translator.solver_translator import (
+    gravity_entity_info_serializer,
+    gravity_translator,
+)
+from flow360.component.simulation.unit_system import SI_unit_system
+from flow360.component.simulation.validation.validation_context import (
+    ParamsValidationInfo,
+    ValidationContext,
+)
+
+
+class MockSimulationParams:
+    """Mock class for testing gravity_translator"""
+
+    def __init__(self, base_length_m=1.0, base_velocity_ms=340.0):
+        self._base_length = base_length_m * u.m
+        self._base_velocity = base_velocity_ms * u.m / u.s
+
+    @property
+    def base_length(self):
+        return self._base_length
+
+    @property
+    def base_velocity(self):
+        return self._base_velocity
+
+
+def test_gravity_creation_basic():
+    """Test basic Gravity model creation with default units."""
+    gravity = Gravity(
+        direction=(0, 0, -1),
+        magnitude=9.81 * u.m / u.s**2,
+    )
+    assert gravity.type == "Gravity"
+    assert gravity.name == "Gravity"
+    assert gravity.entities is None  # Default to all zones
+    # Direction should be normalized (already is in this case)
+    assert tuple(gravity.direction) == (0.0, 0.0, -1.0)
+
+
+def test_gravity_direction_normalization():
+    """Test that direction is normalized automatically."""
+    gravity = Gravity(
+        direction=(0, 3, -4),  # magnitude is 5, will be normalized
+        magnitude=9.81 * u.m / u.s**2,
+    )
+    # Expected: (0, 0.6, -0.8)
+    assert math.isclose(gravity.direction[0], 0.0, abs_tol=1e-10)
+    assert math.isclose(gravity.direction[1], 0.6, rel_tol=1e-10)
+    assert math.isclose(gravity.direction[2], -0.8, rel_tol=1e-10)
+
+
+def test_gravity_zero_direction_raises():
+    """Test that zero direction vector raises an error."""
+    with pytest.raises(ValueError, match=re.escape("Axis cannot be (0, 0, 0)")):
+        Gravity(
+            direction=(0, 0, 0),
+            magnitude=9.81 * u.m / u.s**2,
+        )
+
+
+def test_gravity_with_entities():
+    """Test Gravity model with specific entities."""
+    mock_context = ValidationContext(
+        levels=None, info=ParamsValidationInfo(param_as_dict={}, referenced_expressions=[])
+    )
+
+    with mock_context:
+        gravity = Gravity(
+            entities=[GenericVolume(name="fluid_zone")],
+            direction=(0, 0, -1),
+            magnitude=9.81 * u.m / u.s**2,
+        )
+        assert gravity.entities is not None
+        assert len(gravity.entities.stored_entities) == 1
+
+
+def test_gravity_different_units():
+    """Test Gravity model with different acceleration units."""
+    # Using ft/s^2
+    gravity_imperial = Gravity(
+        direction=(0, 0, -1),
+        magnitude=32.174 * u.ft / u.s**2,  # ~9.8 m/s^2
+    )
+    # Verify magnitude is stored correctly
+    assert gravity_imperial.magnitude.to("m/s**2").value > 9.7
+    assert gravity_imperial.magnitude.to("m/s**2").value < 9.9
+
+
+def test_gravity_custom_name():
+    """Test Gravity model with custom name."""
+    gravity = Gravity(
+        name="Earth Gravity",
+        direction=(0, 0, -1),
+        magnitude=9.81 * u.m / u.s**2,
+    )
+    assert gravity.name == "Earth Gravity"
+
+
+def test_gravity_arbitrary_direction():
+    """Test Gravity model with arbitrary direction."""
+    gravity = Gravity(
+        direction=(1, 1, 1),
+        magnitude=9.81 * u.m / u.s**2,
+    )
+    # Should be normalized to unit vector
+    norm = math.sqrt(sum(d**2 for d in gravity.direction))
+    assert math.isclose(norm, 1.0, rel_tol=1e-10)
+
+
+def test_gravity_translator_nondimensionalization():
+    """Test that gravity_translator correctly non-dimensionalizes the gravity vector.
+
+    Non-dimensionalization: g* = g * L_ref / a_∞²
+    """
+    gravity = Gravity(
+        direction=(0, 0, -1),
+        magnitude=9.81 * u.m / u.s**2,
+    )
+
+    # With L_ref = 1 m and a_∞ = 340 m/s (speed of sound at sea level)
+    mock_params = MockSimulationParams(base_length_m=1.0, base_velocity_ms=340.0)
+
+    result = gravity_translator(gravity, mock_params)
+
+    # Expected: g* = 9.81 * 1 / 340^2 = 9.81 / 115600 ≈ 8.49e-5
+    expected_nondim = 9.81 / (340.0**2)
+
+    assert "gravityVector" in result
+    assert len(result["gravityVector"]) == 3
+    assert math.isclose(result["gravityVector"][0], 0.0, abs_tol=1e-15)
+    assert math.isclose(result["gravityVector"][1], 0.0, abs_tol=1e-15)
+    assert math.isclose(result["gravityVector"][2], -expected_nondim, rel_tol=1e-5)
+
+
+def test_gravity_translator_direction():
+    """Test that gravity_translator preserves direction correctly."""
+    gravity = Gravity(
+        direction=(1, 0, 0),  # Gravity in +x direction
+        magnitude=10.0 * u.m / u.s**2,
+    )
+
+    mock_params = MockSimulationParams(base_length_m=1.0, base_velocity_ms=100.0)
+    result = gravity_translator(gravity, mock_params)
+
+    # Expected: g* = 10 * 1 / 100^2 = 10 / 10000 = 0.001
+    expected_nondim = 10.0 / (100.0**2)
+
+    assert math.isclose(result["gravityVector"][0], expected_nondim, rel_tol=1e-5)
+    assert math.isclose(result["gravityVector"][1], 0.0, abs_tol=1e-15)
+    assert math.isclose(result["gravityVector"][2], 0.0, abs_tol=1e-15)
+
+
+def test_gravity_entity_info_serializer_global():
+    """Test gravity entity serializer for global gravity."""
+    result = gravity_entity_info_serializer(None)
+    assert result == {"zoneType": "global"}
+
+
+def test_gravity_entity_info_serializer_zone():
+    """Test gravity entity serializer for zone-specific gravity."""
+    mock_context = ValidationContext(
+        levels=None, info=ParamsValidationInfo(param_as_dict={}, referenced_expressions=[])
+    )
+
+    with mock_context:
+        volume = GenericVolume(name="test_zone")
+        result = gravity_entity_info_serializer(volume)
+
+        assert result["zoneType"] == "mesh"
+        assert result["zoneName"] == volume.full_name
+
+
+def test_gravity_very_small_direction():
+    """Test that very small but non-zero direction is normalized correctly."""
+    gravity = Gravity(
+        direction=(1e-10, 1e-10, 1e-10),
+        magnitude=9.81 * u.m / u.s**2,
+    )
+    norm = math.sqrt(sum(d**2 for d in gravity.direction))
+    assert math.isclose(norm, 1.0, rel_tol=1e-10)
+
+
+def test_gravity_negative_direction_components():
+    """Test Gravity model with all negative direction components."""
+    gravity = Gravity(
+        direction=(-1, -1, -1),
+        magnitude=9.81 * u.m / u.s**2,
+    )
+    # All components should be negative and equal after normalization
+    expected = -1.0 / math.sqrt(3)
+    assert math.isclose(gravity.direction[0], expected, rel_tol=1e-10)
+    assert math.isclose(gravity.direction[1], expected, rel_tol=1e-10)
+    assert math.isclose(gravity.direction[2], expected, rel_tol=1e-10)
+
+
+def test_gravity_large_magnitude():
+    """Test Gravity model with large magnitude (e.g., Jupiter-like)."""
+    gravity = Gravity(
+        direction=(0, 0, -1),
+        magnitude=24.79 * u.m / u.s**2,  # Jupiter's surface gravity
+    )
+    assert gravity.magnitude.to("m/s**2").value > 24.0
+    assert gravity.magnitude.to("m/s**2").value < 25.0
+
+
+def test_gravity_small_magnitude():
+    """Test Gravity model with small magnitude (e.g., Moon-like)."""
+    gravity = Gravity(
+        direction=(0, 0, -1),
+        magnitude=1.62 * u.m / u.s**2,  # Moon's surface gravity
+    )
+    assert gravity.magnitude.to("m/s**2").value > 1.6
+    assert gravity.magnitude.to("m/s**2").value < 1.7
+
+
+def test_single_global_gravity_is_valid():
+    """A single Gravity model with entities=None (global) should be accepted."""
+    gravity = Gravity(
+        direction=(0, 0, -1),
+        magnitude=9.81 * u.m / u.s**2,
+    )
+    with SI_unit_system:
+        params = SimulationParams(models=[gravity])
+    assert params
+
+
+def test_multiple_global_gravity_raises():
+    """Two Gravity models with entities=None should raise a conflict error."""
+    gravity1 = Gravity(
+        direction=(0, 0, -1),
+        magnitude=9.81 * u.m / u.s**2,
+    )
+    gravity2 = Gravity(
+        direction=(1, 0, 0),
+        magnitude=5.0 * u.m / u.s**2,
+    )
+    with SI_unit_system, pytest.raises(
+        ValueError,
+        match=re.escape(
+            "Multiple Gravity models with unspecified entities (applying to all zones) "
+            "are not allowed."
+        ),
+    ):
+        SimulationParams(models=[gravity1, gravity2])
+
+
+def test_global_gravity_with_zone_specific_raises():
+    """A global Gravity (entities=None) mixed with zone-specific Gravity should raise."""
+    mock_context = ValidationContext(
+        levels=None, info=ParamsValidationInfo(param_as_dict={}, referenced_expressions=[])
+    )
+    global_gravity = Gravity(
+        direction=(0, 0, -1),
+        magnitude=9.81 * u.m / u.s**2,
+    )
+    with mock_context:
+        zone_gravity = Gravity(
+            entities=[GenericVolume(name="zone1")],
+            direction=(1, 0, 0),
+            magnitude=5.0 * u.m / u.s**2,
+        )
+    with SI_unit_system, pytest.raises(
+        ValueError,
+        match=re.escape(
+            "A Gravity model that applies to all zones (entities not specified) "
+            "cannot coexist with other Gravity models."
+        ),
+    ):
+        SimulationParams(models=[global_gravity, zone_gravity])
+
+
+def test_multiple_zone_specific_gravity_is_valid():
+    """Multiple Gravity models with distinct entities should be accepted."""
+    mock_context = ValidationContext(
+        levels=None, info=ParamsValidationInfo(param_as_dict={}, referenced_expressions=[])
+    )
+    with mock_context:
+        gravity1 = Gravity(
+            entities=[GenericVolume(name="zone1")],
+            direction=(0, 0, -1),
+            magnitude=9.81 * u.m / u.s**2,
+        )
+        gravity2 = Gravity(
+            entities=[GenericVolume(name="zone2")],
+            direction=(1, 0, 0),
+            magnitude=5.0 * u.m / u.s**2,
+        )
+    with SI_unit_system:
+        params = SimulationParams(models=[gravity1, gravity2])
+    assert params