diff --git a/autotest/test_netcdf_gwe_cnd.py b/autotest/test_netcdf_gwe_cnd.py
index 95a47210343..f7eaf2759b8 100644
--- a/autotest/test_netcdf_gwe_cnd.py
+++ b/autotest/test_netcdf_gwe_cnd.py
@@ -24,6 +24,7 @@
 
 xa = pytest.importorskip("xarray")
 xu = pytest.importorskip("xugrid")
+netcdf4 = pytest.importorskip("netCDF4")
 
 
 def build_models(idx, test, export, gridded_input):
@@ -57,6 +58,10 @@ def check_output(idx, test, export, gridded_input):
 
     name = "gwe-" + test.name
 
+    # verify format of generated netcdf file
+    with netcdf4.Dataset(test.workspace / f"{name}.nc") as ds:
+        assert ds.data_model == "NETCDF4"
+
     if gridded_input == "netcdf":
         # re-run the simulation with model netcdf input
         input_fname = f"{name}.nc"
diff --git a/autotest/test_netcdf_gwf_disv.py b/autotest/test_netcdf_gwf_disv.py
index 0a0392a5cbd..6d41f5c526a 100644
--- a/autotest/test_netcdf_gwf_disv.py
+++ b/autotest/test_netcdf_gwf_disv.py
@@ -14,6 +14,7 @@
 
 xa = pytest.importorskip("xarray")
 xu = pytest.importorskip("xugrid")
+netcdf4 = pytest.importorskip("netCDF4")
 
 wkt = (
     'PROJCS["NAD83 / UTM zone 18N", '
@@ -72,6 +73,10 @@ def check_output(idx, test, export, gridded_input):
 
     name = test.name
 
+    # verify format of generated netcdf file
+    with netcdf4.Dataset(test.workspace / f"{name}.nc") as ds:
+        assert ds.data_model == "NETCDF4"
+
     if gridded_input == "netcdf":
         # re-run the simulation with model netcdf input
         input_fname = f"{name}.nc"
diff --git a/autotest/test_netcdf_gwf_lak_wetlakbedarea02.py b/autotest/test_netcdf_gwf_lak_wetlakbedarea02.py
index 025cf886cd3..c0641292b03 100644
--- a/autotest/test_netcdf_gwf_lak_wetlakbedarea02.py
+++ b/autotest/test_netcdf_gwf_lak_wetlakbedarea02.py
@@ -24,6 +24,7 @@
 
 xa = pytest.importorskip("xarray")
 xu = pytest.importorskip("xugrid")
+netcdf4 = pytest.importorskip("netCDF4")
 
 
 def build_models(idx, test, export, gridded_input):
@@ -61,6 +62,10 @@ def check_output(idx, test, export, gridded_input):
     name = cases[idx]
     gwfname = "gwf-" + name
 
+    # verify format of generated netcdf file
+    with netcdf4.Dataset(test.workspace / f"{gwfname}.nc") as ds:
+        assert ds.data_model == "NETCDF4"
+
     if gridded_input == "netcdf":
         # re-run the simulation with model netcdf input
         input_fname = f"{gwfname}.nc"
diff --git a/autotest/test_netcdf_gwf_rch01.py b/autotest/test_netcdf_gwf_rch01.py
index c54ab869334..76f02091717 100644
--- a/autotest/test_netcdf_gwf_rch01.py
+++ b/autotest/test_netcdf_gwf_rch01.py
@@ -24,6 +24,7 @@
 
 xa = pytest.importorskip("xarray")
 xu = pytest.importorskip("xugrid")
+netcdf4 = pytest.importorskip("netCDF4")
 
 
 def build_models(idx, test, export, gridded_input):
@@ -58,6 +59,10 @@ def check_output(idx, test, export, gridded_input):
 
     name = "rch"
 
+    # verify format of generated netcdf file
+    with netcdf4.Dataset(test.workspace / f"{name}.nc") as ds:
+        assert ds.data_model == "NETCDF4"
+
     if gridded_input == "netcdf":
         # re-run the simulation with model netcdf input
         input_fname = f"{name}.nc"
diff --git a/autotest/test_netcdf_gwf_rch03.py b/autotest/test_netcdf_gwf_rch03.py
index e473fdc5133..fa91ccece64 100644
--- a/autotest/test_netcdf_gwf_rch03.py
+++ b/autotest/test_netcdf_gwf_rch03.py
@@ -24,6 +24,7 @@
 
 xa = pytest.importorskip("xarray")
 xu = pytest.importorskip("xugrid")
+netcdf4 = pytest.importorskip("netCDF4")
 
 
 def build_models(idx, test, export, gridded_input):
@@ -59,6 +60,10 @@ def check_output(idx, test, export, gridded_input):
 
     name = "rch"
 
+    # verify format of generated netcdf file
+    with netcdf4.Dataset(test.workspace / f"{name}.nc") as ds:
+        assert ds.data_model == "NETCDF4"
+
     if gridded_input == "netcdf":
         # re-run the simulation with model netcdf input
         input_fname = f"{name}.nc"
diff --git a/autotest/test_netcdf_gwf_sto01.py b/autotest/test_netcdf_gwf_sto01.py
index d1939b31452..b591f91abce 100644
--- a/autotest/test_netcdf_gwf_sto01.py
+++ b/autotest/test_netcdf_gwf_sto01.py
@@ -12,6 +12,7 @@
 
 xa = pytest.importorskip("xarray")
 xu = pytest.importorskip("xugrid")
+netcdf4 = pytest.importorskip("netCDF4")
 
 htol = [None for _ in range(len(cases))]
 
@@ -68,6 +69,10 @@ def build_models(idx, test, export, gridded_input):
 def check_output(idx, test, export, gridded_input):
     from test_gwf_sto01 import check_output as check
 
+    # verify format of generated netcdf file
+    with netcdf4.Dataset(test.workspace / "gwf_sto01.nc") as ds:
+        assert ds.data_model == "NETCDF4"
+
     if gridded_input == "netcdf":
         # re-run the simulation with model netcdf input
         input_fname = "gwf_sto01.nc"
diff --git a/autotest/test_netcdf_gwf_vsc03_sfr.py b/autotest/test_netcdf_gwf_vsc03_sfr.py
index 50be2d19098..8ee0d0fe95a 100644
--- a/autotest/test_netcdf_gwf_vsc03_sfr.py
+++ b/autotest/test_netcdf_gwf_vsc03_sfr.py
@@ -24,6 +24,7 @@
 
 xa = pytest.importorskip("xarray")
 xu = pytest.importorskip("xugrid")
+netcdf4 = pytest.importorskip("netCDF4")
 
 
 def build_models(idx, test, export, gridded_input):
@@ -61,6 +62,10 @@ def check_output(idx, test, export, gridded_input):
 
     name = "gwf-" + test.name
 
+    # verify format of generated netcdf file
+    with netcdf4.Dataset(test.workspace / f"{name}.nc") as ds:
+        assert ds.data_model == "NETCDF4"
+
     if gridded_input == "netcdf":
         # re-run the simulation with model netcdf input
         input_fname = f"{name}.nc"
diff --git a/autotest/test_netcdf_gwt_dsp01.py b/autotest/test_netcdf_gwt_dsp01.py
index fd37f1f8c09..410d7c4d6f1 100644
--- a/autotest/test_netcdf_gwt_dsp01.py
+++ b/autotest/test_netcdf_gwt_dsp01.py
@@ -12,6 +12,7 @@
 
 xa = pytest.importorskip("xarray")
 xu = pytest.importorskip("xugrid")
+netcdf4 = pytest.importorskip("netCDF4")
 
 
 def build_models(idx, test, export, gridded_input):
@@ -56,6 +57,10 @@ def check_output(idx, test, export, gridded_input):
     name = cases[idx]
     gwtname = "gwt_" + name
 
+    # verify format of generated netcdf file
+    with netcdf4.Dataset(test.workspace / f"{gwtname}.nc") as ds:
+        assert ds.data_model == "NETCDF4"
+
     if gridded_input == "netcdf":
         # re-run the simulation with model netcdf input
         input_fname = f"{gwtname}.nc"
diff --git a/autotest/test_netcdf_gwt_prudic2004t2.py b/autotest/test_netcdf_gwt_prudic2004t2.py
index a4772c47c79..0f6c2b56a01 100644
--- a/autotest/test_netcdf_gwt_prudic2004t2.py
+++ b/autotest/test_netcdf_gwt_prudic2004t2.py
@@ -12,6 +12,7 @@
 
 xa = pytest.importorskip("xarray")
 xu = pytest.importorskip("xugrid")
+netcdf4 = pytest.importorskip("netCDF4")
 
 
 def build_models(idx, test, export, gridded_input):
@@ -47,6 +48,10 @@ def check_output(idx, test, export, gridded_input):
     name = test.name
     gwtname = "gwt_" + name
 
+    # verify format of generated netcdf file
+    with netcdf4.Dataset(test.workspace / f"{gwtname}.nc") as ds:
+        assert ds.data_model == "NETCDF4"
+
     if gridded_input == "netcdf":
         # re-run the simulation with model netcdf input
         input_fname = f"{gwtname}.nc"
diff --git a/doc/mf6io/mf6ivar/dfn/utl-ncf.dfn b/doc/mf6io/mf6ivar/dfn/utl-ncf.dfn
index cb90c69d1f5..5df020e4ac4 100644
--- a/doc/mf6io/mf6ivar/dfn/utl-ncf.dfn
+++ b/doc/mf6io/mf6ivar/dfn/utl-ncf.dfn
@@ -35,7 +35,7 @@ in_record true
 reader urword
 optional false
 longname chunking parameter for the time dimension
-description is the keyword used to provide a user defined data chunk size for the time dimension in a NETCDF\_MESH2D or NETCDF\_STRUCTURED type netcdf output file. Must be used in combination with the the chunk_face parameter (NETCDF\_MESH2D) or the chunk_z, chunk_y, and chunk_x parameter set (NETCDF\_STRUCTURED) to have an effect.
+description is the keyword used to provide a user defined data chunk size for the time dimension in a NETCDF\_MESH2D or NETCDF\_STRUCTURED type netcdf output file. Must be used in combination with the the chunk\_face parameter (NETCDF\_MESH2D) or the chunk\_z, chunk\_y, and chunk\_x parameter set (NETCDF\_STRUCTURED) to have an effect.
 
 block options
 name chunk_face
@@ -44,7 +44,7 @@ in_record true
 reader urword
 optional true
 longname chunking parameter for the mesh face dimension
-description is the keyword used to provide a user defined data chunk size for the face dimension in a NETCDF\_MESH2D type netcdf output file. Must be used in combination with the the chunk_time parameter to have an effect.
+description is the keyword used to provide a user defined data chunk size for the face dimension in a NETCDF\_MESH2D type netcdf output file. Must be used in combination with the the chunk\_time parameter to have an effect.
 
 block options
 name chunk_z
@@ -53,7 +53,7 @@ in_record true
 reader urword
 optional true
 longname chunking parameter for structured z
-description is the keyword used to provide a user defined data chunk size for the z dimension in a NETCDF\_STRUCTURED type netcdf output file. Must be used in combination with the the chunk_time, chunk_x and chunk_y parameter set to have an effect.
+description is the keyword used to provide a user defined data chunk size for the z dimension in a NETCDF\_STRUCTURED type netcdf output file. Must be used in combination with the the chunk\_time, chunk\_x and chunk\_y parameter set to have an effect.
 
 block options
 name chunk_y
@@ -62,7 +62,7 @@ in_record true
 reader urword
 optional true
 longname chunking parameter for structured y
-description is the keyword used to provide a user defined data chunk size for the y dimension in a NETCDF\_STRUCTURED type netcdf output file. Must be used in combination with the the chunk_time, chunk_x and chunk_z parameter set to have an effect.
+description is the keyword used to provide a user defined data chunk size for the y dimension in a NETCDF\_STRUCTURED type netcdf output file. Must be used in combination with the the chunk\_time, chunk\_x and chunk\_z parameter set to have an effect.
 
 block options
 name chunk_x
@@ -71,7 +71,7 @@ in_record true
 reader urword
 optional true
 longname chunking parameter for structured x
-description is the keyword used to provide a user defined data chunk size for the x dimension in a NETCDF\_STRUCTURED type netcdf output file. Must be used in combination with the the chunk_time, chunk_y and chunk_z parameter set to have an effect.
+description is the keyword used to provide a user defined data chunk size for the x dimension in a NETCDF\_STRUCTURED type netcdf output file. Must be used in combination with the the chunk\_time, chunk\_y and chunk\_z parameter set to have an effect.
 
 block options
 name modflow6_attr_off