diff --git a/nlmod/cache.py b/nlmod/cache.py
index 6ccd717f..b06f9995 100644
--- a/nlmod/cache.py
+++ b/nlmod/cache.py
@@ -692,12 +692,6 @@ def ds_contains(
             datavars.append("yv")
             datavars.append("icvert")
 
-            # TODO: temporary fix until delr/delc are completely removed
-            if "delr" in datavars:
-                datavars.remove("delr")
-            if "delc" in datavars:
-                datavars.remove("delc")
-
         if "angrot" in ds.attrs:
             attrs.append("angrot")
 
@@ -721,14 +715,16 @@ def ds_contains(
         if "time" not in ds.coords:
             msg = "time not in dataset. Run nlmod.time.set_ds_time() first."
             raise ValueError(msg)
-        
+
         # Check if time-coord is complete
         time_attrs_required = ["start", "time_units"]
 
         for t_attr in time_attrs_required:
             if t_attr not in ds["time"].attrs:
-                msg = f"{t_attr} not in dataset['time'].attrs. " +\
-                    "Run nlmod.time.set_ds_time() to set time."
+                msg = (
+                    f"{t_attr} not in dataset['time'].attrs. "
+                    + "Run nlmod.time.set_ds_time() to set time."
+                )
                 raise ValueError(msg)
 
     # User-unfriendly error messages
diff --git a/nlmod/dims/base.py b/nlmod/dims/base.py
index 3ec3f946..5de5b007 100644
--- a/nlmod/dims/base.py
+++ b/nlmod/dims/base.py
@@ -167,8 +167,6 @@ def to_model_ds(
         )
     elif isinstance(delr, np.ndarray) and isinstance(delc, np.ndarray):
         ds["area"] = ("y", "x"), np.outer(delc, delr)
-        ds["delr"] = ("x"), delr
-        ds["delc"] = ("y"), delc
     else:
         raise TypeError("unexpected type for delr and/or delc")
 
@@ -366,11 +364,6 @@ def _get_structured_grid_ds(
         coords=coords,
         attrs=attrs,
     )
-    # set delr and delc
-    delr = np.diff(xedges)
-    ds["delr"] = ("x"), delr
-    delc = -np.diff(yedges)
-    ds["delc"] = ("y"), delc
 
     if crs is not None:
         ds.rio.set_crs(crs)
diff --git a/nlmod/dims/grid.py b/nlmod/dims/grid.py
index 68a4c9b9..42ed87c6 100644
--- a/nlmod/dims/grid.py
+++ b/nlmod/dims/grid.py
@@ -40,6 +40,8 @@
     get_affine_mod_to_world,
     get_affine_world_to_mod,
     structured_da_to_ds,
+    get_delr,
+    get_delc,
 )
 
 logger = logging.getLogger(__name__)
@@ -211,17 +213,9 @@ def modelgrid_from_ds(ds, rotated=True, nlay=None, top=None, botm=None, **kwargs
             raise TypeError(
                 f"extent should be a list, tuple or numpy array, not {type(ds.extent)}"
             )
-        if "delc" in ds:
-            delc = ds["delc"].values
-        else:
-            raise KeyError("delc not in dataset")
-        if "delr" in ds:
-            delr = ds["delr"].values
-        else:
-            raise KeyError("delr not in dataset")
         modelgrid = StructuredGrid(
-            delc=delc,
-            delr=delr,
+            delc=get_delc(ds),
+            delr=get_delr(ds),
             **kwargs,
         )
     elif ds.gridtype == "vertex":
@@ -455,8 +449,8 @@ def refine(
             gwf,
             nrow=len(ds.y),
             ncol=len(ds.x),
-            delr=ds.delr,
-            delc=ds.delc,
+            delr=get_delr(ds),
+            delc=get_delc(ds),
             xorigin=ds.extent[0],
             yorigin=ds.extent[2],
         )
@@ -1940,28 +1934,8 @@ def polygons_from_model_ds(model_ds):
     """
 
     if model_ds.gridtype == "structured":
-        # TODO: update with new delr/delc calculation
-        # check if coordinates are consistent with delr/delc values
-        delr_x = np.unique(model_ds.x.values[1:] - model_ds.x.values[:-1])
-        delc_y = np.unique(model_ds.y.values[:-1] - model_ds.y.values[1:])
-
-        delr = np.unique(model_ds.delr)
-        if len(delr) > 1:
-            raise ValueError("delr is variable")
-        else:
-            delr = delr.item()
-
-        delc = np.unique(model_ds.delc)
-        if len(delc) > 1:
-            raise ValueError("delc is variable")
-        else:
-            delc = delc.item()
-
-        if not ((delr_x == delr) and (delc_y == delc)):
-            raise ValueError(
-                "delr and delc attributes of model_ds inconsistent "
-                "with x and y coordinates"
-            )
+        delr = get_delr(model_ds)
+        delc = get_delc(model_ds)
 
         xmins = model_ds.x - (delr * 0.5)
         xmaxs = model_ds.x + (delr * 0.5)
diff --git a/nlmod/dims/resample.py b/nlmod/dims/resample.py
index 73a13427..f3ff7fe1 100644
--- a/nlmod/dims/resample.py
+++ b/nlmod/dims/resample.py
@@ -92,6 +92,58 @@ def get_xy_mid_structured(extent, delr, delc, descending_y=True):
         raise TypeError("unexpected type for delr and/or delc")
 
 
+def get_delr(ds):
+    """
+    Get the distance along rows (delr) from the x-coordinate of a structured model
+    dataset.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        A model dataset containing an x-coordinate and an attribute 'extent'.
+
+    Returns
+    -------
+    delr : np.ndarray
+        The cell-size along rows (of length ncol).
+
+    """
+    assert ds.gridtype == "structured"
+    x = (ds.x - ds.extent[0]).values
+    delr = _get_delta_along_axis(x)
+    return delr
+
+
+def get_delc(ds):
+    """
+    Get the distance along columns (delc) from the y-coordinate of a structured model
+    dataset.
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        A model dataset containing an y-coordinate and an attribute 'extent'.
+
+    Returns
+    -------
+    delc : np.ndarray
+        The cell-size along columns (of length nrow).
+
+    """
+    assert ds.gridtype == "structured"
+    y = (ds.extent[3] - ds.y).values
+    delc = _get_delta_along_axis(y)
+    return delc
+
+
+def _get_delta_along_axis(x):
+    """Internal method to determine delr or delc from x or y relative to xmin or ymax"""
+    delr = [x[0] * 2]
+    for xi in x[1:]:
+        delr.append((xi - np.sum(delr)) * 2)
+    return np.array(delr)
+
+
 def ds_to_structured_grid(
     ds_in,
     extent,
@@ -163,17 +215,11 @@ def ds_to_structured_grid(
             x=x, y=y, method=method, kwargs={"fill_value": "extrapolate"}
         )
         ds_out.attrs = attrs
-        # ds_out = ds_out.expand_dims({"ncol": len(x), "nrow": len(y)})
-        ds_out["delr"] = ("ncol",), delr
-        ds_out["delc"] = ("nrow",), delc
         return ds_out
 
     ds_out = xr.Dataset(coords={"y": y, "x": x, "layer": ds_in.layer.data}, attrs=attrs)
     for var in ds_in.data_vars:
         ds_out[var] = structured_da_to_ds(ds_in[var], ds_out, method=method)
-    # ds_out = ds_out.expand_dims({"ncol": len(x), "nrow": len(y)})
-    ds_out["delr"] = ("x",), delr
-    ds_out["delc"] = ("y",), delc
     return ds_out
 
 
@@ -658,9 +704,13 @@ def get_affine(ds, sx=None, sy=None):
     """Get the affine-transformation, from pixel to real-world coordinates."""
     attrs = _get_attrs(ds)
     if sx is None:
-        sx = attrs["delr"]
+        sx = get_delr(ds)
+        assert len(np.unique(sx)) == 1, "Affine-transformation needs a constant delr"
+        sx = sx[0]
     if sy is None:
-        sy = -attrs["delc"]
+        sy = get_delc(ds)
+        assert len(np.unique(sy)) == 1, "Affine-transformation needs a constant delc"
+        sy = sy[0]
 
     if "angrot" in attrs:
         xorigin = attrs["xorigin"]
diff --git a/nlmod/gis.py b/nlmod/gis.py
index 3998791f..17d4709c 100644
--- a/nlmod/gis.py
+++ b/nlmod/gis.py
@@ -143,9 +143,6 @@ def struc_da_to_gdf(model_ds, data_variables, polygons=None, dealing_with_time="
                 raise ValueError(
                     f"expected dimensions ('layer', 'y', 'x'), got {da.dims}"
                 )
-        # TODO: remove when delr/delc are removed as data vars
-        elif da_name in ["delr", "delc"]:
-            continue
         else:
             raise NotImplementedError(
                 f"expected two or three dimensions got {no_dims} for data variable {da_name}"
@@ -289,9 +286,7 @@ def ds_to_vector_file(
         if model_ds.gridtype == "structured":
             gdf = struc_da_to_gdf(model_ds, (da_name,), polygons=polygons)
         elif model_ds.gridtype == "vertex":
-            # TODO: remove when delr/dec are removed
-            if da_name not in ["delr", "delc"]:
-                gdf = vertex_da_to_gdf(model_ds, (da_name,), polygons=polygons)
+            gdf = vertex_da_to_gdf(model_ds, (da_name,), polygons=polygons)
         if driver == "GPKG":
             gdf.to_file(fname_gpkg, layer=da_name, driver=driver)
         else:
diff --git a/nlmod/gwf/gwf.py b/nlmod/gwf/gwf.py
index 393bb4a7..ce45e29c 100644
--- a/nlmod/gwf/gwf.py
+++ b/nlmod/gwf/gwf.py
@@ -8,10 +8,10 @@
 import warnings
 
 import flopy
-import numpy as np
 import xarray as xr
 
 from ..dims import grid
+from ..dims.resample import get_delr, get_delc
 from ..dims.layers import get_idomain
 from ..sim import ims, sim, tdis
 from ..util import _get_value_from_ds_attr, _get_value_from_ds_datavar
@@ -109,11 +109,6 @@ def _dis(ds, model, length_units="METERS", pname="dis", **kwargs):
         return disv(ds, model, length_units=length_units)
 
     # check attributes
-    for att in ["delr", "delc"]:
-        if att in ds.attrs:
-            if isinstance(ds.attrs[att], np.float32):
-                ds.attrs[att] = float(ds.attrs[att])
-
     if "angrot" in ds.attrs and ds.attrs["angrot"] != 0.0:
         xorigin = ds.attrs["xorigin"]
         yorigin = ds.attrs["yorigin"]
@@ -135,8 +130,8 @@ def _dis(ds, model, length_units="METERS", pname="dis", **kwargs):
             nlay=ds.sizes["layer"],
             nrow=ds.sizes["y"],
             ncol=ds.sizes["x"],
-            delr=ds["delr"].values if "delr" in ds else ds.delr,
-            delc=ds["delc"].values if "delc" in ds else ds.delc,
+            delr=get_delr(ds),
+            delc=get_delc(ds),
             top=ds["top"].data,
             botm=ds["botm"].data,
             idomain=idomain,
@@ -154,8 +149,8 @@ def _dis(ds, model, length_units="METERS", pname="dis", **kwargs):
             nlay=ds.sizes["layer"],
             nrow=ds.sizes["y"],
             ncol=ds.sizes["x"],
-            delr=ds["delr"].values if "delr" in ds else ds.delr,
-            delc=ds["delc"].values if "delc" in ds else ds.delc,
+            delr=get_delr(ds),
+            delc=get_delc(ds),
             top=ds["top"].data,
             botm=ds["botm"].data,
             idomain=idomain,
diff --git a/nlmod/gwf/surface_water.py b/nlmod/gwf/surface_water.py
index 07d7c60c..1cf66719 100644
--- a/nlmod/gwf/surface_water.py
+++ b/nlmod/gwf/surface_water.py
@@ -12,7 +12,7 @@
 
 from ..dims.grid import gdf_to_grid
 from ..dims.layers import get_idomain
-from ..dims.resample import get_extent_polygon, extent_to_polygon
+from ..dims.resample import get_extent_polygon, extent_to_polygon, get_delr, get_delc
 from ..read import bgt, waterboard
 from ..cache import cache_pickle
 
@@ -147,7 +147,11 @@ def agg_de_lange(group, cid, ds, c1=0.0, c0=1.0, N=1e-3, crad_positive=True):
     # correction if group contains multiple shapes
     # but covers whole cell
     if group.area.sum() == A:
-        li = A / np.max([ds.delr, ds.delc])
+        delr = get_delr(ds)
+        assert len(np.unique(delr)) == 1, "Variable grid size is not yet supported"
+        delc = get_delc(ds)
+        assert len(np.unique(delc)) == 1, "Variable grid size is not yet supported"
+        li = A / np.max([delr[0], delc[0]])
 
     # width
     B = group.area.sum(skipna=True) / li
diff --git a/nlmod/read/rws.py b/nlmod/read/rws.py
index 4b00e3c3..bc2f36e3 100644
--- a/nlmod/read/rws.py
+++ b/nlmod/read/rws.py
@@ -37,8 +37,7 @@ def get_gdf_surface_water(ds):
     return gdf_swater
 
 
-# TODO: temporary fix until delr/delc are removed completely
-@cache.cache_netcdf(coords_3d=True, datavars=["delr", "delc"])
+@cache.cache_netcdf(coords_3d=True)
 def get_surface_water(ds, da_basename):
     """create 3 data-arrays from the shapefile with surface water:
 
@@ -92,8 +91,7 @@ def get_surface_water(ds, da_basename):
     return ds_out
 
 
-# TODO: temporary fix until delr/delc are removed completely
-@cache.cache_netcdf(coords_2d=True, datavars=["delc", "delr"])
+@cache.cache_netcdf(coords_2d=True)
 def get_northsea(ds, da_name="northsea"):
     """Get Dataset which is 1 at the northsea and 0 everywhere else. Sea is
     defined by rws surface water shapefile.