diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py
index 3402057e..82106992 100644
--- a/nlmod/dims/grid.py
+++ b/nlmod/dims/grid.py
@@ -44,6 +44,54 @@
 logger = logging.getLogger(__name__)
 
 
+def snap_extent(extent, delr, delc):
+    """
+    snap the extent in such a way that an integer number of columns and rows fit
+    in the extent. The new extent is always equal to, or bigger than the
+    original extent.
+
+    Parameters
+    ----------
+    extent : list, tuple or np.array
+        original extent (xmin, xmax, ymin, ymax)
+    delr : int or float,
+        cell size along rows, equal to dx
+    delc : int or float,
+        cell size along columns, equal to dy
+
+    Returns
+    -------
+    extent : list, tuple or np.array
+        adjusted extent
+    """
+    extent = list(extent).copy()
+
+    logger.debug(f"redefining extent: {extent}")
+
+    if delr <= 0 or delc <= 0:
+        raise ValueError("delr and delc should be positive values")
+
+    # if xmin can be divided by delr do nothing, otherwise rescale xmin
+    if not extent[0] % delr == 0:
+        extent[0] -= extent[0] % delr
+
+    # get number of columns and xmax
+    ncol = int(np.ceil((extent[1] - extent[0]) / delr))
+    extent[1] = extent[0] + (ncol * delr)  # round xmax up to close grid
+
+    # if ymin can be divided by delc do nothing, otherwise rescale ymin
+    if not extent[2] % delc == 0:
+        extent[2] -= extent[2] % delc
+
+    # get number of rows and ymax
+    nrow = int(np.ceil((extent[3] - extent[2]) / delc))
+    extent[3] = extent[2] + (nrow * delc)  # round ymax up to close grid
+
+    logger.debug(f"new extent is {extent} and has {nrow} rows and {ncol} columns")
+
+    return extent
+
+
 def xy_to_icell2d(xy, ds):
     """get the icell2d value of a point defined by its x and y coordinates.
 
diff --git a/pyproject.toml b/pyproject.toml
index 8464577d..548c6c22 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -23,7 +23,7 @@ requires-python = ">= 3.8"
 dependencies = [
     "flopy>=3.3.6",
     "xarray>=0.16.1",
-    "netcdf4>=1.5.7",
+    "netcdf4>=1.6.3",
     "rasterio>=1.1.0",
     "rioxarray",
     "affine>=0.3.1",
diff --git a/tests/test_001_model.py b/tests/test_001_model.py
index fde77874..d73130a0 100644
--- a/tests/test_001_model.py
+++ b/tests/test_001_model.py
@@ -17,6 +17,16 @@ def test_model_directories(tmpdir):
     figdir, cachedir = nlmod.util.get_model_dirs(model_ws)
 
 
+def test_snap_extent():
+    extent = (0.22, 1056.12, 7.43, 1101.567)
+    new_extent = nlmod.dims.snap_extent(extent, 10, 20)
+    assert new_extent == [0.0, 1060.0, 0.0, 1120.0]
+
+    extent = (1000, 2000, 8000, 10000)
+    new_extent = nlmod.dims.snap_extent(extent, 250, 55)
+    assert new_extent == [1000.0, 2000.0, 7975.0, 10010.0]
+
+
 def get_ds_time_steady(tmpdir, modelname="test"):
     model_ws = os.path.join(tmpdir, "test_model")
     ds = nlmod.base.set_ds_attrs(xr.Dataset(), modelname, model_ws)