diff --git a/CHANGES.rst b/CHANGES.rst
index d58150c5..c506eec6 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -1,6 +1,7 @@
 0.22.0 (unreleased)
 -------------------
 
+- replace usages of ``copy_arrays`` with ``memmap`` [#503]
 
 0.21.0 (2024-03-10)
 -------------------
diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py
index 1ff78bae..f07fe9f9 100644
--- a/gwcs/tests/test_wcs.py
+++ b/gwcs/tests/test_wcs.py
@@ -14,6 +14,7 @@
 from astropy import wcs as astwcs
 from astropy.wcs import wcsapi
 from astropy.time import Time
+from astropy.utils.introspection import minversion
 import asdf
 
 from .. import wcs
@@ -50,6 +51,14 @@
 y = np.linspace(0, 1, ny)
 xv, yv = np.meshgrid(x, y)
 
+
+def asdf_open_memory_mapping_kwarg(memmap: bool) -> dict:
+    if minversion("asdf", "3.1.0"):
+        return {"memmap": memmap}
+    else :
+        return {"copy_arrays": not memmap}
+
+
 # Test initializing a WCS
 
 
@@ -648,7 +657,7 @@ def test_to_fits_sip():
     xflat = np.ravel(x[1:-1, 1:-1])
     yflat = np.ravel(y[1:-1, 1:-1])
     fn = os.path.join(data_path, 'miriwcs.asdf')
-    with asdf.open(fn, lazy_load=False, copy_arrays=True, ignore_missing_extensions=True) as af:
+    with asdf.open(fn, lazy_load=False, ignore_missing_extensions=True, **asdf_open_memory_mapping_kwarg(memmap=False)) as af:
         miriwcs = af.tree['wcs']
     bounding_box = ((0, 1024), (0, 1024))
     mirisip = miriwcs.to_fits_sip(bounding_box, max_inv_pix_error=0.1, verbose=True)
@@ -1009,7 +1018,7 @@ def test_to_fits_tab_time_cube(gwcs_cube_with_separable_time):
 def test_to_fits_tab_miri_image():
     # gWCS:
     fn = os.path.join(data_path, 'miriwcs.asdf')
-    with asdf.open(fn, copy_arrays=True, lazy_load=False, ignore_missing_extensions=True) as af:
+    with asdf.open(fn, lazy_load=False, ignore_missing_extensions=True, **asdf_open_memory_mapping_kwarg(memmap=False)) as af:
         w = af.tree['wcs']
 
     # FITS WCS -TAB:
@@ -1033,7 +1042,7 @@ def test_to_fits_tab_miri_image():
 
 def test_to_fits_tab_miri_lrs():
     fn = os.path.join(data_path, 'miri_lrs_wcs.asdf')
-    with asdf.open(fn, copy_arrays=True, lazy_load=False, ignore_missing_extensions=True) as af:
+    with asdf.open(fn, lazy_load=False, ignore_missing_extensions=True, **asdf_open_memory_mapping_kwarg(memmap=False)) as af:
         w = af.tree['wcs']
 
     # FITS WCS -TAB:
@@ -1102,7 +1111,7 @@ def test_in_image():
 
 def test_iter_inv():
     fn = os.path.join(data_path, 'nircamwcs.asdf')
-    with asdf.open(fn, lazy_load=False, copy_arrays=False, ignore_missing_extensions=True) as af:
+    with asdf.open(fn, lazy_load=False, ignore_missing_extensions=True, **asdf_open_memory_mapping_kwarg(memmap=True)) as af:
         w = af.tree['wcs']
     # remove analytic/user-supplied inverse:
     w.pipeline[0].transform.inverse = None
@@ -1129,7 +1138,7 @@ def test_iter_inv():
     )
     assert np.allclose((x, y), (xp, yp))
 
-    with asdf.open(fn, lazy_load=False, copy_arrays=False, ignore_missing_extensions=True) as af:
+    with asdf.open(fn, lazy_load=False, ignore_missing_extensions=True, **asdf_open_memory_mapping_kwarg(memmap=True)) as af:
         w = af.tree['wcs']
 
     # test single point
diff --git a/gwcs/wcs.py b/gwcs/wcs.py
index 0432e541..62173532 100644
--- a/gwcs/wcs.py
+++ b/gwcs/wcs.py
@@ -679,7 +679,7 @@ def numerical_inverse(self, *args, tolerance=1e-5, maxiter=50, adaptive=True,
         >>> import numpy as np
 
         >>> filename = get_pkg_data_filename('data/nircamwcs.asdf', package='gwcs.tests')
-        >>> with asdf.open(filename, copy_arrays=True, lazy_load=False, ignore_missing_extensions=True) as af:
+        >>> with asdf.open(filename, memmap=False, lazy_load=False, ignore_missing_extensions=True) as af:
         ...    w = af.tree['wcs']
 
         >>> ra, dec = w([1,2,3], [1,1,1])