diff --git a/nlmod/gwf/surface_water.py b/nlmod/gwf/surface_water.py
index 29d459ac..4a95d2a9 100644
--- a/nlmod/gwf/surface_water.py
+++ b/nlmod/gwf/surface_water.py
@@ -85,6 +85,8 @@ def get_surfacewater_params(group, method, cid=None, ds=None, delange_params=Non
         rbot = group["botm"].min()
 
     elif method == "de_lange":
+        if ds is None:
+            raise ValueError("Please supply model dataset (ds) when method=='de_lange'")
         # get additional requisite parameters
         if delange_params is None:
             delange_params = {}
@@ -120,11 +122,12 @@ def agg_area_weighted(gdf, col):
 
 
 def agg_de_lange(group, cid, ds, c1=0.0, c0=1.0, N=1e-3, crad_positive=True):
-    (A, laytop, laybot, kh, kv, thickness) = get_subsurface_params_by_cellid(ds, cid)
+    (A, laytop, laybot, kh, kv) = get_subsurface_params_by_cellid(ds, cid)
 
     rbot = group["botm"].min()
 
     # select active layers
+    thickness = -np.diff(np.hstack((laytop, laybot)))
     active = thickness > 0
     laybot = laybot[active]
     kh = kh[active]
@@ -132,6 +135,8 @@ def agg_de_lange(group, cid, ds, c1=0.0, c0=1.0, N=1e-3, crad_positive=True):
     thickness = thickness[active]
 
     # layer thickn.
+    if np.isnan(rbot):
+        raise ValueError(f"rbot is NaN in cell {cid}")
     H0 = laytop - laybot[laybot < rbot][0]
     ilay = 0
     rlay = np.where(laybot < rbot)[0][0]
@@ -171,13 +176,12 @@ def agg_de_lange(group, cid, ds, c1=0.0, c0=1.0, N=1e-3, crad_positive=True):
 
 def get_subsurface_params_by_cellid(ds, cid):
     r, c = cid
-    A = ds.delr * ds.delc  # cell area
+    A = ds.area.isel(x=c, y=r).data
     laytop = ds["top"].isel(x=c, y=r).data
-    laybot = ds["bot"].isel(x=c, y=r).data
+    laybot = ds["botm"].isel(x=c, y=r).data
     kv = ds["kv"].isel(x=c, y=r).data
     kh = ds["kh"].isel(x=c, y=r).data
-    thickness = ds["thickness"].isel(x=c, y=r).data
-    return A, laytop, laybot, kh, kv, thickness
+    return A, laytop, laybot, kh, kv
 
 
 def de_lange_eqns(A, H0, kv, kh, c1, li, Bin, c0, p, N, crad_positive=True):
diff --git a/nlmod/read/ahn.py b/nlmod/read/ahn.py
index 8922bf50..4f9d4970 100644
--- a/nlmod/read/ahn.py
+++ b/nlmod/read/ahn.py
@@ -624,6 +624,8 @@ def _get_ahn_ellipsis(extent, identifier="AHN5_5M_M", **kwargs):
     da = merge_arrays(das)
     if da.dims[0] == "band":
         da = da[0].drop_vars("band")
+    if "_FillValue" in da.attrs:
+        del da.attrs["_FillValue"]
     return da
 
 
diff --git a/tests/test_012_plot.py b/tests/test_012_plot.py
index 9eb3279d..3668b244 100644
--- a/tests/test_012_plot.py
+++ b/tests/test_012_plot.py
@@ -8,6 +8,11 @@ def test_plot_modelgrid():
     nlmod.plot.modelgrid(ds)
 
 
+def test_plot_modelextent():
+    ds = util.get_ds_structured()
+    nlmod.plot.modelextent(ds)
+
+
 def test_plot_surface_water_empty():
     ds = util.get_ds_structured()
     nlmod.plot.surface_water(ds)
diff --git a/tests/test_013_surface_water.py b/tests/test_013_surface_water.py
index 8cdc43c7..cfc77507 100644
--- a/tests/test_013_surface_water.py
+++ b/tests/test_013_surface_water.py
@@ -2,28 +2,39 @@
 
 import geopandas as gpd
 import pandas as pd
+import flopy
 
 import nlmod
 
 
-def test_gdf_to_seasonal_pkg():
+def get_ds_and_gdf():
     model_name = "sw"
     model_ws = os.path.join("data", model_name)
     extent = [119000, 120000, 523000, 524000]
     ds = nlmod.get_ds(extent, model_ws=model_ws, model_name=model_name)
     ds = nlmod.time.set_ds_time(ds, time=[365.0], start=pd.Timestamp.today())
-    gdf = nlmod.gwf.surface_water.get_gdf(ds)
+    fname = os.path.join(ds.model_ws, "sw_gdf.gpkg")
+    if not os.path.isfile(fname):
+        gdf = nlmod.gwf.surface_water.get_gdf(ds)
+        gdf.to_file(fname)
+    gdf = gpd.read_file(fname)
+    gdf["cellid"] = [eval(x) for x in gdf["cellid"]]
+    gdf = gdf.set_index("cellid")
+    return ds, gdf
+
+
+def test_gdf_to_seasonal_pkg():
+    ds, gdf = get_ds_and_gdf()
 
     sim = nlmod.sim.sim(ds)
     nlmod.sim.tdis(ds, sim)
-    nlmod.sim.ims(sim)
     gwf = nlmod.gwf.gwf(ds, sim)
     nlmod.gwf.dis(ds, gwf)
-    nlmod.gwf.npf(ds, gwf)
-    nlmod.gwf.ic(ds, gwf, starting_head=1.0)
-    nlmod.gwf.oc(ds, gwf)
 
-    nlmod.gwf.surface_water.gdf_to_seasonal_pkg(gdf, gwf, ds, pkg="DRN")
+    for layer_method in ["lay_of_rbot" and "distribute_cond_over_lays"]:
+        nlmod.gwf.surface_water.gdf_to_seasonal_pkg(
+            gdf, gwf, ds, pkg="DRN", layer_method=layer_method
+        )
 
 
 def test_get_seaonal_timeseries():
@@ -90,3 +101,23 @@ def test_gdf_lake():
     )
 
     nlmod.gwf.lake_from_gdf(gwf, gdf_lake, ds, boundname_column="name")
+
+
+def test_aggregate():
+    ds, gdf = get_ds_and_gdf()
+    gdf = gdf.reset_index()
+    gdf["stage"] = gdf[["summer_stage", "winter_stage"]].mean(1)
+    gdf["botm"] = gdf["bottom_height"]
+    mask = gdf["botm"].isna()
+    gdf.loc[mask, "botm"] = gdf.loc[mask, "stage"] - 0.5
+    gdf["c0"] = 1.0
+
+    sim = nlmod.sim.sim(ds)
+    nlmod.sim.tdis(ds, sim)
+    gwf = nlmod.gwf.gwf(ds, sim)
+    nlmod.gwf.dis(ds, gwf)
+    for method in ["area_weighted", "max_area", "de_lange"]:
+        celldata = nlmod.gwf.aggregate(gdf, method, ds=ds)
+        assert not celldata.isna().any(axis=None)
+        riv_spd = nlmod.gwf.surface_water.build_spd(celldata, "RIV", ds)
+        flopy.mf6.ModflowGwfriv(gwf, stress_period_data=riv_spd)
diff --git a/tests/test_022_gwt.py b/tests/test_022_gwt.py
index 44cf1e16..24cfc2c9 100644
--- a/tests/test_022_gwt.py
+++ b/tests/test_022_gwt.py
@@ -42,7 +42,7 @@ def test_gwt_model():
     ds["sea"] = nlmod.read.rws.calculate_sea_coverage(ahn, ds=ds, method="average")
 
     # download knmi recharge data
-    knmi_ds = nlmod.read.knmi.get_recharge(ds)
+    knmi_ds = nlmod.read.knmi.get_recharge(ds, method="separate")
 
     # update model dataset
     ds.update(knmi_ds)
@@ -90,6 +90,9 @@ def test_gwt_model():
     # create recharge package
     nlmod.gwf.rch(ds, gwf, mask=ds["sea"] == 0)
 
+    # create evapotranspiration package
+    nlmod.gwf.evt(ds, gwf, mask=ds["sea"] == 0)
+
     # BUY: buoyancy package for GWF model
     nlmod.gwf.buy(ds, gwf)