diff --git a/compass/ocean/mesh/cull.py b/compass/ocean/mesh/cull.py
index 452ecd97e6..c6e4c62261 100644
--- a/compass/ocean/mesh/cull.py
+++ b/compass/ocean/mesh/cull.py
@@ -305,6 +305,17 @@ def _cull_mesh_with_logging(logger, with_cavities, with_critical_passages,
                                 geojson_filename='land_coverage.geojson',
                                 mask_filename='land_mask.nc')
 
+    # Hook point: allow modification of land mask before culling
+    # This enables regional culling and other mask modifications
+    if os.path.exists('modify_land_mask_hook.py'):
+        logger.info('Running land mask modification hook...')
+        import importlib.util
+        spec = importlib.util.spec_from_file_location("hook", "modify_land_mask_hook.py")
+        hook = importlib.util.module_from_spec(spec)
+        spec.loader.exec_module(hook)
+        if hasattr(hook, 'modify_land_mask'):
+            hook.modify_land_mask(logger=logger)
+
     dsBaseMesh = xr.open_dataset('base_mesh.nc')
     dsLandMask = xr.open_dataset('land_mask.nc')
 
diff --git a/compass/ocean/tests/global_ocean/mesh/__init__.py b/compass/ocean/tests/global_ocean/mesh/__init__.py
index 764848c7de..f5be52105b 100644
--- a/compass/ocean/tests/global_ocean/mesh/__init__.py
+++ b/compass/ocean/tests/global_ocean/mesh/__init__.py
@@ -12,6 +12,8 @@
 from compass.ocean.tests.global_ocean.mesh.fris02to60 import FRIS02to60BaseMesh
 from compass.ocean.tests.global_ocean.mesh.fris04to60 import FRIS04to60BaseMesh
 from compass.ocean.tests.global_ocean.mesh.fris08to60 import FRIS08to60BaseMesh
+from compass.ocean.tests.global_ocean.mesh.thwaites01to60 import Thwaites01to60BaseMesh
+from compass.ocean.tests.global_ocean.mesh.thwaites01to60.cull_mesh import ThwaitesCullMeshStep
 from compass.ocean.tests.global_ocean.mesh.kuroshio import KuroshioBaseMesh
 from compass.ocean.tests.global_ocean.mesh.qu import (
     IcosMeshFromConfigStep,
@@ -138,6 +140,8 @@ def __init__(self, test_group, mesh_name,  # noqa: C901
             base_mesh_step = FRIS04to60BaseMesh(self, name=name, subdir=subdir)
         elif mesh_name in ['FRIS08to60', 'FRISwISC08to60']:
             base_mesh_step = FRIS08to60BaseMesh(self, name=name, subdir=subdir)
+        elif mesh_name in ['Thwaites01to60', 'ThwaitesWISC01to60']:
+            base_mesh_step = Thwaites01to60BaseMesh(self, name=name, subdir=subdir)
         elif mesh_name.startswith('Kuroshio'):
             base_mesh_step = KuroshioBaseMesh(self, name=name, subdir=subdir)
         elif mesh_name in ['WC14', 'WCwISC14']:
@@ -191,10 +195,17 @@ def __init__(self, test_group, mesh_name,  # noqa: C901
 
             self.add_step(smoothed_topo)
 
-        self.add_step(CullMeshStep(
-            test_case=self, base_mesh_step=base_mesh_step,
-            with_ice_shelf_cavities=self.with_ice_shelf_cavities,
-            unsmoothed_topo=unsmoothed_topo, smoothed_topo=smoothed_topo))
+        # Use custom cull step for Thwaites meshes (supports regional culling)
+        if mesh_name.startswith('Thwaites'):
+            self.add_step(ThwaitesCullMeshStep(
+                test_case=self, base_mesh_step=base_mesh_step,
+                with_ice_shelf_cavities=self.with_ice_shelf_cavities,
+                unsmoothed_topo=unsmoothed_topo, smoothed_topo=smoothed_topo))
+        else:
+            self.add_step(CullMeshStep(
+                test_case=self, base_mesh_step=base_mesh_step,
+                with_ice_shelf_cavities=self.with_ice_shelf_cavities,
+                unsmoothed_topo=unsmoothed_topo, smoothed_topo=smoothed_topo))
 
     def configure(self, config=None):
         """
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/README.md b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/README.md
new file mode 100644
index 0000000000..75c4b7c707
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/README.md
@@ -0,0 +1,141 @@
+# Thwaites01to60 Mesh
+
+Variable-resolution global ocean mesh with 200m refinement at Thwaites grounding zone.
+
+## Current Status
+
+### ✅ Implemented:
+1. **Base mesh generation** (`__init__.py`)
+   - `Thwaites01to60BaseMesh` class
+   - **Grounding line extraction from BedMachine** (automatic, no geojson needed)
+   - Calculates flotation criterion: thickness > -bed × (ρ_ocean/ρ_ice)
+   - Identifies grounding line as boundary between grounded/floating ice
+   - Creates signed distance field for mesh refinement
+   - Falls back to Gaussian approximation if BedMachine unavailable
+   - Configurable resolutions: 200m GZ, 1km cavity, 3km shelf, 8km far-field
+
+2. **Configuration** (`thwaites01to60.cfg`)
+   - Resolution parameters
+   - Regional domain settings (enabled by default for Amundsen sector)
+   - Thin film parameters (disabled by default)
+
+3. **Regional culling** ✅ **FULLY INTEGRATED**
+   - `ThwaitesCullMeshStep` custom cull step (`cull_mesh.py`)
+   - Hook-based integration into standard culling workflow
+   - Modifies land mask to mark cells outside domain as land
+   - Supports geojson polygon or lat-lon bounding box
+   - **Working and ready to test**
+
+4. **Registration**
+   - Mesh registered in `compass/ocean/tests/global_ocean/mesh/__init__.py`
+   - Custom CullMeshStep wired in for Thwaites meshes
+   - Ready for `compass list` and `compass setup`
+
+### 🚧 TODO:
+
+#### High Priority:
+1. ~~**Integrate regional culling into CullMeshStep**~~ ✅ **DONE**
+   - ✅ Added hook point in base `cull.py`
+   - ✅ Created `ThwaitesCullMeshStep` with hook file generation
+   - ✅ Wired into Mesh test case for Thwaites meshes
+   - ✅ Ready for testing
+
+2. ~~**Replace Gaussian refinement with grounding-line-based refinement**~~ ✅ **DONE**
+   - ✅ Extracts grounding line directly from BedMachine
+   - ✅ Calculates flotation criterion automatically
+   - ✅ Uses `signed_distance_from_geojson` for accurate GZ band
+   - ✅ No manual geojson file required
+   - ⚠️ Requires BedMachine file in bathymetry database
+
+#### Medium Priority:
+3. **Add thin film support**
+   - Modify land mask to keep cells beneath grounded ice
+   - Based on height above flotation from BedMachine
+   - Requires integration with `remap_topography` step
+
+4. **Testing**
+   - Set up and run basic mesh generation
+   - Validate cell widths and resolution
+   - Test regional culling with lat-lon bounds
+   - Test with/without thin film
+
+#### Low Priority:
+5. **Clean up copied FRIS files**
+   - Remove unused geojson files (atlantic.geojson, fris_*.geojson, etc.)
+   - Remove or adapt FRIS-specific namelists if needed
+
+## Usage
+
+### Basic mesh (no regional culling):
+```bash
+compass list | grep Thwaites
+compass setup -t global_ocean/mesh/Thwaites01to60 -w $WORK
+compass run $WORK
+```
+
+### With regional domain (enabled by default):
+Regional culling is **enabled by default** for Amundsen sector (76-73°S, 116-98°W).
+
+To change the domain, edit `thwaites01to60.cfg`:
+```ini
+[thwaites01to60]
+# Option A: lat-lon bounds (currently enabled)
+lat_min = -76.0
+lat_max = -73.0
+lon_min = -116.0
+lon_max = -98.0
+
+# Option B: geojson polygon
+# regional_domain_geojson = amundsen_domain.geojson
+```
+
+To disable regional culling (full global mesh), comment out all domain options.
+
+## Implementation Notes
+
+### Mesh Refinement:
+The implementation now extracts the grounding line directly from BedMachine:
+
+1. **Load BedMachine data** (ice thickness, bed elevation)
+2. **Calculate flotation criterion**:
+   ```python
+   # Ice is grounded when: thickness > -bed × (ρ_ocean/ρ_ice)
+   thickness_flotation = -bed * (1028.0 / 918.0)  # where bed < 0
+   grounded = thickness > thickness_flotation
+   ```
+3. **Identify grounding line**: Boundary between grounded and floating ice
+4. **Create feature collection**: Convert GL cells to geojson polygon
+5. **Calculate signed distance**: Use `signed_distance_from_geojson`
+6. **Apply refinement**:
+   - GZ band: ±15 km from GL → 200m resolution
+   - Cavity: 100 km into cavity → 1 km resolution  
+   - Shelf: 200 km into cavity → 3 km resolution
+   - Far-field: everywhere else → 8 km resolution
+
+**Fallback**: If BedMachine is unavailable, uses Gaussian approximation around Thwaites location (75°S, 106°W).
+
+### Regional Culling:
+**Fully implemented and integrated!**
+
+The `ThwaitesCullMeshStep` (in `cull_mesh.py`) creates a hook file that:
+1. Is executed by the modified `cull.py` after land mask creation
+2. Loads base mesh and land mask
+3. Determines cells inside/outside regional domain
+4. Marks outside cells as "land"
+5. Standard culling then removes these cells
+
+The hook runs between lines 307 and 308 of `compass/ocean/mesh/cull.py`,
+after `land_mask.nc` is created but before culling begins.
+
+### Thin Film:
+Will require:
+1. BedMachine topography (from `remap_topography` step)
+2. Compute height above flotation for grounded ice
+3. Cells with HAF < threshold (e.g., 30m) kept as ocean
+4. Add `thinFilmMask` field for diagnostics
+
+## References
+
+- Design document: `/Users/trhille/Documents/Antarctica/thwaites_grounding_zone_intrusion/Grounding_zone_intrusion_simulation_design.pdf`
+- Implementation plan: `/Users/trhille/Documents/Antarctica/thwaites_grounding_zone_intrusion/QUICK_START.md`
+- FRIS mesh (template): `compass/ocean/tests/global_ocean/mesh/fris01to60/`
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/__init__.py b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/__init__.py
new file mode 100644
index 0000000000..040bd62f26
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/__init__.py
@@ -0,0 +1,308 @@
+import mpas_tools.mesh.creation.mesh_definition_tools as mdt
+import numpy as np
+import xarray as xr
+from geometric_features import FeatureCollection, GeometricFeatures
+from mpas_tools.cime.constants import constants
+from mpas_tools.mesh.creation.signed_distance import (
+    signed_distance_from_geojson,
+)
+from scipy import ndimage
+
+from compass.mesh import QuasiUniformSphericalMeshStep
+
+
+class Thwaites01to60BaseMesh(QuasiUniformSphericalMeshStep):
+    """
+    A step for creating Thwaites 200m-to-60km variable resolution mesh
+
+    Attributes
+    ----------
+    cell_width : numpy.ndarray
+        m x n array of cell width in km
+
+    x, y, z : numpy.ndarray
+        m x n arrays defining the sphere
+    """
+
+    def setup(self):
+        """
+        Add geojson files as inputs
+        """
+        # TODO: Add thwaites_grounding_line.geojson when available
+        # self.add_input_file(
+        #     filename='thwaites_grounding_line.geojson',
+        #     package=self.__module__)
+
+        super().setup()
+
+    def build_cell_width_lat_lon(self):
+        """
+        Create cell width array for this mesh on a regular latitude-longitude grid
+
+        Returns
+        -------
+        cellWidth : numpy.array
+            m x n array of cell width in km
+
+        lon : numpy.array
+            longitude in degrees (length n and between -180 and 180)
+
+        lat : numpy.array
+            latitude in degrees (length m and between -90 and 90)
+        """
+        # Get config parameters
+        config = self.config
+        section = config['thwaites01to60']
+
+        res_gz = section.getfloat('res_gz')
+        res_cavity = section.getfloat('res_cavity')
+        res_shelf = section.getfloat('res_shelf')
+        res_far = section.getfloat('res_far')
+        gz_band_halfwidth = section.getfloat('gz_band_halfwidth')
+
+        print('\nCreating Thwaites01to60 mesh with:')
+        print(f'  GZ resolution: {res_gz} km')
+        print(f'  Cavity resolution: {res_cavity} km')
+        print(f'  Shelf resolution: {res_shelf} km')
+        print(f'  Far-field resolution: {res_far} km')
+        print(f'  GZ band halfwidth: {gz_band_halfwidth} km')
+
+        # Create lat-lon grid for cellWidth
+        dlon = 0.1
+        dlat = dlon
+        earth_radius = constants['SHR_CONST_REARTH']
+        nlon = int(360. / dlon) + 1
+        nlat = int(180. / dlat) + 1
+        lon = np.linspace(-180., 180., nlon)
+        lat = np.linspace(-90., 90., nlat)
+
+        # Start with far-field resolution everywhere
+        cellWidth = res_far * np.ones((nlat, nlon))
+
+        # Extract and use grounding line from BedMachine
+        print('\n  Loading BedMachine to extract grounding line...')
+        gz_geojson = self._extract_grounding_line_geojson()
+
+        if gz_geojson is not None:
+            print('  Calculating signed distance from grounding line...')
+            # Calculate signed distance from grounding line
+            lon_grid, lat_grid = np.meshgrid(lon, lat)
+            gz_signed_distance = signed_distance_from_geojson(
+                gz_geojson, lon_grid, lat_grid, earth_radius,
+                max_length=0.25)
+
+            # Convert from meters to km
+            gz_signed_distance_km = gz_signed_distance / 1000.0
+
+            # Apply refinement in grounding zone band (±gz_band_halfwidth km)
+            # Use smooth tanh transition
+            transition_width = 5.0  # km
+
+            # GZ band refinement
+            print(f'  Applying GZ refinement (±{gz_band_halfwidth} km band)...')
+            gz_mask = 0.5 * (1 + np.tanh(
+                (np.abs(gz_signed_distance_km) - gz_band_halfwidth) / transition_width))
+            cellWidth = res_gz * (1 - gz_mask) + cellWidth * gz_mask
+
+            # Cavity refinement (ocean side, signed distance < 0)
+            # Refine 100km into cavity from grounding line
+            cavity_extent = 100.0  # km
+            print(f'  Applying cavity refinement ({cavity_extent} km from GL)...')
+            cavity_mask = 0.5 * (1 + np.tanh(
+                (-gz_signed_distance_km - cavity_extent) / transition_width))
+            # Only apply on ocean side
+            cavity_mask = np.where(gz_signed_distance_km < 0, cavity_mask, 1.0)
+            cellWidth = np.minimum(cellWidth,
+                                  res_cavity * (1 - cavity_mask) + cellWidth * cavity_mask)
+
+            # Ice shelf refinement (farther into cavity)
+            shelf_extent = 200.0  # km
+            print(f'  Applying shelf refinement ({shelf_extent} km from GL)...')
+            shelf_mask = 0.5 * (1 + np.tanh(
+                (-gz_signed_distance_km - shelf_extent) / transition_width))
+            shelf_mask = np.where(gz_signed_distance_km < 0, shelf_mask, 1.0)
+            cellWidth = np.minimum(cellWidth,
+                                  res_shelf * (1 - shelf_mask) + cellWidth * shelf_mask)
+
+        else:
+            print('  WARNING: Could not extract grounding line from BedMachine')
+            print('  Falling back to Gaussian approximation around Thwaites')
+            # Fall back to simple Gaussian refinement
+            cellWidth = self._apply_gaussian_refinement(
+                cellWidth, lon, lat, res_gz, res_cavity, res_shelf)
+
+        print(f'  CellWidth range: {cellWidth.min():.2f} - {cellWidth.max():.2f} km')
+
+        return cellWidth, lon, lat
+
+    def _extract_grounding_line_geojson(self):
+        """
+        Extract grounding line from BedMachine and create a geojson
+
+        Returns
+        -------
+        fc : geometric_features.FeatureCollection or None
+            Feature collection containing grounding line polygon(s)
+        """
+        try:
+            # Try to find BedMachine in the database
+            # The file should be specified in the config
+            config = self.config
+            bedmachine_file = None
+
+            # Check if we can find BedMachine file from compass database
+            # For now, try a few standard locations
+            import os
+            possible_paths = [
+                '/global/cfs/cdirs/e3sm/mpas_standalonedata/mpas-ocean/bathymetry_database/BedMachineAntarctica-v3.nc',
+                '/usr/projects/climate/SHARED_CLIMATE/mpas_standalonedata/mpas-ocean/bathymetry_database/BedMachineAntarctica-v3.nc',
+            ]
+
+            for path in possible_paths:
+                if os.path.exists(path):
+                    bedmachine_file = path
+                    break
+
+            if bedmachine_file is None:
+                print('    Could not find BedMachine file in standard locations')
+                return None
+
+            print(f'    Loading BedMachine from: {bedmachine_file}')
+            ds_bed = xr.open_dataset(bedmachine_file)
+
+            # Extract ice thickness and bed elevation
+            # BedMachine convention: bed < 0 is below sea level
+            thickness = ds_bed.thickness.values
+            bed = ds_bed.bed.values
+
+            # Get coordinates
+            x_bed = ds_bed.x.values
+            y_bed = ds_bed.y.values
+
+            # Convert polar stereographic to lat-lon
+            # BedMachine uses EPSG:3031 (Antarctic Polar Stereographic)
+            from pyproj import Transformer
+            transformer = Transformer.from_crs("EPSG:3031", "EPSG:4326")
+
+            # Create coordinate grids (subsample for efficiency)
+            skip = 5  # Use every 5th point
+            x_grid, y_grid = np.meshgrid(x_bed[::skip], y_bed[::skip])
+            lat_bed, lon_bed = transformer.transform(x_grid.flatten(), y_grid.flatten())
+            lat_bed = lat_bed.reshape(x_grid.shape)
+            lon_bed = lon_bed.reshape(x_grid.shape)
+
+            # Subsample thickness and bed
+            thickness_sub = thickness[::skip, ::skip]
+            bed_sub = bed[::skip, ::skip]
+
+            # Calculate flotation criterion
+            # Ice is grounded when: thickness > -bed * (ρ_ocean / ρ_ice)
+            # (only where bed < 0)
+            rho_ocean = 1028.0
+            rho_ice = 918.0
+
+            print(f'    Calculating flotation (ρ_ocean={rho_ocean}, ρ_ice={rho_ice})...')
+            thickness_flotation = np.where(bed_sub < 0,
+                                          -bed_sub * (rho_ocean / rho_ice),
+                                          np.inf)  # Always grounded if bed > 0
+
+            # Mask for grounded ice
+            grounded = thickness_sub > thickness_flotation
+
+            # Find grounding line as edge of grounded region
+            # Use morphological gradient to find boundaries
+            print('    Identifying grounding line cells...')
+            from scipy.ndimage import binary_dilation, binary_erosion
+
+            # Grounding line = grounded cells adjacent to floating cells
+            grounded_dilated = binary_dilation(grounded)
+            grounded_eroded = binary_erosion(grounded)
+            gl_mask = grounded_dilated & ~grounded_eroded
+
+            # Focus on Thwaites/Amundsen region
+            # Thwaites: roughly 75°S, 106°W to 74°S, 104°W
+            amundsen_mask = ((lat_bed >= -77) & (lat_bed <= -73) &
+                            (lon_bed >= -120) & (lon_bed <= -95))
+            gl_mask = gl_mask & amundsen_mask
+
+            if np.sum(gl_mask) == 0:
+                print('    No grounding line cells found in Amundsen region')
+                return None
+
+            print(f'    Found {np.sum(gl_mask)} grounding line cells')
+
+            # Extract grounding line coordinates
+            gl_lons = lon_bed[gl_mask]
+            gl_lats = lat_bed[gl_mask]
+
+            # Create a simple polygon/multipoint feature
+            # For simplicity, create a buffered multipoint
+            from shapely.geometry import MultiPoint, Point
+            import geojson
+
+            points = [Point(lon, lat) for lon, lat in zip(gl_lons, gl_lats)]
+            multipoint = MultiPoint(points)
+
+            # Buffer by ~2km to create polygon
+            # At ~75°S, 1 degree ≈ 30 km, so 2km ≈ 0.067 degrees
+            buffer_deg = 0.1
+            gl_polygon = multipoint.buffer(buffer_deg)
+
+            # Convert to geojson
+            feature = geojson.Feature(geometry=gl_polygon, properties={'name': 'Thwaites_GL'})
+            feature_collection = geojson.FeatureCollection([feature])
+
+            # Write to temp file and read back with geometric_features
+            import tempfile
+            with tempfile.NamedTemporaryFile(mode='w', suffix='.geojson', delete=False) as f:
+                geojson.dump(feature_collection, f)
+                temp_geojson = f.name
+
+            from geometric_features import read_feature_collection
+            fc = read_feature_collection(temp_geojson)
+
+            # Clean up temp file
+            os.remove(temp_geojson)
+
+            print('    Successfully created grounding line feature collection')
+            return fc
+
+        except Exception as e:
+            print(f'    Error extracting grounding line: {e}')
+            import traceback
+            traceback.print_exc()
+            return None
+
+    def _apply_gaussian_refinement(self, cellWidth, lon, lat,
+                                   res_gz, res_cavity, res_shelf):
+        """
+        Apply simple Gaussian refinement around Thwaites location (fallback)
+        """
+        thwaites_lat = -75.0
+        thwaites_lon = -106.0
+
+        lon_grid, lat_grid = np.meshgrid(lon, lat)
+
+        # Simple distance-based refinement
+        dist_lat = (lat_grid - thwaites_lat) ** 2
+        dist_lon = ((lon_grid - thwaites_lon) * np.cos(np.radians(lat_grid))) ** 2
+        approx_dist_deg = np.sqrt(dist_lat + dist_lon)
+
+        gz_dist = 0.2
+        cavity_dist = 1.0
+        shelf_dist = 3.0
+        transition_width = 0.1
+
+        # Apply refinements
+        mask_gz = 0.5 * (1 + np.tanh((approx_dist_deg - gz_dist) / transition_width))
+        cellWidth = res_gz * (1 - mask_gz) + cellWidth * mask_gz
+
+        mask_cavity = 0.5 * (1 + np.tanh((approx_dist_deg - cavity_dist) / transition_width))
+        cellWidth = np.minimum(cellWidth,
+                              res_cavity * (1 - mask_cavity) + cellWidth * mask_cavity)
+
+        mask_shelf = 0.5 * (1 + np.tanh((approx_dist_deg - shelf_dist) / transition_width))
+        cellWidth = np.minimum(cellWidth,
+                              res_shelf * (1 - mask_shelf) + cellWidth * mask_shelf)
+
+        return cellWidth
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/atlantic.geojson b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/atlantic.geojson
new file mode 100644
index 0000000000..227a57f066
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/atlantic.geojson
@@ -0,0 +1,97 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "properties": {
+        "name": "Atlantic region",
+        "component": "ocean",
+        "object": "region",
+        "author": "Xylar Asay-Davis"
+      },
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -97.3828125,
+              85.05112877979998
+            ],
+            [
+              -102.3046875,
+              40.17887331434696
+            ],
+            [
+              -102.3046875,
+              23.241346102386135
+            ],
+            [
+              -93.1640625,
+              15.623036831528264
+            ],
+            [
+              -85.78125,
+              13.581920900545844
+            ],
+            [
+              -83.583984375,
+              9.535748998133627
+            ],
+            [
+              -81.2109375,
+              8.059229627200192
+            ],
+            [
+              -79.013671875,
+              9.795677582829743
+            ],
+            [
+              -75.9375,
+              5.61598581915534
+            ],
+            [
+              -77.6953125,
+              0
+            ],
+            [
+              16.171875,
+              0
+            ],
+            [
+              27.773437499999996,
+              26.745610382199022
+            ],
+            [
+              37.96875,
+              32.24997445586331
+            ],
+            [
+              39.7265625,
+              39.36827914916014
+            ],
+            [
+              32.6953125,
+              53.9560855309879
+            ],
+            [
+              37.6171875,
+              61.438767493682825
+            ],
+            [
+              25.664062500000004,
+              68.26938680456564
+            ],
+            [
+              24.609375,
+              85.05112877979998
+            ],
+            [
+              -97.3828125,
+              85.05112877979998
+            ]
+          ]
+        ]
+      }
+    }
+  ]
+}
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/cull_mesh.py b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/cull_mesh.py
new file mode 100644
index 0000000000..c281454436
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/cull_mesh.py
@@ -0,0 +1,176 @@
+import os
+import xarray as xr
+import numpy as np
+from geometric_features import read_feature_collection
+from mpas_tools.mesh.creation.signed_distance import mask_from_geojson
+
+from compass.ocean.mesh.cull import CullMeshStep
+
+
+class ThwaitesCullMeshStep(CullMeshStep):
+    """
+    Custom cull mesh step for Thwaites mesh with regional domain support
+
+    This subclass extends the standard CullMeshStep to add regional culling.
+    It creates a hook file that the modified cull.py will execute to modify
+    the land mask after creation but before culling.
+    """
+
+    def setup(self):
+        """
+        Set up the step by adding input and output files
+        """
+        # Call parent setup
+        super().setup()
+
+        # Create hook file for land mask modification
+        self._create_hook_file()
+
+    def _create_hook_file(self):
+        """
+        Create the modify_land_mask_hook.py file that will be called
+        during the culling process
+        """
+        hook_content = '''"""
+Hook for modifying land mask before culling - auto-generated
+"""
+import os
+import xarray as xr
+import numpy as np
+from geometric_features import read_feature_collection
+from mpas_tools.mesh.creation.signed_distance import mask_from_geojson
+
+
+def modify_land_mask(logger):
+    """
+    Modify land_mask.nc to add regional domain boundaries
+
+    This is called from compass.ocean.mesh.cull._cull_mesh_with_logging
+    after land_mask.nc is created but before culling begins.
+    """
+    # Read config - need to get it from the step's config
+    # This hook runs in the cull_mesh working directory
+    from compass.config import CompassConfigParser
+    config = CompassConfigParser()
+    config.read('../compass.cfg')  # Read from test case level
+
+    # Check config section
+    if not config.has_section('thwaites01to60'):
+        logger.info('No thwaites01to60 config section')
+        return
+
+    section = config['thwaites01to60']
+
+    # Check if regional domain is specified
+    has_geojson = section.has_option('regional_domain_geojson')
+    has_bounds = (section.has_option('lat_min') and
+                  section.has_option('lon_min'))
+
+    if not (has_geojson or has_bounds):
+        logger.info('No regional domain specified - using standard land mask')
+        return
+
+    logger.info('====================================================')
+    logger.info('APPLYING REGIONAL MASK MODIFICATION FOR THWAITES')
+    logger.info('====================================================')
+
+    # Load base mesh to get cell coordinates
+    ds_mesh = xr.open_dataset('base_mesh.nc')
+
+    # Load land mask
+    ds_mask = xr.open_dataset('land_mask.nc')
+
+    # Find the mask variable
+    if 'regionCellMasks' in ds_mask:
+        mask_var = 'regionCellMasks'
+    elif 'landIceMask' in ds_mask:
+        mask_var = 'landIceMask'
+    else:
+        mask_vars = [v for v in ds_mask.variables if 'mask' in v.lower()]
+        if mask_vars:
+            mask_var = mask_vars[0]
+            logger.warning(f'Using mask variable: {mask_var}')
+        else:
+            logger.error('No mask variable found in land_mask.nc')
+            return
+
+    land_mask = ds_mask[mask_var]
+
+    # Create mask for cells INSIDE the regional domain
+    if has_geojson:
+        geojson_file = section.get('regional_domain_geojson')
+        logger.info(f'  Using domain from {geojson_file}')
+        fc = read_feature_collection(geojson_file)
+        inside_region = mask_from_geojson(
+            fc,
+            np.degrees(ds_mesh.lonCell.values),
+            np.degrees(ds_mesh.latCell.values))
+    else:
+        # Use lat-lon bounds
+        lat_min = section.getfloat('lat_min')
+        lat_max = section.getfloat('lat_max')
+        lon_min = section.getfloat('lon_min')
+        lon_max = section.getfloat('lon_max')
+
+        logger.info(f'  Using domain bounds:')
+        logger.info(f'    lat: [{lat_min}, {lat_max}]')
+        logger.info(f'    lon: [{lon_min}, {lon_max}]')
+
+        lat_deg = np.degrees(ds_mesh.latCell.values)
+        lon_deg = np.degrees(ds_mesh.lonCell.values)
+
+        # Handle longitude wrapping
+        lon_deg = np.where(lon_deg > 180, lon_deg - 360, lon_deg)
+
+        inside_region = ((lat_deg >= lat_min) & (lat_deg <= lat_max) &
+                        (lon_deg >= lon_min) & (lon_deg <= lon_max))
+
+    # Cells OUTSIDE the region get added to land mask
+    outside_region = ~inside_region
+
+    # Count cells
+    ncells_total = len(inside_region)
+    ncells_outside = int(outside_region.sum())
+    ncells_inside = int(inside_region.sum())
+    ncells_originally_land = int((land_mask.values > 0).sum())
+
+    logger.info(f'  Mesh statistics:')
+    logger.info(f'    Total cells: {ncells_total}')
+    logger.info(f'    Cells inside region: {ncells_inside}')
+    logger.info(f'    Cells outside region: {ncells_outside}')
+    logger.info(f'    Originally land: {ncells_originally_land}')
+
+    # Modify land mask
+    if np.issubdtype(land_mask.dtype, np.integer):
+        modified_mask = np.where(outside_region, 1, land_mask.values)
+    else:
+        modified_mask = np.logical_or(outside_region, land_mask.values > 0)
+
+    ncells_new_land = int((modified_mask > 0).sum())
+    logger.info(f'    After modification: {ncells_new_land} land cells')
+    logger.info(f'    Added {ncells_new_land - ncells_originally_land} cells to land mask')
+
+    # Backup original
+    if os.path.exists('land_mask.nc'):
+        os.rename('land_mask.nc', 'land_mask_original.nc')
+
+    # Update dataset
+    ds_mask[mask_var] = (land_mask.dims, modified_mask)
+    ds_mask[mask_var].attrs.update(land_mask.attrs)
+    ds_mask[mask_var].attrs['comment'] = (
+        'Land mask modified to include regional domain boundaries for Thwaites mesh')
+
+    # Write modified mask
+    ds_mask.to_netcdf('land_mask.nc')
+
+    logger.info('  Modified land mask written successfully')
+    logger.info('  (Original backed up to land_mask_original.nc)')
+    logger.info('====================================================')
+'''
+
+        # Write hook file to the step's path
+        # This will be in the cull_mesh subdirectory when the step runs
+        hook_file = os.path.join(self.path, 'modify_land_mask_hook.py')
+        os.makedirs(os.path.dirname(hook_file), exist_ok=True)
+        with open(hook_file, 'w') as f:
+            f.write(hook_content)
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/dynamic_adjustment.yaml b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/dynamic_adjustment.yaml
new file mode 100644
index 0000000000..5f46356813
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/dynamic_adjustment.yaml
@@ -0,0 +1,44 @@
+dynamic_adjustment:
+  land_ice_flux_mode: pressure_only
+  get_dt_from_min_res: False
+
+  steps:
+    damped_adjustment_1:
+      run_duration: 02_00:00:00
+      output_interval: 01_00:00:00
+      restart_interval: 02_00:00:00
+      dt: 00:00:30
+      btr_dt: 00:00:01
+      Rayleigh_damping_coeff: 1.0e-4
+
+    damped_adjustment_2:
+      run_duration: 08_00:00:00
+      output_interval: 10_00:00:00
+      restart_interval: 02_00:00:00
+      dt: 00:00:30
+      btr_dt: 00:00:01
+      Rayleigh_damping_coeff: 1.0e-5
+
+    damped_adjustment_3:
+      run_duration: 10_00:00:00
+      output_interval: 10_00:00:00
+      restart_interval: 10_00:00:00
+      dt: 00:00:30
+      btr_dt: 00:00:01
+      Rayleigh_damping_coeff: 1.0e-6
+
+    damped_adjustment_4:
+      run_duration: 20_00:00:00
+      output_interval: 10_00:00:00
+      restart_interval: 10_00:00:00
+      dt: 00:00:30
+      btr_dt: 00:00:01
+      Rayleigh_damping_coeff: None
+
+    simulation:
+      run_duration: 10_00:00:00
+      output_interval: 10_00:00:00
+      restart_interval: 10_00:00:00
+      dt: 00:00:50
+      btr_dt: 00:00:01.7
+      Rayleigh_damping_coeff: None
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris01to60.cfg b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris01to60.cfg
new file mode 100644
index 0000000000..c9f3cadb15
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris01to60.cfg
@@ -0,0 +1,86 @@
+# Options related to the vertical grid
+[vertical_grid]
+
+# the type of vertical grid
+grid_type = index_tanh_dz
+
+# Number of vertical levels
+vert_levels = 64
+
+# Depth of the bottom of the ocean
+bottom_depth = 5500.0
+
+# The minimum layer thickness
+min_layer_thickness = 10.0
+
+# The maximum layer thickness
+max_layer_thickness = 250.0
+
+# The characteristic number of levels over which the transition between
+# the min and max occurs
+transition_levels = 28
+
+# Anvil runs out of memory so do the following
+# config options related to remapping topography to an MPAS-Ocean mesh
+[remap_topography]
+
+# the target and minimum number of MPI tasks to use in remapping
+ntasks = 4096
+min_tasks = 360
+
+# options for global ocean testcases
+[global_ocean]
+
+## config options related to the initial_state step
+# number of cores to use
+init_ntasks = 512
+# minimum of cores, below which the step fails
+init_min_tasks = 64
+# The number of cores per task in init mode -- used to avoid running out of
+# memory where needed
+init_cpus_per_task = 4
+# whether to update PIO tasks and stride
+init_update_pio = False
+
+# whether to update PIO tasks and stride
+forward_update_pio = False
+
+# the approximate number of cells in the mesh
+approx_cell_count = 2200000# time step per resolution (s/km), since dt is proportional to resolution
+
+dt_per_km = 5
+# barotropic time step per resolution (s/km)
+btr_dt_per_km = 0.2
+
+## metadata related to the mesh
+# the prefix (e.g. QU, EC, WC, SO)
+prefix = FRIS
+# a description of the mesh and initial condition
+mesh_description = MPAS Southern Ocean regionally refined mesh for E3SM version
+                   ${e3sm_version} with enhanced resolution (12 km) around
+                   Antarctica and even more enhanced around the Filchner-Ronne Ice
+                   Shelf (${min_res} km), 45-km resolution in the mid southern latitudes,
+                   30-km resolution in a 15-degree band around the equator, 60-km
+                   resolution in northern mid latitudes, 30 km in the north
+                   Atlantic and 35 km in the Arctic.  This mesh has <<<levels>>>
+                   vertical levels and includes cavities under the ice shelves
+                   around Antarctica.
+# E3SM version that the mesh is intended for
+e3sm_version = 3
+# The revision number of the mesh, which should be incremented each time the
+# mesh is revised
+mesh_revision = 1
+# the minimum (finest) resolution in the mesh
+min_res = 1
+# the maximum (coarsest) resolution in the mesh, can be the same as min_res
+max_res = 60
+# The URL of the pull request documenting the creation of the mesh
+pull_request = N/A
+
+
+# config options related to initial condition and diagnostics support files
+# for E3SM
+[files_for_e3sm]
+
+# CMIP6 grid resolution
+cmip6_grid_res = 180x360
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1.geojson b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1.geojson
new file mode 100644
index 0000000000..51e767e649
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1.geojson
@@ -0,0 +1 @@
+{"type":"FeatureCollection","features":[{"type":"Feature","properties":{"name":"FRIS high res region","component":"ocean","object":"region","author":"Irena Vankova"},"geometry":{"coordinates":[[[-64.7438723001555,-72.8054024693636],[-69.76898821706332,-74.19251323259206],[-74.02484092117425,-74.47843628310608],[-92.43607426605669,-76.75259948500253],[-101.5379406548874,-80.18910132153749],[-97.25418139197761,-83.23536684168135],[-56.69746887655441,-85.05112877980659],[-39.90477420089397,-85],[-14.020185809699171,-83.36577201994825],[-20.52300577998716,-79.88944960088696],[-17.202195820003567,-76.59018199209325],[-14.53780624392931,-74.27113220357664],[-14.621251472019907,-72.71417212450604],[-15.457525929908314,-72.11805949893571],[-16.735628416005596,-71.3506375575714],[-18.81770866866242,-70.89188604300722],[-21.024357086561025,-70.9256187691209],[-24.78630070550588,-71.36235128369906],[-27.81675331568485,-72.1132700200943],[-31.314772897052507,-72.29684652688603],[-35.143774445838886,-72.21433963956319],[-39.97348523222354,-71.75366563644273],[-44.89032326598664,-70.87082796536023],[-50.40043324305111,-70.2508442400395],[-54.354093738939596,-70.20075756173516],[-58.7763090757868,-71.27289439521846],[-64.7438723001555,-72.8054024693636]]],"type":"Polygon"}}]}
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_correction_peninsula.geojson b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_correction_peninsula.geojson
new file mode 100644
index 0000000000..01859be40c
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_correction_peninsula.geojson
@@ -0,0 +1,73 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "properties": {
+        "name": "FRIS - peninsula correction",
+        "component": "ocean",
+        "object": "region",
+        "author": "Irena Vankova"
+      },
+      "geometry": {
+        "coordinates": [
+          [
+            [
+              -57.51168188592774,
+              -63.34469553787141
+            ],
+            [
+              -75.7308231480595,
+              -67.71524986365019
+            ],
+            [
+              -82.23302312713369,
+              -71.9805712785419
+            ],
+            [
+              -81.63135111895728,
+              -75.59567609845988
+            ],
+            [
+              -75.9055834859653,
+              -75.34148074887518
+            ],
+            [
+              -72.47014717927561,
+              -74.7673143097562
+            ],
+            [
+              -68.00554191414759,
+              -73.77984973214085
+            ],
+            [
+              -65.3550099151992,
+              -71.77205397981956
+            ],
+            [
+              -66.27933986206787,
+              -70.52162831122558
+            ],
+            [
+              -64.58699491448971,
+              -69.47546835230543
+            ],
+            [
+              -65.97822782804397,
+              -67.77544209184867
+            ],
+            [
+              -64.15647518604584,
+              -65.7193069906435
+            ],
+            [
+              -57.51168188592774,
+              -63.34469553787141
+            ]
+          ]
+        ],
+        "type": "Polygon"
+      }
+    }
+  ]
+}
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_correction_peninsula_v2.geojson b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_correction_peninsula_v2.geojson
new file mode 100644
index 0000000000..22d512589f
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_correction_peninsula_v2.geojson
@@ -0,0 +1,97 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "properties": {
+        "name": "FRIS - peninsula correction",
+        "component": "ocean",
+        "object": "region",
+        "author": "Irena Vankova"
+      },
+      "geometry": {
+        "coordinates": [
+          [
+            [
+              -57.51168188592774,
+              -63.34469553787141
+            ],
+            [
+              -53.05987330787791,
+              -60.9791790833165
+            ],
+            [
+              -53.06885379594741,
+              -59.10724141600171
+            ],
+            [
+              -69.3124237875646,
+              -61.21460931439017
+            ],
+            [
+              -83.38421745911926,
+              -62.77272333644538
+            ],
+            [
+              -97.90929281723952,
+              -66.221448196283
+            ],
+            [
+              -105.3263003641532,
+              -68.89292813601311
+            ],
+            [
+              -106.3307246723799,
+              -72.31928061089278
+            ],
+            [
+              -96.38678419860847,
+              -74.86775750756654
+            ],
+            [
+              -80.96950547249958,
+              -74.8917003828374
+            ],
+            [
+              -75.9055834859653,
+              -75.34148074887518
+            ],
+            [
+              -72.33775353811286,
+              -75.0153658885297
+            ],
+            [
+              -66.56769783304944,
+              -74.21440843603395
+            ],
+            [
+              -64.79776921678202,
+              -72.25395416577506
+            ],
+            [
+              -65.48170522640746,
+              -70.69919310630064
+            ],
+            [
+              -64.58699491448971,
+              -69.47546835230543
+            ],
+            [
+              -65.5606241702666,
+              -67.78677626662034
+            ],
+            [
+              -64.15647518604584,
+              -65.7193069906435
+            ],
+            [
+              -57.51168188592774,
+              -63.34469553787141
+            ]
+          ]
+        ],
+        "type": "Polygon"
+      }
+    }
+  ]
+}
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_correction_peninsula_v3.geojson b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_correction_peninsula_v3.geojson
new file mode 100644
index 0000000000..3c79f8ae21
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_correction_peninsula_v3.geojson
@@ -0,0 +1,97 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "properties": {
+        "name": "FRIS - peninsula correction",
+        "component": "ocean",
+        "object": "region",
+        "author": "Irena Vankova"
+      },
+      "geometry": {
+        "coordinates": [
+          [
+            [
+              -57.51168188592774,
+              -63.34469553787141
+            ],
+            [
+              -53.05987330787791,
+              -60.9791790833165
+            ],
+            [
+              -53.06885379594741,
+              -59.10724141600171
+            ],
+            [
+              -69.14223788466518,
+              -62.14576624144387
+            ],
+            [
+              -77.16348450278696,
+              -64.13666710372263
+            ],
+            [
+              -85.897068042515,
+              -66.25640705824505
+            ],
+            [
+              -94.59878299061488,
+              -69.72254279603585
+            ],
+            [
+              -100.24805622281325,
+              -72.26566396175254
+            ],
+            [
+              -96.38678419860847,
+              -74.86775750756654
+            ],
+            [
+              -80.96950547249958,
+              -74.8917003828374
+            ],
+            [
+              -75.9055834859653,
+              -75.34148074887518
+            ],
+            [
+              -72.33775353811286,
+              -75.0153658885297
+            ],
+            [
+              -66.56769783304944,
+              -74.21440843603395
+            ],
+            [
+              -64.79776921678202,
+              -72.25395416577506
+            ],
+            [
+              -65.48170522640746,
+              -70.69919310630064
+            ],
+            [
+              -64.58699491448971,
+              -69.47546835230543
+            ],
+            [
+              -65.5606241702666,
+              -67.78677626662034
+            ],
+            [
+              -64.15647518604584,
+              -65.7193069906435
+            ],
+            [
+              -57.51168188592774,
+              -63.34469553787141
+            ]
+          ]
+        ],
+        "type": "Polygon"
+      }
+    }
+  ]
+}
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_peninsula_12km.geojson b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_peninsula_12km.geojson
new file mode 100644
index 0000000000..7a4aec28d1
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_peninsula_12km.geojson
@@ -0,0 +1,93 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "properties": {
+        "name": "FRIS - peninsula correction",
+        "component": "ocean",
+        "object": "region",
+        "author": "Irena Vankova"
+      },
+      "geometry": {
+        "coordinates": [
+          [
+            [
+              -57.69353078184756,
+              -61.96708124373723
+            ],
+            [
+              -61.80022030800379,
+              -63.008208452662345
+            ],
+            [
+              -67.32103132588126,
+              -64.04079340739024
+            ],
+            [
+              -78.16701758864495,
+              -66.88371277926564
+            ],
+            [
+              -87.41765762163665,
+              -70.17410005159073
+            ],
+            [
+              -92.77451946184793,
+              -72.57628695883871
+            ],
+            [
+              -90.82481328167913,
+              -74.34659424903947
+            ],
+            [
+              -81.2027800330412,
+              -74.8636788514809
+            ],
+            [
+              -75.9055834859653,
+              -75.34148074887518
+            ],
+            [
+              -72.22449654558793,
+              -75.015667940154
+            ],
+            [
+              -66.58805223365164,
+              -74.20709557717424
+            ],
+            [
+              -64.5942902255466,
+              -72.25137878330746
+            ],
+            [
+              -65.51019680255283,
+              -70.6800717849537
+            ],
+            [
+              -64.58699491448971,
+              -69.47546835230543
+            ],
+            [
+              -65.55143595844069,
+              -67.79157062052742
+            ],
+            [
+              -64.15647518604584,
+              -65.7193069906435
+            ],
+            [
+              -58.35195871411565,
+              -63.31306460865669
+            ],
+            [
+              -57.69353078184756,
+              -61.96708124373723
+            ]
+          ]
+        ],
+        "type": "Polygon"
+      }
+    }
+  ]
+}
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_peninsula_12km_transition.geojson b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_peninsula_12km_transition.geojson
new file mode 100644
index 0000000000..b8aaed8edb
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_peninsula_12km_transition.geojson
@@ -0,0 +1,89 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "properties": {
+        "name": "FRIS - peninsula correction",
+        "component": "ocean",
+        "object": "region",
+        "author": "Irena Vankova"
+      },
+      "geometry": {
+        "coordinates": [
+          [
+            [
+              -54.54564191055567,
+              -63.1725569390185
+            ],
+            [
+              -53.05987330787791,
+              -60.9791790833165
+            ],
+            [
+              -53.06885379594741,
+              -59.10724141600171
+            ],
+            [
+              -69.14223788466518,
+              -62.14576624144387
+            ],
+            [
+              -77.16348450278696,
+              -64.13666710372263
+            ],
+            [
+              -85.897068042515,
+              -66.25640705824505
+            ],
+            [
+              -94.59878299061488,
+              -69.72254279603585
+            ],
+            [
+              -100.24805622281325,
+              -72.26566396175254
+            ],
+            [
+              -96.38678419860847,
+              -74.86775750756654
+            ],
+            [
+              -89.07987805140462,
+              -73.74748954581011
+            ],
+            [
+              -83.42453978625146,
+              -72.01434576499264
+            ],
+            [
+              -78.97666329388788,
+              -70.01288831317304
+            ],
+            [
+              -73.84932654865166,
+              -68.18364067984133
+            ],
+            [
+              -68.87563501264184,
+              -66.90585850491902
+            ],
+            [
+              -64.15647518604584,
+              -65.7193069906435
+            ],
+            [
+              -59.3178981972859,
+              -64.13380186083748
+            ],
+            [
+              -54.54564191055567,
+              -63.1725569390185
+            ]
+          ]
+        ],
+        "type": "Polygon"
+      }
+    }
+  ]
+}
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_peninsula_12km_v2.geojson b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_peninsula_12km_v2.geojson
new file mode 100644
index 0000000000..aa2a173208
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_peninsula_12km_v2.geojson
@@ -0,0 +1,97 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "properties": {
+        "name": "FRIS - peninsula correction",
+        "component": "ocean",
+        "object": "region",
+        "author": "Irena Vankova"
+      },
+      "geometry": {
+        "coordinates": [
+          [
+            [
+              -57.69353078184756,
+              -61.96708124373723
+            ],
+            [
+              -60.60980216459852,
+              -61.55452131445545
+            ],
+            [
+              -69.55211766107989,
+              -63.152080989803714
+            ],
+            [
+              -80.85673116951071,
+              -66.05555597379733
+            ],
+            [
+              -88.13681707469551,
+              -68.58641430243729
+            ],
+            [
+              -95.69478642361177,
+              -71.59338385976186
+            ],
+            [
+              -96.62401026588171,
+              -72.86186735741381
+            ],
+            [
+              -90.82481328167913,
+              -74.34659424903947
+            ],
+            [
+              -81.2027800330412,
+              -74.8636788514809
+            ],
+            [
+              -75.9055834859653,
+              -75.34148074887518
+            ],
+            [
+              -72.22449654558793,
+              -75.015667940154
+            ],
+            [
+              -66.58805223365164,
+              -74.20709557717424
+            ],
+            [
+              -64.5942902255466,
+              -72.25137878330746
+            ],
+            [
+              -65.51019680255283,
+              -70.6800717849537
+            ],
+            [
+              -64.58699491448971,
+              -69.47546835230543
+            ],
+            [
+              -65.55143595844069,
+              -67.79157062052742
+            ],
+            [
+              -64.15647518604584,
+              -65.7193069906435
+            ],
+            [
+              -58.35195871411565,
+              -63.31306460865669
+            ],
+            [
+              -57.69353078184756,
+              -61.96708124373723
+            ]
+          ]
+        ],
+        "type": "Polygon"
+      }
+    }
+  ]
+}
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_transition.geojson b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_transition.geojson
new file mode 100644
index 0000000000..20ad8b8d9b
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/fris_v1_transition.geojson
@@ -0,0 +1,93 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "properties": {
+        "name": "FRIS - transition region",
+        "component": "ocean",
+        "object": "region",
+        "author": "Irena Vankova"
+      },
+      "geometry": {
+        "coordinates": [
+          [
+            [
+              -57.51168188592774,
+              -63.34469553787141
+            ],
+            [
+              -75.7308231480595,
+              -67.71524986365019
+            ],
+            [
+              -82.23302312713369,
+              -71.9805712785419
+            ],
+            [
+              -81.63135111895728,
+              -75.59567609845988
+            ],
+            [
+              -72.43392556848906,
+              -77.89255200360182
+            ],
+            [
+              -58.71407611692999,
+              -77.78724511207038
+            ],
+            [
+              -46.208104421831024,
+              -79.00092023048762
+            ],
+            [
+              -18.338759042014715,
+              -79.45787446576274
+            ],
+            [
+              -2.7857142568166466,
+              -77.44116982233017
+            ],
+            [
+              6.813437365949312,
+              -73.42744971284144
+            ],
+            [
+              5.00921262812011,
+              -68.98423794486757
+            ],
+            [
+              -0.7923051110294352,
+              -65.9398557197854
+            ],
+            [
+              -10.402032107542851,
+              -64.99932710553239
+            ],
+            [
+              -18.521043388561708,
+              -64.3196282813596
+            ],
+            [
+              -29.861824936615392,
+              -65.13348750621226
+            ],
+            [
+              -40.06168534685344,
+              -64.54083717367145
+            ],
+            [
+              -50.97065707925094,
+              -63.2302794039998
+            ],
+            [
+              -57.51168188592774,
+              -63.34469553787141
+            ]
+          ]
+        ],
+        "type": "Polygon"
+      }
+    }
+  ]
+}
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/high_res_region.geojson b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/high_res_region.geojson
new file mode 100644
index 0000000000..a536ebaba9
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/high_res_region.geojson
@@ -0,0 +1,73 @@
+{
+  "type": "FeatureCollection",
+  "features": [
+    {
+      "type": "Feature",
+      "properties": {
+        "name": "SO12to60 high res region",
+        "component": "ocean",
+        "object": "region",
+        "author": "Xylar Asay-Davis"
+      },
+      "geometry": {
+        "type": "Polygon",
+        "coordinates": [
+          [
+            [
+              -75.5859375,
+              -48.92249926375823
+            ],
+            [
+              -96.50390625,
+              -54.67383096593114
+            ],
+            [
+              -124.45312499999999,
+              -53.95608553098789
+            ],
+            [
+              -180,
+              -53
+            ],
+            [
+              -180,
+              -90
+            ],
+            [
+              180,
+              -90
+            ],
+            [
+              180,
+              -53
+            ],
+            [
+              168.3984375,
+              -51.17934297928927
+            ],
+            [
+              121.640625,
+              -45.82879925192133
+            ],
+            [
+              22.148437499999996,
+              -37.99616267972812
+            ],
+            [
+              -61.17187499999999,
+              -34.30714385628803
+            ],
+            [
+              -68.90625,
+              -40.97989806962013
+            ],
+            [
+              -75.5859375,
+              -48.92249926375823
+            ]
+          ]
+        ]
+      }
+    }
+  ]
+}
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/modify_land_mask.py b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/modify_land_mask.py
new file mode 100644
index 0000000000..c72703a2ee
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/modify_land_mask.py
@@ -0,0 +1,188 @@
+from compass.step import Step
+import xarray as xr
+from geometric_features import read_feature_collection
+from mpas_tools.mesh.creation.signed_distance import mask_from_geojson
+import numpy as np
+
+
+class ModifyLandMask(Step):
+    """
+    Modify land mask to add regional domain boundaries
+
+    This step runs BEFORE the standard cull_mesh step. If a regional domain
+    is specified (via geojson or lat-lon bounds), cells outside the domain
+    are marked as "land" in the mask. This causes them to be culled in the
+    standard cull_mesh step that follows.
+
+    If no regional domain is specified, this step passes the land mask
+    through unchanged.
+    """
+
+    def __init__(self, test_case, base_mesh_step, name='modify_land_mask',
+                 subdir=None):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        test_case : compass.TestCase
+            The test case this step belongs to
+
+        base_mesh_step : compass.Step
+            The base mesh generation step
+
+        name : str, optional
+            The name of the step
+
+        subdir : str, optional
+            The subdirectory for the step
+        """
+        super().__init__(test_case=test_case, name=name, subdir=subdir)
+
+        self.base_mesh_step = base_mesh_step
+
+    def setup(self):
+        """
+        Set up the step by adding input and output files
+        """
+        # We need the base mesh to get cell coordinates
+        base_mesh_path = self.base_mesh_step.path
+        base_mesh_filename = self.base_mesh_step.config.get(
+            'spherical_mesh', 'mpas_mesh_filename')
+        self.add_input_file(
+            filename='base_mesh.nc',
+            work_dir_target=f'{base_mesh_path}/{base_mesh_filename}')
+
+        # Input: land_mask.nc from geometric_features processing
+        # This should exist from the standard global_ocean workflow
+        self.add_input_file(filename='land_mask.nc',
+                           target='land_mask.nc')
+
+        # Output: modified land mask (or same if no regional domain)
+        self.add_output_file(filename='land_mask_with_region.nc')
+
+    def run(self):
+        """
+        Run the step - modify land mask if regional domain is specified
+        """
+        config = self.config
+        logger = self.logger
+
+        # Check config section - try both specific and general sections
+        if config.has_section('thwaites01to60'):
+            section_name = 'thwaites01to60'
+        elif config.has_section('spherical_mesh'):
+            section_name = 'spherical_mesh'
+        else:
+            section_name = None
+
+        # Check if regional domain is specified
+        has_geojson = (section_name and
+                      config.has_option(section_name, 'regional_domain_geojson'))
+        has_bounds = (section_name and
+                     config.has_option(section_name, 'lat_min') and
+                     config.has_option(section_name, 'lon_min'))
+
+        if not (has_geojson or has_bounds):
+            # No regional domain - pass through unchanged
+            logger.info('No regional domain specified - using standard land mask')
+            import os
+            os.symlink('land_mask.nc', 'land_mask_with_region.nc')
+            return
+
+        logger.info('Modifying land mask to add regional boundaries...')
+
+        # Load base mesh to get cell coordinates
+        ds_mesh = xr.open_dataset('base_mesh.nc')
+
+        # Load land mask
+        ds_mask = xr.open_dataset('land_mask.nc')
+
+        # Determine which mask variable to use
+        # Different global_ocean meshes may use different variable names
+        if 'regionCellMasks' in ds_mask:
+            mask_var = 'regionCellMasks'
+        elif 'landIceMask' in ds_mask:
+            mask_var = 'landIceMask'
+        else:
+            # Find any mask-like variable
+            mask_vars = [v for v in ds_mask.variables if 'mask' in v.lower()]
+            if mask_vars:
+                mask_var = mask_vars[0]
+                logger.warning(f'Using mask variable: {mask_var}')
+            else:
+                raise ValueError('No mask variable found in land_mask.nc')
+
+        land_mask = ds_mask[mask_var]
+
+        section = config[section_name]
+
+        # Create mask for cells INSIDE the regional domain
+        if has_geojson:
+            geojson_file = section.get('regional_domain_geojson')
+            logger.info(f'  Using domain from {geojson_file}')
+            fc = read_feature_collection(geojson_file)
+            inside_region = mask_from_geojson(
+                fc,
+                np.degrees(ds_mesh.lonCell.values),
+                np.degrees(ds_mesh.latCell.values))
+        else:
+            # Use lat-lon bounds
+            lat_min = section.getfloat('lat_min')
+            lat_max = section.getfloat('lat_max')
+            lon_min = section.getfloat('lon_min')
+            lon_max = section.getfloat('lon_max')
+
+            logger.info(f'  Using bounds: '
+                       f'lat [{lat_min}, {lat_max}], '
+                       f'lon [{lon_min}, {lon_max}]')
+
+            lat_deg = np.degrees(ds_mesh.latCell.values)
+            lon_deg = np.degrees(ds_mesh.lonCell.values)
+
+            # Handle longitude wrapping (convert to -180 to 180 range)
+            lon_deg = np.where(lon_deg > 180, lon_deg - 360, lon_deg)
+
+            inside_region = ((lat_deg >= lat_min) & (lat_deg <= lat_max) &
+                           (lon_deg >= lon_min) & (lon_deg <= lon_max))
+
+        # Cells OUTSIDE the region get added to land mask
+        # (i.e., they will be culled by the standard cull_mesh step)
+        outside_region = ~inside_region
+
+        # Count how many cells will be affected
+        ncells_total = len(inside_region)
+        ncells_outside = outside_region.sum()
+        ncells_inside = inside_region.sum()
+        ncells_already_land = (land_mask > 0).sum()
+
+        logger.info(f'  Total cells: {ncells_total}')
+        logger.info(f'  Cells inside region: {ncells_inside}')
+        logger.info(f'  Cells outside region: {ncells_outside}')
+        logger.info(f'  Cells already land: {ncells_already_land}')
+
+        # Modify land mask: original land + cells outside region
+        # Use logical OR to combine masks
+        if np.issubdtype(land_mask.dtype, np.integer):
+            # Integer mask: 1 = land, 0 = ocean
+            modified_mask = np.where(outside_region, 1, land_mask.values)
+        else:
+            # Boolean or float mask
+            modified_mask = np.logical_or(outside_region, land_mask.values > 0)
+
+        ncells_new_land = (modified_mask > 0).sum()
+        logger.info(f'  Cells marked as land after modification: {ncells_new_land}')
+
+        # Update dataset with modified mask
+        ds_mask[mask_var] = (land_mask.dims, modified_mask)
+
+        # Copy attributes
+        ds_mask[mask_var].attrs.update(land_mask.attrs)
+        ds_mask[mask_var].attrs['comment'] = (
+            'Land mask modified to include regional domain boundaries. '
+            'Cells outside the specified domain are marked as land.')
+
+        # Write output
+        ds_mask.to_netcdf('land_mask_with_region.nc')
+
+        logger.info('  Modified land mask written to land_mask_with_region.nc')
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/namelist.init b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/namelist.init
new file mode 100644
index 0000000000..fe4098d055
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/namelist.init
@@ -0,0 +1,8 @@
+config_rx1_inner_iter_count = 20
+config_rx1_horiz_smooth_weight = 10.0
+config_rx1_vert_smooth_weight = 10.0
+config_rx1_slope_weight = 1e-1
+config_rx1_zstar_weight = 10.0
+config_rx1_horiz_smooth_open_ocean_cells = 240
+config_rx1_min_levels = 5
+config_rx1_min_layer_thickness = 2
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/namelist.split_explicit_ab2 b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/namelist.split_explicit_ab2
new file mode 100644
index 0000000000..9193d6b176
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/namelist.split_explicit_ab2
@@ -0,0 +1,12 @@
+config_time_integrator = 'split_explicit_ab2'
+config_dt = '00:00:05'
+config_btr_dt = '00:00:00.17'
+config_run_duration = '0000_01:00:00'
+config_use_mom_del2 = .true.
+config_mom_del2 = 38.5
+config_use_mom_del4 = .true.
+config_mom_del4 = 6.83e6
+config_hmix_scaleWithMesh = .true.
+config_use_GM = .true.
+config_GM_closure = 'constant'
+config_GM_constant_kappa = 600.0
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/streams.forward b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/streams.forward
new file mode 100644
index 0000000000..bd42317610
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/streams.forward
@@ -0,0 +1,13 @@
+<streams>
+
+<stream name="output">
+
+    <var_struct name="tracers"/>
+    <var name="xtime"/>
+    <var name="normalVelocity"/>
+    <var name="layerThickness"/>
+    <var name="kineticEnergyCell"/>
+
+</stream>
+
+</streams>
diff --git a/compass/ocean/tests/global_ocean/mesh/thwaites01to60/thwaites01to60.cfg b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/thwaites01to60.cfg
new file mode 100644
index 0000000000..d98f18162b
--- /dev/null
+++ b/compass/ocean/tests/global_ocean/mesh/thwaites01to60/thwaites01to60.cfg
@@ -0,0 +1,40 @@
+# config options for thwaites01to60 mesh
+[thwaites01to60]
+
+# OPTIONAL: Regional domain (if not set, keeps full global mesh)
+# Uncomment ONE of these approaches to enable regional culling:
+
+# Option A: Use geojson polygon
+# regional_domain_geojson = amundsen_domain.geojson
+
+# Option B: Use lat-lon bounding box (Amundsen Sea sector)
+# Recommended for testing: covers Thwaites, Pine Island, and surrounding area
+lat_min = -76.0
+lat_max = -73.0
+lon_min = -116.0
+lon_max = -98.0
+
+# OPTIONAL: Thin film beneath grounded ice
+thin_film_present = False
+thin_film_thickness = 0.05
+haf_threshold = 30.0
+film_inland_extent = 30e3
+
+# Resolution parameters for base mesh generation (km)
+res_gz = 0.2
+res_cavity = 1.0
+res_shelf = 3.0
+res_far = 8.0
+gz_band_halfwidth = 15.0
+
+# the name of the topography file in the bathymetry database
+topo_filename = BedMachineAntarctica-v3_GEBCO_2023_ne3000_20250110.nc
+src_scrip_filename = ne3000_20250110.scrip.nc
+
+# description of the bathymetry
+description = Bathymetry is from GEBCO 2023, combined with BedMachine
+              Antarctica v3 (will update to v4), with Thwaites grounding zone refinement
+
+# Use BedMachine v4 when available
+[combine_topo]
+# antarctic_filename = BedMachineAntarctica-v4.nc