diff --git a/workflow/scripts/generate_velocity_model_parameters.py b/workflow/scripts/generate_velocity_model_parameters.py
index 790240f..175984c 100644
--- a/workflow/scripts/generate_velocity_model_parameters.py
+++ b/workflow/scripts/generate_velocity_model_parameters.py
@@ -426,45 +426,57 @@ def pgv_target(
     )
 
 
-@cli.from_docstring(app)
-@log_utils.log_call()
 def generate_velocity_model_parameters(
-    realisation_ffp: Annotated[Path, typer.Argument()],
-) -> None:
-    """Generate velocity model parameters for a given realisation file.
-
-    This function reads the source and rupture propagation information and computes:
-
-    1. The size of the simulation domain,
-    2. The simulation duration.
-
-    Both of these values are written to the realisation using `VelocityModelParameters`.
+    source_config: SourceConfig,
+    velocity_model_parameters: VelocityModelParameters,
+    rupture_propagation: RupturePropagationConfig,
+    rakes: Rakes,
+    magnitudes: Magnitudes,
+) -> DomainParameters:
+    """Generate model domain extents, maximum depth, and simulation duration for a realisation.
 
     Parameters
     ----------
-    realisation_ffp : Path
-        The path to the realisation file from which to read configurations and to which
-        the generated velocity model parameters will be written.
+    source_config : SourceConfig
+        Source configuration containing a mapping of fault/source geometries.
+    velocity_model_parameters : VelocityModelParameters
+        Velocity-model configuration and controls. Used to compute
+        the PGV target for rrup bounding polygons.
+    rupture_propagation : RupturePropagationConfig
+        Rupture propagation settings.
+    rakes : Rakes
+        Mapping of fault identifiers to rake angles (degrees). Used in
+        the GMPE that estimates domain parameters.
+    magnitudes : Magnitudes
+        Mapping of fault identifiers to magnitudes (Mw). Used to compute the
+        total rupture magnitude and rrup bounding polygons.
 
     Returns
     -------
-    None
-        The function does not return any value. It writes the computed parameters to
-        the specified realisation file.
+    DomainParameters
+        The domain geometry for simulation.
+
+    See Also
+    --------
+    total_magnitude
+    pgv_target
+    get_max_depth
+    find_rrup_bounding_polygon
+    bounding_box.minimum_area_bounding_box_for_polygons_masked
+    estimate_simulation_duration
+
+    Notes
+    -----
+    - Fault polygons are buffered by 2,000 meters and must be fully included
+      in the final domain.
+    - Rrup bounding polygons represent areas within rupture distance of any
+      source-geometry corner; they may be included when they overlap land.
+    - The domain computation is masked by the New Zealand outline `utils.get_nz_outline_polygon`,
+      ensuring only land areas are considered for optional inclusion.
+    - Depth is computed from the total rupture magnitude and the initial fault's
+      bottom depth.
     """
-    metadata = RealisationMetadata.read_from_realisation(realisation_ffp)
-    source_config = SourceConfig.read_from_realisation(realisation_ffp)
-    velocity_model_parameters = (
-        VelocityModelParameters.read_from_realisation_or_defaults(
-            realisation_ffp, metadata.defaults_version
-        )
-    )
 
-    rupture_propagation = RupturePropagationConfig.read_from_realisation(
-        realisation_ffp
-    )
-    magnitudes = Magnitudes.read_from_realisation(realisation_ffp)
-    rakes = Rakes.read_from_realisation(realisation_ffp)
     rupture_magnitude = total_magnitude(np.array(list(magnitudes.magnitudes.values())))
     realisation_pgv_target = pgv_target(magnitudes, velocity_model_parameters)
 
@@ -483,7 +495,6 @@ def generate_velocity_model_parameters(
         shapely.buffer(fault.geometry, 2000)
         for fault in source_config.source_geometries.values()
     ]
-    rakes = Rakes.read_from_realisation(realisation_ffp)
     # This polygon includes all areas within rrup distance of any
     # corner in the source geometries.
     # These may be in the domain where they are over land.
@@ -516,10 +527,60 @@ def generate_velocity_model_parameters(
         velocity_model_parameters.s_wave_velocity,
     )
 
-    domain_parameters = DomainParameters(
+    return DomainParameters(
         domain=model_domain,
         depth=max_depth,
         duration=sim_duration,
     )
+
+
+@cli.from_docstring(app)
+@log_utils.log_call()
+def main(
+    realisation_ffp: Annotated[Path, typer.Argument()],
+) -> None:
+    """Generate velocity model parameters for a given realisation file.
+
+    This function reads the source and rupture propagation information and computes:
+
+    1. The size of the simulation domain,
+    2. The simulation duration.
+
+    Both of these values are written to the realisation using `DomainParameters`.
+
+    Parameters
+    ----------
+    realisation_ffp : Path
+        The path to the realisation file from which to read configurations and to which
+        the generated velocity model parameters will be written.
+
+    Returns
+    -------
+    None
+        The function does not return any value. It writes the computed parameters to
+        the specified realisation file.
+    """
+    metadata = RealisationMetadata.read_from_realisation(realisation_ffp)
+    source_config = SourceConfig.read_from_realisation(realisation_ffp)
+    velocity_model_parameters = (
+        VelocityModelParameters.read_from_realisation_or_defaults(
+            realisation_ffp, metadata.defaults_version
+        )
+    )
+
+    rupture_propagation = RupturePropagationConfig.read_from_realisation(
+        realisation_ffp
+    )
+    magnitudes = Magnitudes.read_from_realisation(realisation_ffp)
+    rakes = Rakes.read_from_realisation(realisation_ffp)
+
+    domain_parameters = generate_velocity_model_parameters(
+        source_config=source_config,
+        velocity_model_parameters=velocity_model_parameters,
+        rupture_propagation=rupture_propagation,
+        magnitudes=magnitudes,
+        rakes=rakes,
+    )
+
     domain_parameters.write_to_realisation(realisation_ffp)
     realisations.append_log_entry(realisation_ffp)