diff --git a/.docs/Notebooks/modelgrid_examples.py b/.docs/Notebooks/modelgrid_examples.py
index 6fe6ce6f1e..49dac53b43 100644
--- a/.docs/Notebooks/modelgrid_examples.py
+++ b/.docs/Notebooks/modelgrid_examples.py
@@ -159,7 +159,7 @@
 
 # set reference infromation
 modelgrid1.set_coord_info(
-    xoff=xoff, yoff=yoff, angrot=angrot, epsg=epsg, proj4=proj4
+    xoff=xoff, yoff=yoff, angrot=angrot, crs=epsg
 )
 
 print("After: {}".format(modelgrid1))
@@ -213,9 +213,8 @@
 #    - `botm` : Array of layer Botm elevations
 #    - `idomain` : An ibound or idomain array that specifies active and inactive cells
 #    - `lenuni` : Model length unit integer
-#    - `epsg` : epsg code of model coordinate system
-#    - `proj4` : proj4 str describining model coordinate system
-#    - `prj` : path to ".prj" projection file that describes the model coordinate system
+#    - `crs` : either an epsg integer code of model coordinate system, a proj4 string or pyproj CRS instance
+#    - `prjfile` : path to ".prj" projection file that describes the model coordinate system
 #    - `xoff` : x-coordinate of the lower-left corner of the modelgrid
 #    - `yoff` : y-coordinate of the lower-left corner of the modelgrid
 #    - `angrot` : model grid rotation
@@ -358,9 +357,8 @@
 #    - `botm` : Array of layer Botm elevations
 #    - `idomain` : An ibound or idomain array that specifies active and inactive cells
 #    - `lenuni` : Model length unit integer
-#    - `epsg` : epsg code of model coordinate system
-#    - `proj4` : proj4 str describining model coordinate system
-#    - `prj` : path to ".prj" projection file that describes the model coordinate system
+#    - `crs` : either an epsg integer code of model coordinate system, a proj4 string or pyproj CRS instance
+#    - `prjfile` : path to ".prj" projection file that describes the model coordinate system
 #    - `xoff` : x-coordinate of the lower-left corner of the modelgrid
 #    - `yoff` : y-coordinate of the lower-left corner of the modelgrid
 #    - `angrot` : model grid rotation
@@ -476,9 +474,8 @@ def load_iverts(fname):
 #    - `idomain` : An ibound or idomain array that specifies active and inactive cells
 #    - `lenuni` : Model length unit integer
 #    - `ncpl` : one dimensional array of number of cells per model layer
-#    - `epsg` : epsg code of model coordinate system
-#    - `proj4` : proj4 str describining model coordinate system
-#    - `prj` : path to ".prj" projection file that describes the model coordinate system
+#    - `crs` : either an epsg integer code of model coordinate system, a proj4 string or pyproj CRS instance
+#    - `prjfile` : path to ".prj" projection file that describes the model coordinate system
 #    - `xoff` : x-coordinate of the lower-left corner of the modelgrid
 #    - `yoff` : y-coordinate of the lower-left corner of the modelgrid
 #    - `angrot` : model grid rotation   
@@ -579,7 +576,7 @@ def load_iverts(fname):
     xoff=622241.1904510253,
     yoff=3343617.741737109,
     angrot=15.0,
-    proj4="+proj=utm +zone=14 +ellps=WGS84 +datum=WGS84 +units=m +no_defs",
+    crs="+proj=utm +zone=14 +ellps=WGS84 +datum=WGS84 +units=m +no_defs",
 )
 # -
 
@@ -620,7 +617,7 @@ def load_iverts(fname):
 #    - `angrot_radians` : returns the angle of rotation of the modelgrid in radians
 #    - `epsg` : returns the modelgrid epsg code if it is set
 #    - `proj4` : returns the modelgrid proj4 string if it is set
-#    - `prj` : returns the path to the modelgrid projection file if it is set
+#    - `prjfile` : returns the path to the modelgrid projection file if it is set
 
 # Access and print some of these properties
 print(
@@ -723,8 +720,8 @@ def load_iverts(fname):
 #    - `xoff` : lower-left corner of modelgrid x-coordinate location
 #    - `yoff` : lower-left corner of modelgrid y-coordinate location
 #    - `angrot` : rotation of model grid in degrees
-#    - `epsg` : epsg code for model grid projection
-#    - `proj4` : proj4 string describing the model grid projection
+#    - `crs` : either an epsg integer code of model coordinate system, a proj4 string or pyproj CRS instance
+#    - `prjfile` : path to ".prj" projection file that describes the model coordinate system
 #    - `merge_coord_info` : boolean flag to either merge changes with the existing coordinate info or clear existing coordinate info before applying changes.
 
 # +
@@ -900,15 +897,15 @@ def load_iverts(fname):
 # Method to write a shapefile of the grid with just the cellid attributes. Input parameters include:
 #
 #    - `filename` : shapefile name
-#    - `epsg` : optional epsg code of the coordinate system projection
-#    - `prj` : optional, input projection file to be used to define the coordinate system projection
+#    - `crs` : either an epsg integer code of model coordinate system, a proj4 string or pyproj CRS instance
+#    - `prjfile` : path to ".prj" projection file that describes the model coordinate system
 
 # +
 # write a shapefile
 shp_name = os.path.join(gridgen_ws, "freyberg-6_grid.shp")
-epsg = "32614"
+epsg = 32614
 
-ml.modelgrid.write_shapefile(shp_name, epsg)
+ml.modelgrid.write_shapefile(shp_name, crs=epsg)
 # -
 
 try:
diff --git a/.docs/Notebooks/shapefile_export_example.py b/.docs/Notebooks/shapefile_export_example.py
index e790d64b75..ac6dc87e10 100644
--- a/.docs/Notebooks/shapefile_export_example.py
+++ b/.docs/Notebooks/shapefile_export_example.py
@@ -63,7 +63,7 @@
 # the coordinate information where the grid is located in a projected coordinate system (e.g. UTM)
 
 grid = m.modelgrid
-grid.set_coord_info(xoff=273170, yoff=5088657, epsg=26916)
+grid.set_coord_info(xoff=273170, yoff=5088657, crs=26916)
 
 grid.extent
 
@@ -142,7 +142,7 @@
 # ##### write the shapefile
 
 fname = "{}/bcs.shp".format(outdir)
-recarray2shp(spd.to_records(), geoms=polygons, shpname=fname, epsg=grid.epsg)
+recarray2shp(spd.to_records(), geoms=polygons, shpname=fname, crs=grid.epsg)
 
 ax = plt.subplot(1, 1, 1, aspect="equal")
 extents = grid.extent
@@ -172,7 +172,7 @@
 geoms = [Point(x, y) for x, y in zip(welldata.x_utm, welldata.y_utm)]
 
 fname = "{}/wel_data.shp".format(outdir)
-recarray2shp(welldata.to_records(), geoms=geoms, shpname=fname, epsg=grid.epsg)
+recarray2shp(welldata.to_records(), geoms=geoms, shpname=fname, crs=grid.epsg)
 # -
 
 ax = plt.subplot(1, 1, 1, aspect="equal")
@@ -206,7 +206,7 @@
     rivdata.to_records(index=False),
     geoms=lines,
     shpname=lines_shapefile,
-    epsg=grid.epsg,
+    crs=grid.epsg,
 )
 # -
 
@@ -241,7 +241,7 @@
     linesdata.drop("geometry", axis=1).to_records(),
     geoms=linesdata.geometry.values,
     shpname=lines_shapefile,
-    epsg=grid.epsg,
+    crs=grid.epsg,
 )
 
 ax = plt.subplot(1, 1, 1, aspect="equal")
diff --git a/autotest/test_export.py b/autotest/test_export.py
index 96f3503a28..a21c060c11 100644
--- a/autotest/test_export.py
+++ b/autotest/test_export.py
@@ -328,8 +328,12 @@ def test_write_gridlines_shapefile(function_tmpdir):
     outshp = function_tmpdir / "gridlines.shp"
     write_gridlines_shapefile(outshp, sg)
 
-    for suffix in [".dbf", ".prj", ".shp", ".shx"]:
+    for suffix in [".dbf", ".shp", ".shx"]:
         assert outshp.with_suffix(suffix).exists()
+    if has_pkg("pyproj"):
+        assert outshp.with_suffix(".prj").exists()
+    else:
+        assert not outshp.with_suffix(".prj").exists()
 
     with shapefile.Reader(str(outshp)) as sf:
         assert sf.shapeType == shapefile.POLYLINE
@@ -991,7 +995,7 @@ def test_polygon_from_ij_with_epsg(function_tmpdir):
     fpth2 = os.path.join(ws, "26715.prj")
     shutil.copy(fpth, fpth2)
     fpth = os.path.join(ws, "test.shp")
-    recarray2shp(recarray, geoms, fpth, prj=fpth2)
+    recarray2shp(recarray, geoms, fpth, prjfile=fpth2)
 
     # test_dtypes
     fpth = os.path.join(ws, "test.shp")
diff --git a/autotest/test_grid.py b/autotest/test_grid.py
index 22b9724a80..246f9e48b9 100644
--- a/autotest/test_grid.py
+++ b/autotest/test_grid.py
@@ -1,4 +1,6 @@
+import re
 import os
+import warnings
 from warnings import warn
 
 import matplotlib
@@ -9,10 +11,9 @@
 from flaky import flaky
 from matplotlib import pyplot as plt
 from modflow_devtools.markers import requires_exe, requires_pkg
-from packaging import version
+from modflow_devtools.misc import has_pkg
 from pytest_cases import parametrize_with_cases
 
-import flopy
 from flopy.discretization import StructuredGrid, UnstructuredGrid, VertexGrid
 from flopy.mf6 import MFSimulation
 from flopy.modflow import Modflow, ModflowDis
@@ -22,6 +23,11 @@
 from flopy.utils.voronoi import VoronoiGrid
 
 
+HAS_PYPROJ = has_pkg("pyproj")
+if HAS_PYPROJ:
+    import pyproj
+
+
 @pytest.fixture
 def minimal_unstructured_grid_info():
     d = {
@@ -549,7 +555,6 @@ def test_unstructured_from_gridspec(example_data_path):
         assert min(grid.botm) == min([xyz[2] for xyz in expected_verts])
 
 
-@requires_pkg("pyproj")
 @pytest.mark.parametrize(
     "crs,expected_srs",
     (
@@ -566,43 +571,59 @@ def test_unstructured_from_gridspec(example_data_path):
 def test_grid_crs(
     minimal_unstructured_grid_info, crs, expected_srs, function_tmpdir
 ):
-    import pyproj
+    expected_epsg = None
+    if match := re.findall(r"epsg:([\d]+)", expected_srs or "", re.IGNORECASE):
+        expected_epsg = int(match[0])
+    if not HAS_PYPROJ and isinstance(crs, str) and "epsg" not in crs.lower():
+        # pyproj needed to derive 'epsg' from PROJ string
+        expected_epsg = None
 
     d = minimal_unstructured_grid_info
     delr = np.ones(10)
     delc = np.ones(10)
-    sg = StructuredGrid(delr=delr, delc=delc, crs=crs)
-    if crs is not None:
-        assert isinstance(sg.crs, pyproj.CRS)
-        assert sg.crs.srs == expected_srs
 
-    usg = UnstructuredGrid(**d, crs=crs)
-    assert getattr(sg.crs, "srs", None) == expected_srs
+    def do_checks(g):
+        if HAS_PYPROJ and crs is not None:
+            assert isinstance(g.crs, pyproj.CRS)
+            assert g.crs.srs == expected_srs
+        else:
+            assert g.crs is None
+        assert g.epsg == expected_epsg
+
+    # test each grid type (class)
+    do_checks(StructuredGrid(delr=delr, delc=delc, crs=crs))
+    do_checks(UnstructuredGrid(**d, crs=crs))
+    do_checks(VertexGrid(vertices=d["vertices"], crs=crs))
 
-    vg = VertexGrid(vertices=d["vertices"], crs=crs)
-    assert getattr(sg.crs, "srs", None) == expected_srs
+    # test deprecated 'epsg' parameter
+    if isinstance(crs, int):
+        with pytest.deprecated_call():
+            do_checks(StructuredGrid(delr=delr, delc=delc, epsg=crs))
 
-    # test input of pyproj.CRS object
-    if crs == 26916:
-        sg2 = StructuredGrid(delr=delr, delc=delc, crs=sg.crs)
+    if HAS_PYPROJ and crs == 26916:
+        crs_obj = get_authority_crs(crs)
 
-        if crs is not None:
-            assert isinstance(sg2.crs, pyproj.CRS)
-        assert getattr(sg2.crs, "srs", None) == expected_srs
+        # test input of pyproj.CRS object
+        do_checks(StructuredGrid(delr=delr, delc=delc, crs=crs_obj))
 
         # test input of projection file
         prjfile = function_tmpdir / "grid_crs.prj"
-        with open(prjfile, "w", encoding="utf-8") as dest:
-            write_text = sg.crs.to_wkt()
-            dest.write(write_text)
+        prjfile.write_text(crs_obj.to_wkt(), encoding="utf-8")
+
+        do_checks(StructuredGrid(delr=delr, delc=delc, prjfile=prjfile))
 
-        sg3 = StructuredGrid(delr=delr, delc=delc, prjfile=prjfile)
-        if crs is not None:
-            assert isinstance(sg3.crs, pyproj.CRS)
-        assert getattr(sg3.crs, "srs", None) == expected_srs
+        # test deprecated 'prj' parameter
+        with pytest.deprecated_call():
+            do_checks(StructuredGrid(delr=delr, delc=delc, prj=prjfile))
+
+        # test deprecated 'proj4' parameter
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")  # pyproj warning about conversion
+            proj4 = crs_obj.to_proj4()
+        with pytest.deprecated_call():
+            do_checks(StructuredGrid(delr=delr, delc=delc, proj4=proj4))
 
 
-@requires_pkg("pyproj")
 @pytest.mark.parametrize(
     "crs,expected_srs",
     (
@@ -618,82 +639,208 @@ def test_grid_crs(
     ),
 )
 def test_grid_set_crs(crs, expected_srs, function_tmpdir):
-    import pyproj
+    expected_epsg = None
+    if match := re.findall(r"epsg:([\d]+)", expected_srs or "", re.IGNORECASE):
+        expected_epsg = int(match[0])
+    if not HAS_PYPROJ and isinstance(crs, str) and "epsg" not in crs.lower():
+        # pyproj needed to derive 'epsg' from PROJ string
+        expected_epsg = None
+    elif HAS_PYPROJ and crs is not None and expected_epsg is None:
+        expected_epsg = pyproj.CRS.from_user_input(crs).to_epsg()
 
     delr = np.ones(10)
     delc = np.ones(10)
-    sg = StructuredGrid(delr=delr, delc=delc)
-    sg.crs = crs
-    if crs is not None:
-        assert isinstance(sg.crs, pyproj.CRS)
-    assert getattr(sg.crs, "srs", None) == expected_srs
 
-    # test setting the crs via set_coord_info
+    def do_checks(g, *, exp_srs=expected_srs, exp_epsg=expected_epsg):
+        if HAS_PYPROJ:
+            if crs is not None:
+                assert isinstance(g.crs, pyproj.CRS)
+            assert getattr(g.crs, "srs", None) == exp_srs
+        else:
+            assert g.crs is None
+        assert g.epsg == exp_epsg
+
+    # test set_coord_info with a grid object
+    sg = StructuredGrid(delr=delr, delc=delc, crs=crs)
+    do_checks(sg)
+
+    # no change
     sg.set_coord_info()
-    if crs is not None:
-        assert isinstance(sg.crs, pyproj.CRS)
-    assert getattr(sg.crs, "srs", None) == expected_srs
+    do_checks(sg)
+
+    # use 'crs' arg
     sg.set_coord_info(crs=crs)
-    if crs is not None:
-        assert isinstance(sg.crs, pyproj.CRS)
-    assert getattr(sg.crs, "srs", None) == expected_srs
+    do_checks(sg)
+
+    # use different 'crs'
     sg.set_coord_info(crs=26915, merge_coord_info=False)
-    assert getattr(sg.crs, "srs", None) == "EPSG:26915"
-    # reset back to test case crs
+    do_checks(sg, exp_srs="EPSG:26915", exp_epsg=26915)
+
+    # test deprecated 'epsg' parameter
+    if isinstance(crs, int):
+        with pytest.deprecated_call():
+            sg.set_coord_info(epsg=crs)
+        do_checks(sg)
+
+    # use 'crs' setter
+    sg = StructuredGrid(delr=delr, delc=delc)
+    if HAS_PYPROJ:
+        sg.crs = crs
+        do_checks(sg)
+    else:
+        if crs is None:
+            sg.crs = crs
+        else:
+            with pytest.warns():
+                # cannot set 'crs' property without pyproj
+                sg.crs = crs
+        do_checks(sg, exp_epsg=None)
+
+    # unset 'crs', and check that 'epsg' is also none (sometimes)
+    sg = StructuredGrid(delr=delr, delc=delc, crs=crs)
+    sg.crs = None
+    assert sg.crs is None
+    if HAS_PYPROJ:
+        # with pyproj, '_epsg' was never populated by getter
+        assert sg.epsg is None
+    else:
+        # without pyproj, '_epsg' is populated from specific 'crs' inputs
+        assert sg.epsg == expected_epsg
+
+    # unset 'crs', but check that 'epsg' is retained
     sg = StructuredGrid(delr=delr, delc=delc, crs=crs)
+    assert sg.epsg == expected_epsg  # populate '_epsg' via getter
+    sg.crs = None
+    assert sg.crs is None
+    assert sg.epsg == expected_epsg
+
+    # unset 'epsg'
+    sg.epsg = None
+    assert sg.epsg is None
+
+    # unset 'proj4'
+    sg.proj4 = None
+    assert sg.proj4 is None
+
+    # set and unset 'prjfile' with a non-existing file
+    prjfile = "test"
+    assert sg.prjfile is None
+    sg.prjfile = prjfile
+    assert sg.prjfile == prjfile
+    sg.prjfile = None
+    assert sg.prjfile is None
+
+    if HAS_PYPROJ and crs is not None:
+        crs_obj = get_authority_crs(crs)
 
-    # test input of projection file
-    if crs is not None:
+        # test input of projection file
         prjfile = function_tmpdir / "grid_crs.prj"
-        with open(prjfile, "w", encoding="utf-8") as dest:
-            write_text = sg.crs.to_wkt()
-            dest.write(write_text)
+        prjfile.write_text(crs_obj.to_wkt(), encoding="utf-8")
+
+        def do_prjfile_checks(g):
+            assert isinstance(g.crs, pyproj.CRS)
+            assert g.crs.srs == expected_srs
+            assert g.epsg == expected_epsg
+            assert g.prjfile == prjfile
+
+        # test with 'prjfile' setter and parameter
         sg = StructuredGrid(delr=delr, delc=delc)
         sg.prjfile = prjfile
-        if crs is not None:
-            assert isinstance(sg.crs, pyproj.CRS)
-        assert getattr(sg.crs, "srs", None) == expected_srs
-        assert sg.prjfile == prjfile
+        do_prjfile_checks(sg)
         sg.set_coord_info(prjfile=prjfile)
-        if crs is not None:
-            assert isinstance(sg.crs, pyproj.CRS)
-            assert getattr(sg.crs, "srs", None) == expected_srs
-        assert sg.prjfile == prjfile
-
-    # test setting another crs
-    sg.crs = 26915
-    assert sg.crs == get_authority_crs(26915)
-
-    if (
-        version.parse(flopy.__version__) < version.parse("3.4")
-        and crs is not None
-    ):
-        pyproj_crs = get_authority_crs(crs)
-        epsg = pyproj_crs.to_epsg()
-        if epsg is not None:
+        do_prjfile_checks(sg)
+
+        # test with deprecated 'prj' getter/setter
+        with pytest.deprecated_call():
+            assert sg.prj == prjfile
+        with pytest.deprecated_call():
+            sg.prj = prjfile
+        do_prjfile_checks(sg)
+
+        # test deprecated 'proj4' parameter
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")  # pyproj warning about conversion
+            proj4 = crs_obj.to_proj4()
+        with pytest.deprecated_call():
+            sg.set_coord_info(proj4=proj4)
+        do_checks(sg)
+
+    if HAS_PYPROJ:
+        # copy previous non-None epsg
+        prev_epsg = sg.epsg
+        # test setting another crs
+        sg.crs = 26915
+        assert sg.crs == get_authority_crs(26915)
+        # note that 'epsg' is not updated by setting 'crs', unless it was None
+        assert sg.epsg == prev_epsg or 26915
+
+    if HAS_PYPROJ and crs is not None:
+        if epsg := crs_obj.to_epsg():
             sg = StructuredGrid(delr=delr, delc=delc, crs=crs)
             sg.epsg = epsg
-            if crs is not None:
-                assert isinstance(sg.crs, pyproj.CRS)
-            assert getattr(sg.crs, "srs", None) == expected_srs
-            assert sg.epsg == epsg
+            do_checks(sg, exp_epsg=epsg)
 
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")  # pyproj warning about conversion
+            proj4 = crs_obj.to_proj4()
         sg = StructuredGrid(delr=delr, delc=delc)
-        sg.proj4 = pyproj_crs.to_proj4()
-        if crs is not None:
-            assert isinstance(sg.crs, pyproj.CRS)
-        assert getattr(sg.crs, "srs", None) == expected_srs
-        assert sg.proj4 == pyproj_crs.to_proj4()
-
-        if crs is not None:
-            sg = StructuredGrid(delr=delr, delc=delc)
+        sg.proj4 = proj4
+        do_checks(sg)
+        assert sg.proj4 == proj4
+
+        sg = StructuredGrid(delr=delr, delc=delc)
+        with pytest.deprecated_call():
             sg.prj = prjfile
-            if crs is not None:
-                assert isinstance(sg.crs, pyproj.CRS)
-            assert getattr(sg.crs, "srs", None) == expected_srs
+        do_checks(sg)
+        with pytest.deprecated_call():
             assert sg.prj == prjfile
 
 
+def test_grid_crs_exceptions():
+    delr = np.ones(10)
+    delc = np.ones(10)
+    sg = StructuredGrid(delr=delr, delc=delc, crs="EPSG:26915")
+
+    # test bad 'epsg' parameter
+    bad_epsg = "EPSG:26915"
+    with pytest.raises(ValueError):
+        StructuredGrid(delr=delr, delc=delc, epsg=bad_epsg)
+    with pytest.raises(ValueError):
+        sg.epsg = bad_epsg
+    with pytest.raises(ValueError):
+        sg.set_coord_info(epsg=bad_epsg)
+
+    # test bad 'proj4' parameter
+    bad_proj4 = 26915
+    with pytest.raises(ValueError):
+        StructuredGrid(delr=delr, delc=delc, proj4=bad_proj4)
+    with pytest.raises(ValueError):
+        sg.proj4 = bad_proj4
+    with pytest.raises(ValueError):
+        sg.set_coord_info(proj4=bad_proj4)
+
+    # test bad 'prjfile' parameter
+    bad_prjfile = 0
+    with pytest.raises(ValueError):
+        StructuredGrid(delr=delr, delc=delc, prjfile=bad_prjfile)
+    with pytest.raises(ValueError):
+        sg.prjfile = bad_prjfile
+
+    # test non-existing file
+    not_a_file = "not-a-file"
+    with pytest.raises(FileNotFoundError):
+        StructuredGrid(delr=delr, delc=delc, prjfile=not_a_file)
+    # note "sg.prjfile = not_a_file" intentionally does not raise anything
+
+    # test unhandled keyword
+    with pytest.raises(TypeError):
+        StructuredGrid(delr=delr, delc=delc, unused_param=None)
+
+    # set_coord_info never had a 'prj' parameter; test it as unhandled keyword
+    with pytest.raises(TypeError):
+        sg.set_coord_info(prj=not_a_file)
+
+
 def test_tocvfd1():
     vertdict = {}
     vertdict[0] = [(0, 0), (100, 0), (100, 100), (0, 100), (0, 0)]
diff --git a/autotest/test_lgr.py b/autotest/test_lgr.py
index 2c72b68740..2a01b8ae8c 100644
--- a/autotest/test_lgr.py
+++ b/autotest/test_lgr.py
@@ -33,7 +33,7 @@ def singleModel(
     steady,
     xul,
     yul,
-    proj4_str,
+    crs,
     mfExe,
     rundir=".",
     welInfo=[],
@@ -77,7 +77,7 @@ def singleModel(
         unitnumber=11 + iLUoffset,
         xul=xul,
         yul=yul,
-        proj4_str=proj4_str,
+        crs=crs,
         start_datetime="28/2/2019",
     )
 
@@ -192,7 +192,7 @@ def test_simple_lgrmodel_from_scratch(function_tmpdir):
     laytyp = 0
     xul_c = 50985.00
     yul_c = 416791.06
-    proj4_str = "EPSG:28992"
+    crs = "EPSG:28992"
     nper = 1
     at = 42
     perlen = [at]
@@ -235,7 +235,7 @@ def test_simple_lgrmodel_from_scratch(function_tmpdir):
         steady,
         xul_c,
         yul_c,
-        proj4_str,
+        crs,
         "mflgr",
         rundir=function_tmpdir,
         welInfo=welInfo,
@@ -265,7 +265,7 @@ def test_simple_lgrmodel_from_scratch(function_tmpdir):
         steady,
         xul_m,
         yul_m,
-        proj4_str,
+        crs,
         "mflgr",
         rundir=function_tmpdir,
         welInfo=welInfo,
diff --git a/autotest/test_modflow.py b/autotest/test_modflow.py
index 72993c8909..3838ef28cf 100644
--- a/autotest/test_modflow.py
+++ b/autotest/test_modflow.py
@@ -261,8 +261,7 @@ def test_sr(function_tmpdir):
         raise AssertionError()
     if extents[3] != 12355:
         raise AssertionError()
-    if mm.modelgrid.crs.srs != "EPSG:26916":
-        raise AssertionError()
+    assert mm.modelgrid.epsg == 26916
 
     mm.dis.top = 5000
 
@@ -543,7 +542,7 @@ def test_read_usgs_model_reference(function_tmpdir, model_reference_path):
     m2 = Modflow.load("junk.nam", model_ws=ws)
     m2.modelgrid.read_usgs_model_reference_file(mrf_path)
 
-    assert m2.modelgrid.crs.to_epsg() == 26916
+    assert m2.modelgrid.epsg == 26916
     # have to delete this, otherwise it will mess up other tests
     to_del = glob.glob(f"{mrf_path}*")
     for f in to_del:
diff --git a/autotest/test_modflowdis.py b/autotest/test_modflowdis.py
index eb36dbe3e5..be51b7c014 100644
--- a/autotest/test_modflowdis.py
+++ b/autotest/test_modflowdis.py
@@ -25,10 +25,12 @@ def test_dis_sr():
         rotation=rotation,
         xul=xul,
         yul=yul,
-        proj4_str="epsg:2243",
+        crs="epsg:2243",
     )
 
     # Use StructuredGrid instead
     x, y = bg.modelgrid.get_coords(0, delc * nrow)
     np.testing.assert_almost_equal(x, xul)
     np.testing.assert_almost_equal(y, yul)
+    assert bg.modelgrid.epsg == 2243
+    assert bg.modelgrid.angrot == rotation
diff --git a/autotest/test_mp6.py b/autotest/test_mp6.py
index a42a07f7f4..4637876e01 100644
--- a/autotest/test_mp6.py
+++ b/autotest/test_mp6.py
@@ -257,8 +257,7 @@ def test_get_destination_data(function_tmpdir, mp6_test_path):
         xoff=mg.xoffset,
         yoff=mg.yoffset,
         angrot=mg.angrot,
-        epsg=mg.epsg,
-        proj4=mg.proj4,
+        crs=mg.epsg,
     )
     ra = shp2recarray(function_tmpdir / "pathlines_1per2.shp")
     p3_2 = ra.geometry[ra.particleid == 4][0]
@@ -298,8 +297,7 @@ def test_get_destination_data(function_tmpdir, mp6_test_path):
         xoff=mg4._xul_to_xll(xul, 0.0),
         yoff=mg4._yul_to_yll(yul, 0.0),
         angrot=0.0,
-        epsg=mg4.epsg,
-        proj4=mg4.proj4,
+        crs=mg4.epsg,
     )
 
     fpth = function_tmpdir / "dis2.shp"
diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py
index 428873846d..27ed0ee3e4 100644
--- a/flopy/discretization/grid.py
+++ b/flopy/discretization/grid.py
@@ -1,5 +1,6 @@
 import copy
 import os
+import re
 import warnings
 from collections import defaultdict
 
@@ -28,6 +29,18 @@ def update_data(self, data):
         self.out_of_date = False
 
 
+def _get_epsg_from_crs_or_proj4(crs, proj4=None):
+    """Try to get EPSG identifier from a crs object."""
+    if isinstance(crs, int):
+        return crs
+    if isinstance(crs, str):
+        if match := re.findall(r"epsg:([\d]+)", crs, re.IGNORECASE):
+            return int(match[0])
+    if proj4 and isinstance(proj4, str):
+        if match := re.findall(r"epsg:([\d]+)", proj4, re.IGNORECASE):
+            return int(match[0])
+
+
 class Grid:
     """
     Base class for a structured or unstructured model grid
@@ -44,7 +57,7 @@ class Grid:
         ibound/idomain value for each cell
     lenuni : int or ndarray
         model length units
-    crs : pyproj.CRS, optional if `prjfile` is specified
+    crs : pyproj.CRS, int, str, optional if `prjfile` is specified
         Coordinate reference system (CRS) for the model grid
         (must be projected; geographic CRS are not supported).
         The value can be anything accepted by
@@ -61,6 +74,15 @@ class Grid:
         in the spatial reference coordinate system
     angrot : float
         rotation angle of model grid, as it is rotated around the origin point
+    **kwargs : dict, optional
+        Support deprecated keyword options.
+
+        .. deprecated:: 3.3.7
+           The following keyword options will be removed for FloPy 3.4:
+
+             - ``prj`` (str or pathlike): use ``prjfile`` instead.
+             - ``epsg`` (int): use ``crs`` instead.
+             - ``proj4`` (str): use ``crs`` instead.
 
     Attributes
     ----------
@@ -72,8 +94,6 @@ class Grid:
         bottom elevations of all cells
     idomain : int or ndarray
         ibound/idomain value for each cell
-    crs : pyproj.CRS
-        Coordinate reference system (CRS) for the model grid
     lenuni : int
         modflow lenuni parameter
     xoffset : float
@@ -148,13 +168,11 @@ def __init__(
         idomain=None,
         lenuni=None,
         crs=None,
-        epsg=None,
-        proj4=None,
-        prj=None,
         prjfile=None,
         xoff=0.0,
         yoff=0.0,
         angrot=0.0,
+        **kwargs,
     ):
         lenunits = {0: "undefined", 1: "feet", 2: "meters", 3: "centimeters"}
         LENUNI = {"u": 0, "f": 1, "m": 2, "c": 3}
@@ -175,12 +193,28 @@ def __init__(
         self._lenuni = lenuni
 
         self._units = lenunits[self._lenuni]
-        self._crs = get_crs(
-            prjfile=prjfile, prj=prj, epsg=epsg, proj4=proj4, crs=crs
-        )
-        self._epsg = epsg
-        self._proj4 = proj4
-        self._prj = prj
+        # Handle deprecated projection kwargs; warnings are raised in crs.py
+        self._crs = None
+        get_crs_args = {"crs": crs, "prjfile": prjfile}
+        if "epsg" in kwargs:
+            self.epsg = get_crs_args["epsg"] = kwargs.pop("epsg")
+        if "proj4" in kwargs:
+            self.proj4 = get_crs_args["proj4"] = kwargs.pop("proj4")
+        if "prj" in kwargs:
+            self.prjfile = get_crs_args["prj"] = kwargs.pop("prj")
+        if kwargs:
+            raise TypeError(f"unhandled keywords: {kwargs}")
+        if prjfile is not None:
+            self.prjfile = prjfile
+        try:
+            self._crs = get_crs(**get_crs_args)
+        except ImportError:
+            # provide some support without pyproj by retaining 'epsg' integer
+            if getattr(self, "_epsg", None) is None:
+                epsg = _get_epsg_from_crs_or_proj4(crs, self.proj4)
+                if epsg is not None:
+                    self.epsg = epsg
+
         self._prjfile = prjfile
         self._xoff = xoff
         self._yoff = yoff
@@ -214,6 +248,10 @@ def __repr__(self):
             ]
         if self.crs is not None:
             items.append(f"crs:{self.crs.srs}")
+        elif self.epsg is not None:
+            items.append(f"crs:EPSG:{self.epsg}")
+        elif self.proj4 is not None:
+            items.append(f"proj4_str:{self.proj4}")
         if self.units is not None:
             items.append(f"units:{self.units}")
         if self.lenuni is not None:
@@ -256,49 +294,123 @@ def angrot_radians(self):
 
     @property
     def crs(self):
+        """Coordinate reference system (CRS) for the model grid.
+
+        If pyproj is not installed this property is always None;
+        see :py:attr:`epsg` for an alternative CRS definition.
+        """
         return self._crs
 
     @crs.setter
     def crs(self, crs):
-        self._crs = get_crs(crs=crs)
+        if crs is None:
+            self._crs = None
+            return
+        try:
+            self._crs = get_crs(crs=crs)
+        except ImportError:
+            warnings.warn(
+                "cannot set 'crs' property without pyproj; "
+                "try setting 'epsg' or 'proj4' instead",
+                UserWarning,
+            )
+            self._crs = None
 
     @property
     def epsg(self):
-        if self._crs is not None:
-            return self._crs.to_epsg()
+        """EPSG integer code registered to a coordinate reference system.
+
+        This property is derived from :py:attr:`crs` if pyproj is installed,
+        otherwise it preserved from the constructor.
+        """
+        epsg = None
+        if hasattr(self, "_epsg"):
+            epsg = self._epsg
+        elif self._crs is not None:
+            epsg = self._crs.to_epsg()
+            if epsg is not None:
+                self._epsg = epsg
+        return epsg
 
     @epsg.setter
     def epsg(self, epsg):
-        self._crs = get_crs(epsg=epsg)
-        self._epsg = self._crs.to_epsg()
+        print(f"epsg.setter: {epsg}")
+        if not (isinstance(epsg, int) or epsg is None):
+            raise ValueError("epsg property must be an int or None")
+        self._epsg = epsg
+        # If crs was previously unset, use EPSG code
+        if self._crs is None and epsg is not None:
+            try:
+                self._crs = get_crs(crs=epsg)
+            except ImportError:
+                pass
 
     @property
     def proj4(self):
-        if self._crs is not None:
-            return self._crs.to_proj4()
+        """PROJ string for a coordinate reference system.
+
+        This property is derived from :py:attr:`crs` if pyproj is installed,
+        otherwise it preserved from the constructor.
+        """
+        proj4 = None
+        if hasattr(self, "_proj4"):
+            proj4 = self._proj4
+        elif self._crs is not None:
+            proj4 = self._crs.to_proj4()
+            if proj4 is not None:
+                self._proj4 = proj4
+        return proj4
 
     @proj4.setter
     def proj4(self, proj4):
-        self._crs = get_crs(proj4=proj4)
-        self._proj4 = self._crs.to_proj4()
+        print(f"proj4.setter: {proj4}")
+        if not (isinstance(proj4, str) or proj4 is None):
+            raise ValueError("proj4 property must be a str or None")
+        self._proj4 = proj4
+        # If crs was previously unset, use lossy PROJ string
+        if self._crs is None and proj4 is not None:
+            try:
+                self._crs = get_crs(crs=proj4)
+            except ImportError:
+                pass
 
     @property
     def prj(self):
-        return self._prj
+        warnings.warn(
+            "prj property is deprecated, use prjfile instead",
+            PendingDeprecationWarning,
+        )
+        return self.prjfile
 
     @prj.setter
     def prj(self, prj):
-        self._crs = get_crs(prj=prj)
-        self._prj = prj
+        warnings.warn(
+            "prj property is deprecated, use prjfile instead",
+            PendingDeprecationWarning,
+        )
+        self.prjfile = prj
 
     @property
     def prjfile(self):
-        return self._prjfile
+        """
+        Path to a .prj file containing WKT for a coordinate reference system.
+        """
+        return getattr(self, "_prjfile", None)
 
     @prjfile.setter
     def prjfile(self, prjfile):
-        self._crs = get_crs(prjfile=prjfile)
+        if prjfile is None:
+            self._prjfile = None
+            return
+        if not isinstance(prjfile, (str, os.PathLike)):
+            raise ValueError("prjfile property must be str, PathLike or None")
         self._prjfile = prjfile
+        # If crs was previously unset, use .prj file input
+        if self._crs is None:
+            try:
+                self._crs = get_crs(prjfile=prjfile)
+            except (ImportError, FileNotFoundError):
+                pass
 
     @property
     def top(self):
@@ -910,11 +1022,66 @@ def set_coord_info(
         angrot=None,
         crs=None,
         prjfile=None,
-        epsg=None,
-        proj4=None,
         merge_coord_info=True,
+        **kwargs,
     ):
-        new_crs = get_crs(prjfile=prjfile, epsg=epsg, proj4=proj4, crs=crs)
+        """Set coordinate information for a grid.
+
+        Parameters
+        ----------
+        xoff, yoff : float, optional
+            X and Y coordinate of the origin point in the spatial reference
+            coordinate system.
+        angrot : float, optional
+            Rotation angle of model grid, as it is rotated around the origin
+            point.
+        crs : pyproj.CRS, int, str, optional
+            Coordinate reference system (CRS) for the model grid
+            (must be projected; geographic CRS are not supported).
+            The value can be anything accepted by
+            :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+            such as an authority string (eg "EPSG:26916") or a WKT string.
+        prjfile : str or pathlike, optional
+            ESRI-style projection file with well-known text defining the CRS
+            for the model grid (must be projected; geographic CRS are not
+            supported).
+        merge_coord_info : bool, default True
+            If True, retaining previous properties.
+            If False, overwrite properties with defaults, unless specified.
+        **kwargs : dict, optional
+            Support deprecated keyword options.
+
+            .. deprecated:: 3.3.7
+               The following keyword options will be removed for FloPy 3.4:
+
+                 - ``epsg`` (int): use ``crs`` instead.
+                 - ``proj4`` (str): use ``crs`` instead.
+
+        """
+        if crs is not None:
+            # Force these to be re-evaluated, if/when needed
+            if hasattr(self, "_epsg"):
+                delattr(self, "_epsg")
+            if hasattr(self, "_proj4"):
+                delattr(self, "_proj4")
+        # Handle deprecated projection kwargs; warnings are raised in crs.py
+        get_crs_args = {"crs": crs, "prjfile": prjfile}
+        if "epsg" in kwargs:
+            self.epsg = get_crs_args["epsg"] = kwargs.pop("epsg")
+        if "proj4" in kwargs:
+            self.proj4 = get_crs_args["proj4"] = kwargs.pop("proj4")
+        if kwargs:
+            raise TypeError(f"unhandled keywords: {kwargs}")
+        try:
+            new_crs = get_crs(**get_crs_args)
+        except ImportError:
+            new_crs = None
+            # provide some support without pyproj by retaining 'epsg' integer
+            if getattr(self, "_epsg", None) is None:
+                epsg = _get_epsg_from_crs_or_proj4(crs, self.proj4)
+                if epsg is not None:
+                    self.epsg = epsg
+
         if merge_coord_info:
             if xoff is None:
                 xoff = self._xoff
@@ -936,7 +1103,7 @@ def set_coord_info(
         self._yoff = yoff
         self._angrot = angrot
         self._prjfile = prjfile
-        self.crs = new_crs
+        self._crs = new_crs
         self._require_cache_updates()
 
     def load_coord_info(self, namefile=None, reffile="usgs.model.reference"):
@@ -959,6 +1126,7 @@ def attribs_from_namfile_header(self, namefile):
         if namefile is None:
             return False
         xul, yul = None, None
+        set_coord_info_args = {}
         header = []
         with open(namefile) as f:
             for line in f:
@@ -994,11 +1162,18 @@ def attribs_from_namfile_header(self, namefile):
                     self._angrot = float(item.split(":")[1])
                 except:
                     pass
+            elif "crs" in item.lower():
+                try:
+                    crs = ":".join(item.split(":")[1:]).strip()
+                    if crs.lower() != "none":
+                        set_coord_info_args["crs"] = crs
+                except:
+                    pass
             elif "proj4_str" in item.lower():
                 try:
-                    self._proj4 = ":".join(item.split(":")[1:]).strip()
-                    if self._proj4.lower() == "none":
-                        self._proj4 = None
+                    proj4 = ":".join(item.split(":")[1:]).strip()
+                    if proj4.lower() != "none":
+                        set_coord_info_args["proj4"] = crs
                 except:
                     pass
             elif "start" in item.lower():
@@ -1010,11 +1185,12 @@ def attribs_from_namfile_header(self, namefile):
         # we need to rotate the modelgrid first, then we can
         # calculate the xll and yll from xul and yul
         if (xul, yul) != (None, None):
-            self.set_coord_info(
-                xoff=self._xul_to_xll(xul),
-                yoff=self._yul_to_yll(yul),
-                angrot=self._angrot,
-            )
+            set_coord_info_args["xoff"] = self._xul_to_xll(xul)
+            set_coord_info_args["yoff"] = self._yul_to_yll(yul)
+            set_coord_info_args["angrot"] = self._angrot
+
+        if set_coord_info_args:
+            self.set_coord_info(**set_coord_info_args)
 
         return True
 
@@ -1042,7 +1218,7 @@ def read_usgs_model_reference_file(self, reffile="usgs.model.reference"):
                                 elif info[0] == "rotation":
                                     self._angrot = float(data)
                                 elif info[0] == "epsg":
-                                    self.crs = int(data)
+                                    self.epsg = int(data)
                                 elif info[0] == "proj4":
                                     self.crs = data
                                 elif info[0] == "start_date":
@@ -1111,21 +1287,32 @@ def _zcoords(self):
         return zbdryelevs, zcenters
 
     # Exporting
-    def write_shapefile(self, filename="grid.shp", epsg=None, prj=None):
+    def write_shapefile(
+        self, filename="grid.shp", crs=None, prjfile=None, **kwargs
+    ):
         """
         Write a shapefile of the grid with just the row and column attributes.
 
         """
         from ..export.shapefile_utils import write_grid_shapefile
 
+        # Handle deprecated projection kwargs; warnings are raised in crs.py
+        write_grid_shapefile_args = {}
+        if "epsg" in kwargs:
+            write_grid_shapefile_args["epsg"] = kwargs.pop("epsg")
+        if "prj" in kwargs:
+            write_grid_shapefile_args["prj"] = kwargs.pop("prj")
+        if kwargs:
+            raise TypeError(f"unhandled keywords: {kwargs}")
+        if crs is None:
+            crs = self.crs
         write_grid_shapefile(
             filename,
             self,
             array_dict={},
             nan_val=-1.0e9,
-            crs=self.crs,
-            epsg=epsg,
-            prj=prj,
+            crs=crs,
+            **write_grid_shapefile_args,
         )
         return
 
diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py
index d1ab5d732c..47576d61cd 100644
--- a/flopy/discretization/structuredgrid.py
+++ b/flopy/discretization/structuredgrid.py
@@ -97,7 +97,7 @@ class for a structured model grid
         ibound/idomain value for each cell
     lenuni : int or ndarray
         model length units
-    crs : pyproj.CRS, optional if `prjfile` is specified
+    crs : pyproj.CRS, int, str, optional if `prjfile` is specified
         Coordinate reference system (CRS) for the model grid
         (must be projected; geographic CRS are not supported).
         The value can be anything accepted by
@@ -114,6 +114,15 @@ class for a structured model grid
         in the spatial reference coordinate system
     angrot : float
         rotation angle of model grid, as it is rotated around the origin point
+    **kwargs : dict, optional
+        Support deprecated keyword options.
+
+        .. deprecated:: 3.3.7
+           The following keyword options will be removed for FloPy 3.4:
+
+             - ``prj`` (str or pathlike): use ``prjfile`` instead.
+             - ``epsg`` (int): use ``crs`` instead.
+             - ``proj4`` (str): use ``crs`` instead.
 
     Properties
     ----------
@@ -146,9 +155,6 @@ def __init__(
         idomain=None,
         lenuni=None,
         crs=None,
-        epsg=None,
-        proj4=None,
-        prj=None,
         prjfile=None,
         xoff=0.0,
         yoff=0.0,
@@ -157,21 +163,20 @@ def __init__(
         nrow=None,
         ncol=None,
         laycbd=None,
+        **kwargs,
     ):
         super().__init__(
             "structured",
-            top,
-            botm,
-            idomain,
-            lenuni,
-            crs,
-            epsg,
-            proj4,
-            prj,
-            prjfile,
-            xoff,
-            yoff,
-            angrot,
+            top=top,
+            botm=botm,
+            idomain=idomain,
+            lenuni=lenuni,
+            crs=crs,
+            prjfile=prjfile,
+            xoff=xoff,
+            yoff=yoff,
+            angrot=angrot,
+            **kwargs,
         )
         if delc is not None:
             self.__nrow = len(delc)
diff --git a/flopy/discretization/unstructuredgrid.py b/flopy/discretization/unstructuredgrid.py
index aca36444b6..5f799e78c6 100644
--- a/flopy/discretization/unstructuredgrid.py
+++ b/flopy/discretization/unstructuredgrid.py
@@ -51,7 +51,7 @@ class UnstructuredGrid(Grid):
         If the model grid defined in verts and iverts applies for all model
         layers, then the length of iverts can be equal to ncpl[0] and there
         is no need to repeat all of the vertex information for cells in layers
-    crs : pyproj.CRS, optional if `prjfile` is specified
+    crs : pyproj.CRS, int, str, optional if `prjfile` is specified
         Coordinate reference system (CRS) for the model grid
         (must be projected; geographic CRS are not supported).
         The value can be anything accepted by
@@ -73,6 +73,15 @@ class UnstructuredGrid(Grid):
         optional number of connections per node array
     ja : list or ndarray
         optional jagged connection array
+    **kwargs : dict, optional
+        Support deprecated keyword options.
+
+        .. deprecated:: 3.3.7
+           The following keyword options will be removed for FloPy 3.4:
+
+             - ``prj`` (str or pathlike): use ``prjfile`` instead.
+             - ``epsg`` (int): use ``crs`` instead.
+             - ``proj4`` (str): use ``crs`` instead.
 
     Properties
     ----------
@@ -117,30 +126,26 @@ def __init__(
         lenuni=None,
         ncpl=None,
         crs=None,
-        epsg=None,
-        proj4=None,
-        prj=None,
         prjfile=None,
         xoff=0.0,
         yoff=0.0,
         angrot=0.0,
         iac=None,
         ja=None,
+        **kwargs,
     ):
         super().__init__(
             "unstructured",
-            top,
-            botm,
-            idomain,
-            lenuni,
-            crs,
-            epsg,
-            proj4,
-            prj,
-            prjfile,
-            xoff,
-            yoff,
-            angrot,
+            top=top,
+            botm=botm,
+            idomain=idomain,
+            lenuni=lenuni,
+            crs=crs,
+            prjfile=prjfile,
+            xoff=xoff,
+            yoff=yoff,
+            angrot=angrot,
+            **kwargs,
         )
 
         # if any of these are None, then the grid is not valid
diff --git a/flopy/discretization/vertexgrid.py b/flopy/discretization/vertexgrid.py
index 340439d3b2..e35da8a366 100644
--- a/flopy/discretization/vertexgrid.py
+++ b/flopy/discretization/vertexgrid.py
@@ -27,7 +27,7 @@ class for a vertex model grid
         ibound/idomain value for each cell
     lenuni : int or ndarray
         model length units
-    crs : pyproj.CRS, optional if `prjfile` is specified
+    crs : pyproj.CRS, int, str, optional if `prjfile` is specified
         Coordinate reference system (CRS) for the model grid
         (must be projected; geographic CRS are not supported).
         The value can be anything accepted by
@@ -44,6 +44,15 @@ class for a vertex model grid
         in the spatial reference coordinate system
     angrot : float
         rotation angle of model grid, as it is rotated around the origin point
+    **kwargs : dict, optional
+        Support deprecated keyword options.
+
+        .. deprecated:: 3.3.7
+           The following keyword options will be removed for FloPy 3.4:
+
+             - ``prj`` (str or pathlike): use ``prjfile`` instead.
+             - ``epsg`` (int): use ``crs`` instead.
+             - ``proj4`` (str): use ``crs`` instead.
 
     Properties
     ----------
@@ -68,9 +77,6 @@ def __init__(
         idomain=None,
         lenuni=None,
         crs=None,
-        epsg=None,
-        proj4=None,
-        prj=None,
         prjfile=None,
         xoff=0.0,
         yoff=0.0,
@@ -78,21 +84,20 @@ def __init__(
         nlay=None,
         ncpl=None,
         cell1d=None,
+        **kwargs,
     ):
         super().__init__(
             "vertex",
-            top,
-            botm,
-            idomain,
-            lenuni,
-            crs,
-            epsg,
-            proj4,
-            prj,
-            prjfile,
-            xoff,
-            yoff,
-            angrot,
+            top=top,
+            botm=botm,
+            idomain=idomain,
+            lenuni=lenuni,
+            crs=crs,
+            prjfile=prjfile,
+            xoff=xoff,
+            yoff=yoff,
+            angrot=angrot,
+            **kwargs,
         )
         self._vertices = vertices
         self._cell1d = cell1d
diff --git a/flopy/export/shapefile_utils.py b/flopy/export/shapefile_utils.py
index 3e31211f3e..e73a550dc2 100644
--- a/flopy/export/shapefile_utils.py
+++ b/flopy/export/shapefile_utils.py
@@ -55,7 +55,10 @@ def write_gridlines_shapefile(filename: Union[str, os.PathLike], mg):
         wr.record(i)
 
     wr.close()
-    write_prj(filename, modelgrid=mg)
+    try:
+        write_prj(filename, modelgrid=mg)
+    except ImportError:
+        pass
     return
 
 
@@ -65,10 +68,9 @@ def write_grid_shapefile(
     array_dict,
     nan_val=np.nan,
     crs=None,
-    prjfile=None,
-    epsg=None,
-    prj: Optional[Union[str, os.PathLike]] = None,
+    prjfile: Optional[Union[str, os.PathLike]] = None,
     verbose=False,
+    **kwargs,
 ):
     """
     Method to write a shapefile of gridded input data
@@ -83,7 +85,7 @@ def write_grid_shapefile(
         dictionary of model input arrays
     nan_val : float
         value to fill nans
-    crs : pyproj.CRS, optional if `prjfile` is specified
+    crs : pyproj.CRS, int, str, optional if `prjfile` is specified
         Coordinate reference system (CRS) for the model grid
         (must be projected; geographic CRS are not supported).
         The value can be anything accepted by
@@ -92,6 +94,14 @@ def write_grid_shapefile(
     prjfile : str or pathlike, optional if `crs` is specified
         ESRI-style projection file with well-known text defining the CRS
         for the model grid (must be projected; geographic CRS are not supported).
+    **kwargs : dict, optional
+        Support deprecated keyword options.
+
+        .. deprecated:: 3.3.7
+           The following keyword options will be removed for FloPy 3.4:
+
+             - ``epsg`` (int): use ``crs`` instead.
+             - ``prj`` (str or PathLike): use ``prjfile`` instead.
 
     Returns
     -------
@@ -213,8 +223,20 @@ def write_grid_shapefile(
     if verbose:
         print(f"wrote {flopy_io.relpath_safe(path)}")
 
+    # handle deprecated projection kwargs; warnings are raised in crs.py
+    write_prj_args = {}
+    if "epsg" in kwargs:
+        write_prj_args["epsg"] = kwargs.pop("epsg")
+    if "prj" in kwargs:
+        write_prj_args["prj"] = kwargs.pop("prj")
+    if kwargs:
+        raise TypeError(f"unhandled keywords: {kwargs}")
     # write the projection file
-    write_prj(path, mg, crs=crs, epsg=epsg, prj=prj, prjfile=prjfile)
+    try:
+        write_prj(path, mg, crs=crs, prjfile=prjfile, **write_prj_args)
+    except ImportError:
+        if verbose:
+            print("projection file not written")
     return
 
 
@@ -248,7 +270,7 @@ def model_attributes_to_shapefile(
         modelgrid : fp.modflow.Grid object
             if modelgrid is supplied, user supplied modelgrid is used in lieu
             of the modelgrid attached to the modflow model object
-        crs : pyproj.CRS, optional if `prjfile` is specified
+        crs : pyproj.CRS, int, str, optional if `prjfile` is specified
             Coordinate reference system (CRS) for the model grid
             (must be projected; geographic CRS are not supported).
             The value can be anything accepted by
@@ -403,7 +425,10 @@ def model_attributes_to_shapefile(
     write_grid_shapefile(path, grid, array_dict)
     crs = kwargs.get("crs", None)
     prjfile = kwargs.get("prjfile", None)
-    write_prj(path, grid, crs=crs, prjfile=prjfile)
+    try:
+        write_prj(path, grid, crs=crs, prjfile=prjfile)
+    except ImportError:
+        pass
 
 
 def shape_attr_name(name, length=6, keep_layer=False):
@@ -547,9 +572,7 @@ def recarray2shp(
     shpname: Union[str, os.PathLike] = "recarray.shp",
     mg=None,
     crs=None,
-    prjfile=None,
-    epsg=None,
-    prj: Optional[Union[str, os.PathLike]] = None,
+    prjfile: Optional[Union[str, os.PathLike]] = None,
     verbose=False,
     **kwargs,
 ):
@@ -571,7 +594,9 @@ def recarray2shp(
         recarray.
     shpname : str or PathLike, default "recarray.shp"
         Path for the output shapefile
-    crs : pyproj.CRS, optional if `prjfile` is specified
+    mg : flopy.discretization.Grid object
+        flopy model grid
+    crs : pyproj.CRS, int, str, optional if `prjfile` is specified
         Coordinate reference system (CRS) for the model grid
         (must be projected; geographic CRS are not supported).
         The value can be anything accepted by
@@ -580,10 +605,18 @@ def recarray2shp(
     prjfile : str or pathlike, optional if `crs` is specified
         ESRI-style projection file with well-known text defining the CRS
         for the model grid (must be projected; geographic CRS are not supported).
+    **kwargs : dict, optional
+        Support deprecated keyword options.
+
+        .. deprecated:: 3.3.7
+           The following keyword options will be removed for FloPy 3.4:
+
+             - ``epsg`` (int): use ``crs`` instead.
+             - ``prj`` (str or PathLike): use ``prjfile`` instead.
 
     Notes
     -----
-    Uses pyshp.
+    Uses pyshp and optionally pyproj.
     """
     from ..utils.geospatial_utils import GeoSpatialCollection
 
@@ -637,8 +670,23 @@ def recarray2shp(
             w.record(*r)
 
     w.close()
-    write_prj(shpname, mg, crs=crs, epsg=epsg, prj=prj, prjfile=prjfile)
-    print(f"wrote {flopy_io.relpath_safe(os.getcwd(), shpname)}")
+    if verbose:
+        print(f"wrote {flopy_io.relpath_safe(os.getcwd(), shpname)}")
+
+    # handle deprecated projection kwargs; warnings are raised in crs.py
+    write_prj_args = {}
+    if "epsg" in kwargs:
+        write_prj_args["epsg"] = kwargs.pop("epsg")
+    if "prj" in kwargs:
+        write_prj_args["prj"] = kwargs.pop("prj")
+    if kwargs:
+        raise TypeError(f"unhandled keywords: {kwargs}")
+    # write the projection file
+    try:
+        write_prj(shpname, mg, crs=crs, prjfile=prjfile, **write_prj_args)
+    except ImportError:
+        if verbose:
+            print("projection file not written")
     return
 
 
@@ -646,26 +694,30 @@ def write_prj(
     shpname,
     modelgrid=None,
     crs=None,
-    epsg=None,
-    prj=None,
     prjfile=None,
-    wkt_string=None,
+    **kwargs,
 ):
     # projection file name
     output_projection_file = Path(shpname).with_suffix(".prj")
 
-    crs = get_crs(
-        prjfile=prjfile, prj=prj, epsg=epsg, crs=crs, wkt_string=wkt_string
-    )
+    # handle deprecated projection kwargs; warnings are raised in crs.py
+    get_crs_args = {}
+    if "epsg" in kwargs:
+        get_crs_args["epsg"] = kwargs.pop("epsg")
+    if "prj" in kwargs:
+        get_crs_args["prj"] = kwargs.pop("prj")
+    if "wkt_string" in kwargs:
+        get_crs_args["wkt_string"] = kwargs.pop("wkt_string")
+    if kwargs:
+        raise TypeError(f"unhandled keywords: {kwargs}")
+    crs = get_crs(prjfile=prjfile, crs=crs, **get_crs_args)
     if crs is None and modelgrid is not None:
         crs = modelgrid.crs
     if crs is not None:
-        with open(output_projection_file, "w", encoding="utf-8") as dest:
-            write_text = crs.to_wkt()
-            dest.write(write_text)
+        output_projection_file.write_text(crs.to_wkt(), encoding="utf-8")
     else:
         print(
             "No CRS information for writing a .prj file.\n"
-            "Supply an valid coordinate system reference to the attached modelgrid object "
-            "or .export() method."
+            "Supply an valid coordinate system reference to the attached "
+            "modelgrid object or .export() method."
         )
diff --git a/flopy/export/utils.py b/flopy/export/utils.py
index 3b94cb10a9..81d02f0771 100644
--- a/flopy/export/utils.py
+++ b/flopy/export/utils.py
@@ -596,7 +596,7 @@ def model_export(
         modelgrid: flopy.discretization.Grid
             user supplied modelgrid object which will supercede the built
             in modelgrid object
-        crs : pyproj.CRS, optional if `prjfile` is specified
+        crs : pyproj.CRS, int, str, optional if `prjfile` is specified
             Coordinate reference system (CRS) for the model grid
             (must be projected; geographic CRS are not supported).
             The value can be anything accepted by
@@ -689,7 +689,7 @@ def package_export(
         modelgrid: flopy.discretization.Grid
             user supplied modelgrid object which will supercede the built
             in modelgrid object
-        crs : pyproj.CRS, optional if `prjfile` is specified
+        crs : pyproj.CRS, int, str, optional if `prjfile` is specified
             Coordinate reference system (CRS) for the model grid
             (must be projected; geographic CRS are not supported).
             The value can be anything accepted by
@@ -888,7 +888,7 @@ def mflist_export(f: Union[str, os.PathLike, NetCdf], mfl, **kwargs):
     **kwargs : keyword arguments
         modelgrid : flopy.discretization.Grid
             model grid instance which will supercede the flopy.model.modelgrid
-        crs : pyproj.CRS, optional if `prjfile` is specified
+        crs : pyproj.CRS, int, str, optional if `prjfile` is specified
             Coordinate reference system (CRS) for the model grid
             (must be projected; geographic CRS are not supported).
             The value can be anything accepted by
@@ -1679,7 +1679,10 @@ def export_array(
     elif filename.lower().endswith(".shp"):
         from ..export.shapefile_utils import write_grid_shapefile
 
-        crs = get_crs(**kwargs)
+        try:
+            crs = get_crs(**kwargs)
+        except ImportError:
+            crs = None
         write_grid_shapefile(
             filename,
             modelgrid,
@@ -1866,7 +1869,7 @@ def export_array_contours(
     ax = plt.subplots()[-1]
     ctr = contour_array(modelgrid, ax, a, levels=levels)
 
-    kwargs["modelgrid"] = modelgrid
+    kwargs["mg"] = modelgrid
     export_contours(filename, ctr, fieldname, **kwargs)
     plt.close()
 
diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py
index a47cf106c4..94f2e9438a 100644
--- a/flopy/mf6/mfmodel.py
+++ b/flopy/mf6/mfmodel.py
@@ -622,10 +622,6 @@ def export(self, f, **kwargs):
             modelgrid: flopy.discretization.Grid
                 User supplied modelgrid object which will supercede the built
                 in modelgrid object
-            epsg : int
-                EPSG projection code
-            prj : str
-                The prj file name
             if fmt is set to 'vtk', parameters of vtk.export_model
 
         """
diff --git a/flopy/modflow/mf.py b/flopy/modflow/mf.py
index 9962f48583..ea1ec9974e 100644
--- a/flopy/modflow/mf.py
+++ b/flopy/modflow/mf.py
@@ -279,6 +279,12 @@ def modelgrid(self):
         else:
             ibound = None
 
+        common_kwargs = {
+            "crs": self._modelgrid.crs or self._modelgrid.epsg,
+            "xoff": self._modelgrid.xoffset,
+            "yoff": self._modelgrid.yoffset,
+            "angrot": self._modelgrid.angrot,
+        }
         if self.get_package("disu") is not None:
             # build unstructured grid
             self._modelgrid = UnstructuredGrid(
@@ -292,12 +298,9 @@ def modelgrid(self):
                 botm=self.disu.bot.array,
                 idomain=ibound,
                 lenuni=self.disu.lenuni,
-                crs=self._modelgrid.crs,
-                xoff=self._modelgrid.xoffset,
-                yoff=self._modelgrid.yoffset,
-                angrot=self._modelgrid.angrot,
                 iac=self.disu.iac.array,
                 ja=self.disu.ja.array,
+                **common_kwargs,
             )
             print(
                 "WARNING: Model grid functionality limited for unstructured "
@@ -312,12 +315,9 @@ def modelgrid(self):
                 self.dis.botm.array,
                 ibound,
                 self.dis.lenuni,
-                crs=self._modelgrid.crs,
-                xoff=self._modelgrid.xoffset,
-                yoff=self._modelgrid.yoffset,
-                angrot=self._modelgrid.angrot,
                 nlay=self.dis.nlay,
                 laycbd=self.dis.laycbd.array,
+                **common_kwargs,
             )
 
         # resolve offsets
@@ -337,7 +337,7 @@ def modelgrid(self):
             xoff,
             yoff,
             self._modelgrid.angrot,
-            self._modelgrid.crs,
+            self._modelgrid.crs or self._modelgrid.epsg,
         )
         self._mg_resync = not self._modelgrid.is_complete
         return self._modelgrid
diff --git a/flopy/modflow/mfdis.py b/flopy/modflow/mfdis.py
index 85625f3af8..9ff5838867 100644
--- a/flopy/modflow/mfdis.py
+++ b/flopy/modflow/mfdis.py
@@ -88,7 +88,7 @@ class ModflowDis(Package):
     rotation : float
         counter-clockwise rotation (in degrees) of the grid about the lower-
         left corner. default is 0.0
-    crs : pyproj.CRS, optional if `prjfile` is specified
+    crs : pyproj.CRS, int, str, optional if `prjfile` is specified
         Coordinate reference system (CRS) for the model grid
         (must be projected; geographic CRS are not supported).
         The value can be anything accepted by
@@ -99,6 +99,11 @@ class ModflowDis(Package):
         for the model grid (must be projected; geographic CRS are not supported).
     start_datetime : str
         starting datetime of the simulation. default is '1/1/1970'
+    **kwargs : dict, optional
+        Support deprecated keyword options.
+
+        .. deprecated:: 3.3.7
+           ``proj4_str`` will be removed for FloPy 3.4, use ``crs`` instead.
 
     Attributes
     ----------
@@ -147,10 +152,10 @@ def __init__(
         xul=None,
         yul=None,
         rotation=None,
-        proj4_str=None,
         crs=None,
         prjfile=None,
         start_datetime=None,
+        **kwargs,
     ):
         # set default unit number of one is not specified
         if unitnumber is None:
@@ -242,30 +247,37 @@ def __init__(
             5: "years",
         }
 
-        if xul is None:
-            xul = model._xul
-        if yul is None:
-            yul = model._yul
-        if rotation is None:
-            rotation = model._rotation
-        crs = get_crs(prjfile=prjfile, proj4=proj4_str, crs=crs)
-        if crs is None:
-            crs = model._crs
-        if start_datetime is None:
-            start_datetime = model._start_datetime
-
         # set the model grid coordinate info
-        xll = None
-        yll = None
         mg = model.modelgrid
+        if rotation is None:
+            rotation = model._rotation
         if rotation is not None:
-            mg.set_coord_info(xoff=None, yoff=None, angrot=rotation)
+            # set rotation before anything else
+            mg.set_coord_info(angrot=rotation)
+        set_coord_info_args = {"crs": crs, "prjfile": prjfile}
+        if "proj4_str" in kwargs:
+            warnings.warn(
+                "the proj4_str argument will be deprecated and will be "
+                "removed in version 3.4. Use crs instead.",
+                PendingDeprecationWarning,
+            )
+            proj4_str = kwargs.pop("proj4_str")
+            if crs is None:
+                set_coord_info_args["crs"] = proj4_str
+        if xul is None:
+            xul = model._xul
         if xul is not None:
-            xll = mg._xul_to_xll(xul)
+            set_coord_info_args["xoff"] = mg._xul_to_xll(xul)
+        if yul is None:
+            yul = model._yul
         if yul is not None:
-            yll = mg._yul_to_yll(yul)
-        mg.set_coord_info(xoff=xll, yoff=yll, angrot=rotation, crs=crs)
+            set_coord_info_args["yoff"] = mg._yul_to_yll(yul)
+
+        # set all other relevant coordinate properties
+        model._modelgrid.set_coord_info(**set_coord_info_args)
 
+        if start_datetime is None:
+            start_datetime = model._start_datetime
         self.tr = TemporalReference(
             itmuni=self.itmuni, start_datetime=start_datetime
         )
diff --git a/flopy/utils/crs.py b/flopy/utils/crs.py
index a41dcaa141..698e204fa6 100644
--- a/flopy/utils/crs.py
+++ b/flopy/utils/crs.py
@@ -13,7 +13,7 @@ def get_authority_crs(crs):
 
     Parameters
     ----------
-    crs : pyproj.CRS
+    crs : pyproj.CRS, int, str
         Coordinate reference system (CRS) for the model grid
         (must be projected; geographic CRS are not supported).
         The value can be anything accepted by
@@ -22,9 +22,9 @@ def get_authority_crs(crs):
 
     Returns
     -------
-    authority_crs : pyproj.CRS instance
+    pyproj.CRS instance
         CRS instance initiallized with the name
-        and authority code (e.g. epsg: 5070) produced by
+        and authority code (e.g. "epsg:5070") produced by
         :meth:`pyproj.crs.CRS.to_authority`
 
     Notes
@@ -53,27 +53,30 @@ def get_shapefile_crs(shapefile):
     Parameters
     ----------
     shapefile : str or pathlike
-        Path to a shapefile or an associated
-        projection (.prj) file.
+        Path to a shapefile or an associated projection (.prj) file.
 
     Returns
     -------
-    crs : pyproj.CRS instance
+    pyproj.CRS instance
 
     """
-    pyproj = import_optional_dependency("pyproj")
     shapefile = Path(shapefile)
     prjfile = shapefile.with_suffix(".prj")
     if prjfile.exists():
-        with open(prjfile, encoding="utf-8") as src:
-            wkt = src.read()
-            crs = pyproj.crs.CRS.from_wkt(wkt)
-            return get_authority_crs(crs)
+        pyproj = import_optional_dependency("pyproj")
+        wkt = prjfile.read_text(encoding="utf-8")
+        crs = pyproj.crs.CRS.from_wkt(wkt)
+        return get_authority_crs(crs)
+    else:
+        raise FileNotFoundError(f"could not find {prjfile}")
 
 
 def get_crs(prjfile=None, crs=None, **kwargs):
-    """Helper function to produce a pyproj.CRS object from
-    various input. Longer-term, this would just handle the ``crs``
+    """Helper function to produce a pyproj.CRS object from various input.
+
+    Notes
+    -----
+    Longer-term, this would just handle the ``crs``
     and ``prjfile`` arguments, but in the near term, we need to
     warn users about deprecating the ``prj``, ``epsg``, ``proj4``
     and ``wkt_string`` inputs.
@@ -81,49 +84,49 @@ def get_crs(prjfile=None, crs=None, **kwargs):
     Parameters
     ----------
     prjfile : str or pathlike, optional
-        _description_, by default None
-    prj : str or pathlike, optional
-        .. deprecated:: 3.4
-        use ``prjfile`` instead.
-    epsg : int, optional
-        .. deprecated:: 3.4
-        use ``crs`` instead.
-    proj4 : str, optional
-        .. deprecated:: 3.4
-        use ``crs`` instead.
-    crs : pyproj.CRS, optional if `prjfile` is specified
+        ESRI-style projection file with well-known text defining the CRS
+        for the model grid (must be projected; geographic CRS are not supported).
+    crs : pyproj.CRS, int, str, optional if `prjfile` is specified
         Coordinate reference system (CRS) for the model grid
         (must be projected; geographic CRS are not supported).
         The value can be anything accepted by
         :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
         such as an authority string (eg "EPSG:26916") or a WKT string.
-    wkt_string : str, optional
-        .. deprecated:: 3.4
-        use ``crs`` instead.
+    **kwargs : dict, optional
+        Support deprecated keyword options.
+
+        .. deprecated:: 3.3.7
+           The following keyword options will be removed for FloPy 3.4:
+
+             - ``prj`` (str or pathlike): use ``prjfile`` instead.
+             - ``epsg`` (int): use ``crs`` instead.
+             - ``proj4`` (str): use ``crs`` instead.
+             - ``wkt_string`` (str): use ``crs`` instead.
 
     Returns
     -------
-    crs : pyproj.CRS instance
+    pyproj.CRS instance
 
     """
     if crs is not None:
         crs = get_authority_crs(crs)
-    if kwargs.get("prj") is not None:
+    if "prj" in kwargs:
         warnings.warn(
-            "the prj argument will be deprecated and will be removed in version "
-            "3.4. Use prjfile instead.",
+            "the prj argument will be deprecated and will be removed in "
+            "version 3.4. Use prjfile instead.",
             PendingDeprecationWarning,
         )
-        prjfile = kwargs.get("prj")
-    if kwargs.get("epsg") is not None:
+        prjfile = kwargs.pop("prj")
+    if "epsg" in kwargs:
         warnings.warn(
-            "the epsg argument will be deprecated and will be removed in version "
-            "3.4. Use crs instead.",
+            "the epsg argument will be deprecated and will be removed in "
+            "version 3.4. Use crs instead.",
             PendingDeprecationWarning,
         )
+        epsg = kwargs.pop("epsg")
         if crs is None:
-            crs = get_authority_crs(kwargs.get("epsg"))
-    elif prjfile is not None:
+            crs = get_authority_crs(epsg)
+    if prjfile is not None:
         prjfile_crs = get_shapefile_crs(prjfile)
         if (crs is not None) and (crs != prjfile_crs):
             raise ValueError(
@@ -133,17 +136,21 @@ def get_crs(prjfile=None, crs=None, **kwargs):
             )
         else:
             crs = prjfile_crs
-    elif kwargs.get("proj4") is not None:
+    if "proj4" in kwargs:
         warnings.warn(
-            "the proj4 argument will be deprecated and will be removed in version "
-            "3.4. Use crs instead.",
+            "the proj4 argument will be deprecated and will be removed in "
+            "version 3.4. Use crs instead.",
             PendingDeprecationWarning,
         )
+        proj4 = kwargs.pop("proj4")
         if crs is None:
-            crs = get_authority_crs(kwargs.get("proj4"))
-    elif kwargs.get("wkt_string") is not None:
+            crs = get_authority_crs(proj4)
+    if "wkt_string" in kwargs:
+        wkt_string = kwargs.pop("wkt_string")
         if crs is None:
-            crs = get_authority_crs(kwargs.get("wkt_string"))
+            crs = get_authority_crs(wkt_string)
+    if kwargs:
+        raise TypeError(f"unhandled keywords: {kwargs}")
     if crs is not None and not crs.is_projected:
         raise ValueError(
             f"Only projected coordinate reference systems are supported.\n{crs}"
diff --git a/flopy/utils/mfreadnam.py b/flopy/utils/mfreadnam.py
index be900cde26..d51854fffe 100644
--- a/flopy/utils/mfreadnam.py
+++ b/flopy/utils/mfreadnam.py
@@ -213,7 +213,6 @@ def attribs_from_namfile_header(namefile):
         "yul": None,
         "rotation": 0.0,
         "crs": None,
-        "proj4_str": None,
     }
     if namefile is None:
         return defaults
@@ -256,21 +255,22 @@ def attribs_from_namfile_header(namefile):
             except:
                 print(f"   could not parse rotation in {namefile}")
         elif "proj4_str" in item.lower():
+            # deprecated, use "crs" instead
             try:
                 proj4 = ":".join(item.split(":")[1:]).strip()
                 if proj4.lower() == "none":
                     proj4 = None
-                defaults["crs"] = proj4
+                defaults["proj4_str"] = proj4
             except:
                 print(f"   could not parse proj4_str in {namefile}")
         elif "crs" in item.lower():
             try:
                 crs = ":".join(item.split(":")[1:]).strip()
                 if crs.lower() == "none":
-                    proj4 = None
+                    crs = None
                 defaults["crs"] = crs
             except:
-                print(f"   could not parse proj4_str in {namefile}")
+                print(f"   could not parse crs in {namefile}")
         elif "start" in item.lower():
             try:
                 start_datetime = item.split(":")[1].strip()
diff --git a/flopy/utils/modpathfile.py b/flopy/utils/modpathfile.py
index 3a1a70193a..6af0023f5e 100644
--- a/flopy/utils/modpathfile.py
+++ b/flopy/utils/modpathfile.py
@@ -272,7 +272,7 @@ def write_shapefile(
         direction="ending",
         shpname="endpoints.shp",
         mg=None,
-        epsg=None,
+        crs=None,
         **kwargs,
     ):
         """
@@ -295,12 +295,19 @@ def write_shapefile(
             File path for shapefile
         mg : flopy.discretization.grid instance
             Used to scale and rotate Global x,y,z values.
-        epsg : int
-            EPSG code for writing projection (.prj) file. If this is not
-            supplied, the proj4 string or epgs code associated with mg will be
-            used.
+        crs : pyproj.CRS, int, str, optional
+            Coordinate reference system (CRS) for the model grid
+            (must be projected; geographic CRS are not supported).
+            The value can be anything accepted by
+            :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+            such as an authority string (eg "EPSG:26916") or a WKT string.
         kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp
 
+          .. deprecated:: 3.3.7
+             The following keyword options will be removed for FloPy 3.4:
+
+               - ``epsg`` (int): use ``crs`` instead.
+
         """
         from ..discretization import StructuredGrid
         from ..export.shapefile_utils import recarray2shp
@@ -325,9 +332,6 @@ def write_shapefile(
         if mg is None:
             raise ValueError("A modelgrid object was not provided.")
 
-        if epsg is None:
-            epsg = mg.epsg
-
         particles = np.unique(series.particleid)
         geoms = []
 
@@ -406,7 +410,7 @@ def write_shapefile(
             sdata[n] += 1
 
         # write the final recarray to a shapefile
-        recarray2shp(sdata, geoms, shpname=shpname, epsg=epsg, **kwargs)
+        recarray2shp(sdata, geoms, shpname=shpname, crs=crs, **kwargs)
 
 
 class PathlineFile(_ModpathSeries):
@@ -744,7 +748,7 @@ def write_shapefile(
         direction="ending",
         shpname="pathlines.shp",
         mg=None,
-        epsg=None,
+        crs=None,
         **kwargs,
     ):
         """
@@ -769,12 +773,19 @@ def write_shapefile(
         mg : flopy.discretization.grid instance
             Used to scale and rotate Global x,y,z values in MODPATH Pathline
             file.
-        epsg : int
-            EPSG code for writing projection (.prj) file. If this is not
-            supplied, the proj4 string or epgs code associated with mg will be
-            used.
+        crs : pyproj.CRS, int, str, optional
+            Coordinate reference system (CRS) for the model grid
+            (must be projected; geographic CRS are not supported).
+            The value can be anything accepted by
+            :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+            such as an authority string (eg "EPSG:26916") or a WKT string.
         kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp
 
+          .. deprecated:: 3.3.7
+             The following keyword options will be removed for FloPy 3.4:
+
+               - ``epsg`` (int): use ``crs`` instead.
+
         """
         super().write_shapefile(
             data=pathline_data,
@@ -782,7 +793,7 @@ def write_shapefile(
             direction=direction,
             shpname=shpname,
             mg=mg,
-            epsg=epsg,
+            crs=crs,
             **kwargs,
         )
 
@@ -1191,7 +1202,7 @@ def write_shapefile(
         shpname="endpoints.shp",
         direction="ending",
         mg=None,
-        epsg=None,
+        crs=None,
         **kwargs,
     ):
         """
@@ -1208,12 +1219,19 @@ def write_shapefile(
         mg : flopy.discretization.grid instance
             Used to scale and rotate Global x,y,z values in MODPATH Endpoint
             file.
-        epsg : int
-            EPSG code for writing projection (.prj) file. If this is not
-            supplied, the proj4 string or epgs code associated with mg will be
-            used.
+        crs : pyproj.CRS, int, str, optional
+            Coordinate reference system (CRS) for the model grid
+            (must be projected; geographic CRS are not supported).
+            The value can be anything accepted by
+            :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+            such as an authority string (eg "EPSG:26916") or a WKT string.
         kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp
 
+          .. deprecated:: 3.3.7
+             The following keyword options will be removed for FloPy 3.4:
+
+               - ``epsg`` (int): use ``crs`` instead.
+
         """
         from ..discretization import StructuredGrid
         from ..export.shapefile_utils import recarray2shp
@@ -1235,8 +1253,6 @@ def write_shapefile(
             )
         if mg is None:
             raise ValueError("A modelgrid object was not provided.")
-        if epsg is None:
-            epsg = mg.epsg
 
         if isinstance(mg, StructuredGrid):
             x, y = geometry.transform(
@@ -1255,7 +1271,7 @@ def write_shapefile(
         for n in self.kijnames:
             if n in epd.dtype.names:
                 epd[n] += 1
-        recarray2shp(epd, geoms, shpname=shpname, epsg=epsg, **kwargs)
+        recarray2shp(epd, geoms, shpname=shpname, crs=crs, **kwargs)
 
 
 class TimeseriesFile(_ModpathSeries):
@@ -1583,7 +1599,7 @@ def write_shapefile(
         direction="ending",
         shpname="pathlines.shp",
         mg=None,
-        epsg=None,
+        crs=None,
         **kwargs,
     ):
         """
@@ -1608,12 +1624,19 @@ def write_shapefile(
         mg : flopy.discretization.grid instance
             Used to scale and rotate Global x,y,z values in MODPATH Timeseries
             file.
-        epsg : int
-            EPSG code for writing projection (.prj) file. If this is not
-            supplied, the proj4 string or epgs code associated with mg will be
-            used.
+        crs : pyproj.CRS, int, str, optional
+            Coordinate reference system (CRS) for the model grid
+            (must be projected; geographic CRS are not supported).
+            The value can be anything accepted by
+            :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`,
+            such as an authority string (eg "EPSG:26916") or a WKT string.
         kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp
 
+          .. deprecated:: 3.3.7
+             The following keyword options will be removed for FloPy 3.4:
+
+               - ``epsg`` (int): use ``crs`` instead.
+
         """
         super().write_shapefile(
             data=timeseries_data,
@@ -1621,6 +1644,6 @@ def write_shapefile(
             direction=direction,
             shpname=shpname,
             mg=mg,
-            epsg=epsg,
+            crs=crs,
             **kwargs,
         )