diff --git a/mfsetup/mf6model.py b/mfsetup/mf6model.py
index 90452d05..a552c458 100644
--- a/mfsetup/mf6model.py
+++ b/mfsetup/mf6model.py
@@ -216,7 +216,7 @@ def _set_idomain(self):
         idomain = idomain.astype(int)
 
         # remove cells that conincide with lakes
-        idomain[self.isbc == 1] = 0.
+        # idomain[self.isbc == 1] = 0.
 
         # remove cells that are above stream cells
         if self.get_package('sfr') is not None:
diff --git a/mfsetup/tests/test_lakes.py b/mfsetup/tests/test_lakes.py
index 25dbd0af..dd450db4 100644
--- a/mfsetup/tests/test_lakes.py
+++ b/mfsetup/tests/test_lakes.py
@@ -7,6 +7,7 @@
 import pandas as pd
 import pytest
 
+from mfsetup import MF6model
 from mfsetup.fileio import load_array
 from mfsetup.lakes import (
     PrismSourceData,
@@ -15,6 +16,7 @@
     setup_lake_connectiondata,
     setup_lake_info,
 )
+from mfsetup.tests.test_pleasant_mf6_inset import pleasant_mf6_cfg
 
 
 @pytest.fixture
@@ -232,3 +234,28 @@ def test_get_flux_variable_from_config(modflow_version, config, nlakes, nper, ex
     results = get_flux_variable_from_config(variable, config, nlakes, nper)
     # repeat values for each lake
     assert results == expected
+
+
+def test_lake_idomain(pleasant_mf6_cfg):
+    """For MODFLOW 6 models (at least), if no bathymetry raster is supplied,
+    the lake should be simulated 'above' the grid, connecting vertically with
+    the top model layer, and cells within the lake footprint should only
+    be inactive if they are pinched out based on the original layer surfaces
+    (cells containing lake shouldn't be inactive).
+    """
+    cfg = pleasant_mf6_cfg.copy()
+    del cfg['lak']['source_data']['bathymetry_raster']
+    m = MF6model(cfg=cfg)
+    m.setup_dis()
+    m.setup_tdis()
+    m.setup_lak()
+
+
+    k, i, j = zip(*m.lak.connectiondata.array['cellid'])
+    # only a small number of cells
+    # should have connections to layers below the top layer
+    # note that this 5% criterian might change
+    # if the test problem changes substantially
+    assert np.sum(np.array(k) > 0)/np.sum(np.array(k) == 0) < 0.05
+    # all connections should be to active cells
+    assert all(m.dis.idomain.array[k, i, j] == 1)