diff --git a/autotest/t018_test.py b/autotest/t018_test.py
index 7410b1f277..bce9061831 100644
--- a/autotest/t018_test.py
+++ b/autotest/t018_test.py
@@ -26,7 +26,7 @@ def test_tpl_constant():
     mfpackage = 'lpf'
     partype = 'hk'
     parname = 'HK_LAYER_1'
-    idx = np.empty((nlay, nrow, ncol), dtype=np.bool)
+    idx = np.empty((nlay, nrow, ncol), dtype=bool)
     idx[0] = True
     idx[1:] = False
 
diff --git a/autotest/t020_test.py b/autotest/t020_test.py
index 2d4d754751..84d522b483 100644
--- a/autotest/t020_test.py
+++ b/autotest/t020_test.py
@@ -9,7 +9,7 @@
 
 
 def analyticalWaterTableSolution(h1, h2, z, R, K, L, x):
-    h = np.zeros((x.shape[0]), np.float)
+    h = np.zeros((x.shape[0]), float)
     b1 = h1 - z
     b2 = h2 - z
     h = np.sqrt(b1 ** 2 - (x / L) * (b1 ** 2 - b2 ** 2) + (R * x / K) * (L - x)) + z
@@ -46,10 +46,10 @@ def test_mfnwt_run():
     h2 = 11.
 
     # ibound
-    ibound = np.ones((nlay, nrow, ncol), dtype=np.int)
+    ibound = np.ones((nlay, nrow, ncol), dtype=int)
 
     # starting heads
-    strt = np.zeros((nlay, nrow, ncol), dtype=np.float)
+    strt = np.zeros((nlay, nrow, ncol), dtype=float)
     strt[0, 0, 0] = h1
     strt[0, 0, -1] = h2
 
diff --git a/autotest/t031_test.py b/autotest/t031_test.py
index 16d864529b..f64adc596c 100644
--- a/autotest/t031_test.py
+++ b/autotest/t031_test.py
@@ -55,17 +55,17 @@ def test_mpsim():
     node_data = np.array(
         [(3, 12, 12, 'well1', 'skin', -1, 0, 0, 0, 1., 2., 5., 6.2),
          (4, 12, 12, 'well1', 'skin', -1, 0, 0, 0, 0.5, 2., 5., 6.2)],
-        dtype=[('k', np.int), ('i', np.int), ('j', np.int),
-               ('wellid', np.object), ('losstype', np.object),
-               ('pumploc', np.int), ('qlimit', np.int),
-               ('ppflag', np.int), ('pumpcap', np.int),
-               ('rw', np.float), ('rskin', np.float),
-               ('kskin', np.float), ('zpump', np.float)]).view(np.recarray)
+        dtype=[('k', int), ('i', int), ('j', int),
+               ('wellid', object), ('losstype', object),
+               ('pumploc', int), ('qlimit', int),
+               ('ppflag', int), ('pumpcap', int),
+               ('rw', float), ('rskin', float),
+               ('kskin', float), ('zpump', float)]).view(np.recarray)
 
     stress_period_data = {0: np.array([(0, 'well1', -150000.0)],
-                                      dtype=[('per', np.int),
-                                             ('wellid', np.object),
-                                             ('qdes', np.float)])}
+                                      dtype=[('per', int),
+                                             ('wellid', object),
+                                             ('qdes', float)])}
     m.remove_package('WEL')
     mnw2 = flopy.modflow.ModflowMnw2(model=m, mnwmax=1,
                                      node_data=node_data,
diff --git a/autotest/t032_test.py b/autotest/t032_test.py
index 24e502b23b..a7c6239436 100644
--- a/autotest/t032_test.py
+++ b/autotest/t032_test.py
@@ -48,7 +48,7 @@ def test_polygon_from_ij():
                          (0, 7, 8, .3, True, 's2')],
                         dtype=[('k', '<i8'), ('i', '<i8'), ('j', '<i8'),
                                ('stuff', '<f4'), ('stuf', '|b1'),
-                               ('stf', np.object)]).view(np.recarray)
+                               ('stf', object)]).view(np.recarray)
 
     # vertices for a model cell
     geoms = [Polygon(m.modelgrid.get_cell_vertices(i, j)) for i, j in
diff --git a/autotest/t039_test.py b/autotest/t039_test.py
index 8971a29bdf..fdee7bf244 100644
--- a/autotest/t039_test.py
+++ b/autotest/t039_test.py
@@ -61,9 +61,9 @@ def read_zonebudget_file(fname):
         vals = [float(i) for i in items[1:-1]]
         row = (totim, kstp, kper, record,) + tuple(v for v in vals)
         rows.append(row)
-    dtype_list = [('totim', np.float),
-                  ('time_step', np.int),
-                  ('stress_period', np.int),
+    dtype_list = [('totim', float),
+                  ('time_step', int),
+                  ('stress_period', int),
                   ('name', '<U50')] + [(z, '<f8') for z in zonenames]
     dtype = np.dtype(dtype_list)
     return np.array(rows, dtype=dtype)
@@ -120,7 +120,7 @@ def test_compare2zonebudget(rtol=1e-2):
 #     mflistrecs = mflistbud.get_data(idx=-1, incremental=True)
 #     print(repr(mflistrecs))
 #
-#     zon = np.ones((3, 40, 20), dtype=np.int)
+#     zon = np.ones((3, 40, 20), dtype=int)
 #     cbc_fname = os.path.join(loadpth, 'freyberg.cbc')
 #     kstp, kper = CellBudgetFile(cbc_fname).get_kstpkper()[-1]
 #     zb = ZoneBudget(cbc_fname, zon, kstpkper=(kstp, kper))
diff --git a/autotest/t040_test.py b/autotest/t040_test.py
index 899050d876..3d02274d2c 100644
--- a/autotest/t040_test.py
+++ b/autotest/t040_test.py
@@ -97,7 +97,7 @@ def test_subwt():
                              vka=[0.4, 0.4, 0.01, 0.4],
                              sy=0.3, ss=1.0e-6)
 
-    # temp_ib = np.ones((ml.nrow,ml.ncol),dtype=np.int)
+    # temp_ib = np.ones((ml.nrow,ml.ncol),dtype=int)
     # np.savetxt('temp_ib.dat',temp_ib,fmt='%1d')
     ibound = np.loadtxt(ibound_path)
     ibound[ibound == 5] = -1
diff --git a/autotest/t041_test.py b/autotest/t041_test.py
index f8fb7b574c..341190ff52 100644
--- a/autotest/t041_test.py
+++ b/autotest/t041_test.py
@@ -34,7 +34,7 @@ def test_hob_simple():
     nlay, nrow, ncol = 1, 11, 11
     shape3d = (nlay, nrow, ncol)
     shape2d = (nrow, ncol)
-    ib = np.ones(shape3d, dtype=np.int)
+    ib = np.ones(shape3d, dtype=int)
     ib[0, 0, 0] = -1
     m = flopy.modflow.Modflow(modelname=modelname, model_ws=pth,
                               verbose=False, exe_name=exe_name, )
diff --git a/autotest/t047_test.py b/autotest/t047_test.py
index b58133ecd9..1fb449e2fa 100644
--- a/autotest/t047_test.py
+++ b/autotest/t047_test.py
@@ -20,7 +20,7 @@ def get_namefile_entries(fpth):
     except:
         print('could not open...{}'.format(fpth))
         return None
-    dtype = [('ftype', '|S12'), ('unit', np.int),
+    dtype = [('ftype', '|S12'), ('unit', int),
              ('filename', '|S128'), ('status', '|S10')]
     lines = f.readlines()
     data = []
diff --git a/autotest/t051_test.py b/autotest/t051_test.py
index ed7bd85a40..776a27088b 100644
--- a/autotest/t051_test.py
+++ b/autotest/t051_test.py
@@ -50,7 +50,7 @@ def test_mfcbc():
     ml = flopy.modflow.Modflow(modelname='t1', model_ws=cpth, verbose=True)
     dis = flopy.modflow.ModflowDis(ml, nlay=nlay, nrow=nrow, ncol=ncol, top=0,
                                    botm=[-1., -2., -3.])
-    ibound = np.ones((nlay, nrow, ncol), dtype=np.int)
+    ibound = np.ones((nlay, nrow, ncol), dtype=int)
     ibound[0, 1, 1] = 0
     ibound[0, 0, -1] = -1
     bas = flopy.modflow.ModflowBas(ml, ibound=ibound)
diff --git a/autotest/t052_test.py b/autotest/t052_test.py
index be8986003f..5b5424d490 100644
--- a/autotest/t052_test.py
+++ b/autotest/t052_test.py
@@ -33,7 +33,7 @@ def test_binary_well():
                                exe_name=exe_name)
     dis = flopy.modflow.ModflowDis(ml, nlay=nlay, nrow=nrow, ncol=ncol, top=0,
                                    botm=[-1., -2., -3.])
-    ibound = np.ones((nlay, nrow, ncol), dtype=np.int)
+    ibound = np.ones((nlay, nrow, ncol), dtype=int)
     ibound[0, 1, 1] = 0
     ibound[0, 0, -1] = -1
     bas = flopy.modflow.ModflowBas(ml, ibound=ibound)
diff --git a/autotest/t063_test_lgrutil.py b/autotest/t063_test_lgrutil.py
index 98400f53a5..528464f9a9 100644
--- a/autotest/t063_test_lgrutil.py
+++ b/autotest/t063_test_lgrutil.py
@@ -18,7 +18,7 @@ def test_lgrutil():
     delcp = 100.
     topp = 100.
     botmp = [-100, -200, -300, -400, -500]
-    idomainp = np.ones((nlayp, nrowp, ncolp), dtype=np.int)
+    idomainp = np.ones((nlayp, nrowp, ncolp), dtype=int)
     idomainp[0:2, 1:4, 1:4] = 0
     ncpp = 3
     ncppl = [1, 1, 0, 0, 0]
diff --git a/autotest/t065_test_gridintersect.py b/autotest/t065_test_gridintersect.py
index 118f957218..717bccc596 100644
--- a/autotest/t065_test_gridintersect.py
+++ b/autotest/t065_test_gridintersect.py
@@ -51,8 +51,8 @@ def get_tri_grid(angrot=0., xyoffset=0., triangle_exe=None):
 
 
 def get_rect_grid(angrot=0., xyoffset=0., top=None, botm=None):
-    delc = 10 * np.ones(2, dtype=np.float)
-    delr = 10 * np.ones(2, dtype=np.float)
+    delc = 10 * np.ones(2, dtype=float)
+    delr = 10 * np.ones(2, dtype=float)
     sgr = fgrid.StructuredGrid(
         delc, delr, top=top, botm=botm, xoff=xyoffset, yoff=xyoffset,
         angrot=angrot)
diff --git a/autotest/t068_test_ssm.py b/autotest/t068_test_ssm.py
index c4ed09483d..dd2abc7a72 100644
--- a/autotest/t068_test_ssm.py
+++ b/autotest/t068_test_ssm.py
@@ -18,10 +18,10 @@
 def test_mt3d_ssm_with_nodata_in_1st_sp():
 
     nlay, nrow, ncol = 3, 5, 5
-    perlen = np.zeros((10), dtype=np.float) + 10
+    perlen = np.zeros((10), dtype=float) + 10
     nper = len(perlen)
     
-    ibound = np.ones((nlay,nrow,ncol), dtype=np.int)
+    ibound = np.ones((nlay,nrow,ncol), dtype=int)
     
     botm = np.arange(-1,-4,-1)
     top = 0.
diff --git a/autotest/t075_ugridtests.py b/autotest/t075_ugridtests.py
index 6f090cbe28..8ed64558b9 100644
--- a/autotest/t075_ugridtests.py
+++ b/autotest/t075_ugridtests.py
@@ -56,7 +56,7 @@ def test_unstructured_minimal_grid():
     g = UnstructuredGrid(
         vertices=vertices, iverts=iverts, xcenters=xcenters, ycenters=ycenters
     )
-    assert np.allclose(g.ncpl, np.array([2], dtype=np.int))
+    assert np.allclose(g.ncpl, np.array([2], dtype=int))
     assert g.nlay == 1
     assert g.nnodes == 2
     assert g.is_valid
@@ -117,7 +117,7 @@ def test_unstructured_complete_grid():
         top=top,
         botm=botm,
     )
-    assert np.allclose(g.ncpl, np.array([1, 1], dtype=np.int))
+    assert np.allclose(g.ncpl, np.array([1, 1], dtype=int))
     assert g.nlay == 2
     assert g.nnodes == 2
     assert g.is_valid
diff --git a/autotest/t505_test.py b/autotest/t505_test.py
index 14888d3043..0c337103ed 100644
--- a/autotest/t505_test.py
+++ b/autotest/t505_test.py
@@ -870,7 +870,7 @@ def test021_twri():
     )
 
     # build stress_period_data for drn package
-    conc = np.ones((15, 15), dtype=np.float) * 35.0
+    conc = np.ones((15, 15), dtype=float) * 35.0
     auxdata = {0: [6, conc]}
 
     stress_period_data = []
diff --git a/examples/Notebooks/flopy3_MT3DMS_examples.ipynb b/examples/Notebooks/flopy3_MT3DMS_examples.ipynb
index 70ab4b5908..16a8dc1300 100644
--- a/examples/Notebooks/flopy3_MT3DMS_examples.ipynb
+++ b/examples/Notebooks/flopy3_MT3DMS_examples.ipynb
@@ -111,10 +111,10 @@
     "    dis = flopy.modflow.ModflowDis(mf, nlay=nlay, nrow=nrow, ncol=ncol,\n",
     "                                   delr=delr, delc=delc, top=0., botm=[0 - delv],\n",
     "                                   perlen=perlen)\n",
-    "    ibound = np.ones((nlay, nrow, ncol), dtype=np.int)\n",
+    "    ibound = np.ones((nlay, nrow, ncol), dtype=int)\n",
     "    ibound[0, 0, 0] = -1\n",
     "    ibound[0, 0, -1] = -1\n",
-    "    strt = np.zeros((nlay, nrow, ncol), dtype=np.float)\n",
+    "    strt = np.zeros((nlay, nrow, ncol), dtype=float)\n",
     "    h1 = q * Lx\n",
     "    strt[0, 0, 0] = h1\n",
     "    bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)\n",
@@ -128,9 +128,9 @@
     "    mt = flopy.mt3d.Mt3dms(modelname=modelname_mt, model_ws=model_ws, \n",
     "                           exe_name=exe_name_mt, modflowmodel=mf)\n",
     "    c0 = 1.\n",
-    "    icbund = np.ones((nlay, nrow, ncol), dtype=np.int)\n",
+    "    icbund = np.ones((nlay, nrow, ncol), dtype=int)\n",
     "    icbund[0, 0, 0] = -1\n",
-    "    sconc = np.zeros((nlay, nrow, ncol), dtype=np.float)\n",
+    "    sconc = np.zeros((nlay, nrow, ncol), dtype=float)\n",
     "    sconc[0, 0, 0] = c0\n",
     "    btn = flopy.mt3d.Mt3dBtn(mt, icbund=icbund, prsity=prsity, sconc=sconc)\n",
     "    dceps = 1.e-5\n",
@@ -358,10 +358,10 @@
     "    dis = flopy.modflow.ModflowDis(mf, nlay=nlay, nrow=nrow, ncol=ncol,\n",
     "                                   delr=delr, delc=delc, top=0., botm=[0 - delv],\n",
     "                                   perlen=perlen_mf)\n",
-    "    ibound = np.ones((nlay, nrow, ncol), dtype=np.int)\n",
+    "    ibound = np.ones((nlay, nrow, ncol), dtype=int)\n",
     "    ibound[0, 0, 0] = -1\n",
     "    ibound[0, 0, -1] = -1\n",
-    "    strt = np.zeros((nlay, nrow, ncol), dtype=np.float)\n",
+    "    strt = np.zeros((nlay, nrow, ncol), dtype=float)\n",
     "    h1 = q * Lx\n",
     "    strt[0, 0, 0] = h1\n",
     "    bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)\n",
@@ -580,10 +580,10 @@
     "    dis = flopy.modflow.ModflowDis(mf, nlay=nlay, nrow=nrow, ncol=ncol,\n",
     "                                   delr=delr, delc=delc, top=0., botm=[0 - delv],\n",
     "                                   perlen=perlen_mf)\n",
-    "    ibound = np.ones((nlay, nrow, ncol), dtype=np.int)\n",
+    "    ibound = np.ones((nlay, nrow, ncol), dtype=int)\n",
     "    ibound[0, :, 0] = -1\n",
     "    ibound[0, :, -1] = -1\n",
-    "    strt = np.zeros((nlay, nrow, ncol), dtype=np.float)\n",
+    "    strt = np.zeros((nlay, nrow, ncol), dtype=float)\n",
     "    h1 = q * Lx\n",
     "    strt[0, :, 0] = h1\n",
     "    bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)\n",
@@ -761,7 +761,7 @@
     "    dis = flopy.modflow.ModflowDis(mf, nlay=nlay, nrow=nrow, ncol=ncol,\n",
     "                                   delr=delr, delc=delc, top=0., botm=[0 - delv],\n",
     "                                   perlen=perlen_mf)\n",
-    "    ibound = np.ones((nlay, nrow, ncol), dtype=np.int) * -1\n",
+    "    ibound = np.ones((nlay, nrow, ncol), dtype=int) * -1\n",
     "    ibound[:, 1:nrow - 1, 1:ncol-1] = 1\n",
     "\n",
     "    # set strt as a linear gradient at a 45 degree angle\n",
@@ -995,7 +995,7 @@
     "    dis = flopy.modflow.ModflowDis(mf, nlay=nlay, nrow=nrow, ncol=ncol,\n",
     "                                   delr=delr, delc=delc, top=0., botm=[0 - delv],\n",
     "                                   perlen=perlen_mf)\n",
-    "    ibound = np.ones((nlay, nrow, ncol), dtype=np.int) * -1\n",
+    "    ibound = np.ones((nlay, nrow, ncol), dtype=int) * -1\n",
     "    ibound[:, 1:nrow - 1, 1:ncol-1] = 1\n",
     "    strt = 0.\n",
     "\n",
@@ -1183,7 +1183,7 @@
     "    dis = flopy.modflow.ModflowDis(mf, nlay=nlay, nrow=nrow, ncol=ncol,\n",
     "                                   delr=delr, delc=delc, top=0., botm=[0 - delv],\n",
     "                                   nper=2, perlen=perlen_mf)\n",
-    "    ibound = np.ones((nlay, nrow, ncol), dtype=np.int) * -1\n",
+    "    ibound = np.ones((nlay, nrow, ncol), dtype=int) * -1\n",
     "    ibound[:, 1:nrow - 1, 1:ncol-1] = 1\n",
     "    strt = 0.\n",
     "\n",
@@ -1337,10 +1337,10 @@
     "                                   delr=delr, delc=delc, top=0., \n",
     "                                   botm=[-delv * k for k in range(1, nlay + 1)],\n",
     "                                   perlen=perlen_mf)\n",
-    "    ibound = np.ones((nlay, nrow, ncol), dtype=np.int)\n",
+    "    ibound = np.ones((nlay, nrow, ncol), dtype=int)\n",
     "    ibound[:, :, 0] = -1\n",
     "    ibound[:, :, -1] = -1\n",
-    "    strt = np.zeros((nlay, nrow, ncol), dtype=np.float)\n",
+    "    strt = np.zeros((nlay, nrow, ncol), dtype=float)\n",
     "    h1 = q * Lx\n",
     "    strt[:, :, 0] = h1\n",
     "    bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)\n",
@@ -1487,7 +1487,7 @@
     "\n",
     "    k1 = 5e-4 / 100. * 86400 # m/d\n",
     "    k2 = 1e-2 / 100. * 86400 # m/d\n",
-    "    hk = k1 * np.ones((nlay, nrow, ncol), dtype=np.float)\n",
+    "    hk = k1 * np.ones((nlay, nrow, ncol), dtype=float)\n",
     "    hk[11:19, :, 0:24] = k2\n",
     "    hk[11:19, :, 36:] = k2\n",
     "    laytyp = 6 * [1] + 21 * [0]\n",
@@ -1500,10 +1500,10 @@
     "                                   botm=[6.75 - delv * k for k in range(1, nlay + 1)],\n",
     "                                   perlen=perlen_mf)\n",
     "    f = open(os.path.join(datadir, 'p08shead.dat'))\n",
-    "    strt = np.empty((nlay * ncol), dtype=np.float)\n",
+    "    strt = np.empty((nlay * ncol), dtype=float)\n",
     "    strt = read1d(f, strt).reshape((nlay, nrow, ncol))\n",
     "    f.close()\n",
-    "    ibound = np.ones((nlay, nrow, ncol), dtype=np.int)\n",
+    "    ibound = np.ones((nlay, nrow, ncol), dtype=int)\n",
     "    ibound[5:, :, -1] = -1\n",
     "    ibound[strt < 0] = 0\n",
     "\n",
@@ -1537,7 +1537,7 @@
     "                             nlsink=nlsink, npsink=npsink, percel=percel, \n",
     "                             itrack=itrack, wd=wd)\n",
     "    dsp = flopy.mt3d.Mt3dDsp(mt, al=al, trpt=trpt, trpv=trpv, dmcoef=dmcoef)\n",
-    "    crch1 = np.zeros((nrow, ncol), dtype=np.float)\n",
+    "    crch1 = np.zeros((nrow, ncol), dtype=float)\n",
     "    crch1[0, 9:18] = 1.\n",
     "    cnc0 = [(0, 0, j, 1, -1) for j in range(8, 16)]\n",
     "    cnc1 = [(0, 0, j, 0., -1) for j in range(8, 16)]\n",
@@ -1662,7 +1662,7 @@
     "    laytyp = 0\n",
     "    k1 = 1.474e-4\n",
     "    k2 = 1.474e-7\n",
-    "    hk = k1 * np.ones((nlay, nrow, ncol), dtype=np.float)\n",
+    "    hk = k1 * np.ones((nlay, nrow, ncol), dtype=float)\n",
     "    hk[:, 5:8, 1:8] = k2\n",
     "\n",
     "    modelname_mf = dirname + '_mf'\n",
@@ -1670,10 +1670,10 @@
     "    dis = flopy.modflow.ModflowDis(mf, nlay=nlay, nrow=nrow, ncol=ncol,\n",
     "                                   delr=delr, delc=delc, top=0., botm=[0 - delv],\n",
     "                                   perlen=perlen_mf)\n",
-    "    ibound = np.ones((nlay, nrow, ncol), dtype=np.int)\n",
+    "    ibound = np.ones((nlay, nrow, ncol), dtype=int)\n",
     "    ibound[0, 0, :] = -1\n",
     "    ibound[0, -1, :] = -1\n",
-    "    strt = np.zeros((nlay, nrow, ncol), dtype=np.float)\n",
+    "    strt = np.zeros((nlay, nrow, ncol), dtype=float)\n",
     "    strt[0, 0, :] = 250.\n",
     "    xc = dis.sr.xcenter\n",
     "    for j in range(ncol):\n",
@@ -1898,7 +1898,7 @@
     "                                   delr=delr, delc=delc, top=top, \n",
     "                                   botm=botm,\n",
     "                                   perlen=perlen_mf)\n",
-    "    ibound = np.ones((nlay, nrow, ncol), dtype=np.int)\n",
+    "    ibound = np.ones((nlay, nrow, ncol), dtype=int)\n",
     "\n",
     "    # left side\n",
     "    ibound[:, :, 0] = -1\n",
@@ -1910,10 +1910,10 @@
     "    ibound[:, -1, :] = -1\n",
     "\n",
     "    f = open(os.path.join(datadir, 'p10shead.dat'))\n",
-    "    s0 = np.empty((nrow * ncol), dtype=np.float)\n",
+    "    s0 = np.empty((nrow * ncol), dtype=float)\n",
     "    s0 = read1d(f, s0).reshape((nrow, ncol))\n",
     "    f.close()\n",
-    "    strt = np.zeros((nlay, nrow, ncol), dtype=np.float)\n",
+    "    strt = np.zeros((nlay, nrow, ncol), dtype=float)\n",
     "    for k in range(nlay):\n",
     "        strt[k] = s0\n",
     "    bas = flopy.modflow.ModflowBas(mf, ibound=ibound, strt=strt)\n",
@@ -1940,10 +1940,10 @@
     "    mt = flopy.mt3d.Mt3dms(modelname=modelname_mt, model_ws=model_ws, \n",
     "                           exe_name=exe_name_mt, modflowmodel=mf)\n",
     "    f = open(os.path.join(datadir, 'p10cinit.dat'))\n",
-    "    c0 = np.empty((nrow * ncol), dtype=np.float)\n",
+    "    c0 = np.empty((nrow * ncol), dtype=float)\n",
     "    c0 = read1d(f, c0).reshape((nrow, ncol))\n",
     "    f.close()\n",
-    "    sconc = np.zeros((nlay, nrow, ncol), dtype=np.float)\n",
+    "    sconc = np.zeros((nlay, nrow, ncol), dtype=float)\n",
     "    sconc[1] = 0.2 * c0\n",
     "    sconc[2] = c0\n",
     "    obs = [[3 - 1, 11 - 1, 29 - 1],\n",
diff --git a/examples/Notebooks/flopy3_PEST.ipynb b/examples/Notebooks/flopy3_PEST.ipynb
index e86dd2a1c4..b6b86b69d5 100644
--- a/examples/Notebooks/flopy3_PEST.ipynb
+++ b/examples/Notebooks/flopy3_PEST.ipynb
@@ -84,7 +84,7 @@
     "mfpackage = 'lpf'\n",
     "partype = 'hk'\n",
     "parname = 'HK_LAYER_1'\n",
-    "idx = np.empty((nlay, nrow, ncol), dtype=np.bool)\n",
+    "idx = np.empty((nlay, nrow, ncol), dtype=bool)\n",
     "idx[0] = True\n",
     "idx[1:] = False\n",
     "\n",
@@ -674,7 +674,7 @@
     "transform = None\n",
     "\n",
     "# For a recharge index parameter, span['idx'] must be a boolean array or tuple of array indices\n",
-    "idx = np.empty((nrow, ncol), dtype=np.bool)\n",
+    "idx = np.empty((nrow, ncol), dtype=bool)\n",
     "idx[0:3, 0:3] = True\n",
     "span = {'kpers': [1], 'idx': idx}\n",
     "p = flopy.pest.Params(mfpackage, partype, parname, startvalue,\n",
diff --git a/examples/Notebooks/flopy3_WatertableRecharge_example.ipynb b/examples/Notebooks/flopy3_WatertableRecharge_example.ipynb
index b999636950..9a64b33915 100644
--- a/examples/Notebooks/flopy3_WatertableRecharge_example.ipynb
+++ b/examples/Notebooks/flopy3_WatertableRecharge_example.ipynb
@@ -91,7 +91,7 @@
    "outputs": [],
    "source": [
     "def analyticalWaterTableSolution(h1, h2, z, R, K, L, x):\n",
-    "    h = np.zeros((x.shape[0]), np.float)\n",
+    "    h = np.zeros((x.shape[0]), float)\n",
     "    #dx = x[1] - x[0]\n",
     "    #x -= dx\n",
     "    b1 = h1 - z\n",
@@ -138,10 +138,10 @@
     "h2 = 11.\n",
     "\n",
     "# ibound\n",
-    "ibound = np.ones((nlay, nrow, ncol), dtype=np.int)\n",
+    "ibound = np.ones((nlay, nrow, ncol), dtype=int)\n",
     "\n",
     "# starting heads\n",
-    "strt = np.zeros((nlay, nrow, ncol), dtype=np.float)\n",
+    "strt = np.zeros((nlay, nrow, ncol), dtype=float)\n",
     "strt[0, 0, 0] = h1\n",
     "strt[0, 0, -1] = h2\n",
     "\n",
diff --git a/examples/Notebooks/flopy3_Zaidel_example.ipynb b/examples/Notebooks/flopy3_Zaidel_example.ipynb
index 2c3b340d13..83522e1be9 100644
--- a/examples/Notebooks/flopy3_Zaidel_example.ipynb
+++ b/examples/Notebooks/flopy3_Zaidel_example.ipynb
@@ -103,12 +103,12 @@
     "x = np.arange(0., float(ncol)*delr, delr) + delr / 2.\n",
     "\n",
     "# ibound\n",
-    "ibound = np.ones((nlay, nrow, ncol), dtype=np.int)\n",
+    "ibound = np.ones((nlay, nrow, ncol), dtype=int)\n",
     "ibound[:, :, 0] = -1\n",
     "ibound[:, :, -1] = -1\n",
     "\n",
     "# bottom of the model\n",
-    "botm = 25 * np.ones((nlay + 1, nrow, ncol), dtype=np.float)\n",
+    "botm = 25 * np.ones((nlay + 1, nrow, ncol), dtype=float)\n",
     "base = 20.\n",
     "for j in range(ncol):\n",
     "    botm[1, :, j] = base\n",
@@ -117,7 +117,7 @@
     "        base -= 5\n",
     "\n",
     "# starting heads\n",
-    "strt = h1 * np.ones((nlay, nrow, ncol), dtype=np.float)\n",
+    "strt = h1 * np.ones((nlay, nrow, ncol), dtype=float)\n",
     "strt[:, :, -1] = h2"
    ]
   },
diff --git a/examples/Notebooks/flopy3_ZoneBudget_example.ipynb b/examples/Notebooks/flopy3_ZoneBudget_example.ipynb
index ad40ec297e..10f17491b3 100644
--- a/examples/Notebooks/flopy3_ZoneBudget_example.ipynb
+++ b/examples/Notebooks/flopy3_ZoneBudget_example.ipynb
@@ -1137,7 +1137,7 @@
     }
    ],
    "source": [
-    "zon = np.ones((nlay, nrow, ncol), np.int)\n",
+    "zon = np.ones((nlay, nrow, ncol), int)\n",
     "zon[1, :, :] = 2\n",
     "zon[2, :, :] = 3\n",
     "\n",
diff --git a/examples/Notebooks/flopy3_demo_of_modelgrid_classes.ipynb b/examples/Notebooks/flopy3_demo_of_modelgrid_classes.ipynb
index c03007187b..48cce5da55 100644
--- a/examples/Notebooks/flopy3_demo_of_modelgrid_classes.ipynb
+++ b/examples/Notebooks/flopy3_demo_of_modelgrid_classes.ipynb
@@ -435,8 +435,8 @@
     "dis = flopy.modflow.ModflowDis(ml,nlay=5,nrow=5,ncol=5,delr=[50,50,50,50,50],\n",
     "                               delc=[50,50,50,50,50],top=0.0,botm=-50.0,laycbd=0) \n",
     "\n",
-    "ibound = np.ones((5, 5, 5), dtype=np.int)\n",
-    "start = np.zeros((5, 5, 5), dtype=np.float)\n",
+    "ibound = np.ones((5, 5, 5), dtype=int)\n",
+    "start = np.zeros((5, 5, 5), dtype=float)\n",
     "bas = flopy.modflow.ModflowBas(ml,ibound=ibound,strt=start)"
    ]
   },
diff --git a/examples/Notebooks/flopy3_grid_intersection_demo.ipynb b/examples/Notebooks/flopy3_grid_intersection_demo.ipynb
index 728f2c9c27..d780a94ea0 100644
--- a/examples/Notebooks/flopy3_grid_intersection_demo.ipynb
+++ b/examples/Notebooks/flopy3_grid_intersection_demo.ipynb
@@ -124,8 +124,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "delc = 10*np.ones(10, dtype=np.float)\n",
-    "delr = 10*np.ones(10, dtype=np.float)"
+    "delc = 10*np.ones(10, dtype=float)\n",
+    "delr = 10*np.ones(10, dtype=float)"
    ]
   },
   {
diff --git a/examples/Notebooks/flopy3_gridgen.ipynb b/examples/Notebooks/flopy3_gridgen.ipynb
index f48710fa14..c758cee2a1 100644
--- a/examples/Notebooks/flopy3_gridgen.ipynb
+++ b/examples/Notebooks/flopy3_gridgen.ipynb
@@ -472,7 +472,7 @@
     }
    ],
    "source": [
-    "a = np.zeros((g.nodes), dtype=np.int)\n",
+    "a = np.zeros((g.nodes), dtype=int)\n",
     "a[adpoly_intersect.nodenumber] = 1\n",
     "a[well_intersect.nodenumber] = 2\n",
     "a[river_intersect.nodenumber] = 3\n",
diff --git a/examples/Notebooks/flopy3_lake_example.ipynb b/examples/Notebooks/flopy3_lake_example.ipynb
index 74a889650e..880e710ac6 100644
--- a/examples/Notebooks/flopy3_lake_example.ipynb
+++ b/examples/Notebooks/flopy3_lake_example.ipynb
@@ -125,7 +125,7 @@
    "outputs": [],
    "source": [
     "Nhalf = int((N-1) / 2) \n",
-    "ibound = np.ones((Nlay, N, N), dtype=np.int) \n",
+    "ibound = np.ones((Nlay, N, N), dtype=int) \n",
     "ibound[:,0,:] = -1; ibound[:,-1,:] = -1; ibound[:,:,0] = -1; ibound[:,:,-1] = -1 \n",
     "ibound[0,Nhalf,Nhalf] = -1 \n",
     "start = h1 * np.ones((N,N)) \n",
diff --git a/examples/Notebooks/flopy3_mf6_A_simple-model.ipynb b/examples/Notebooks/flopy3_mf6_A_simple-model.ipynb
index 232b549bbc..d0fe58f7b0 100644
--- a/examples/Notebooks/flopy3_mf6_A_simple-model.ipynb
+++ b/examples/Notebooks/flopy3_mf6_A_simple-model.ipynb
@@ -263,7 +263,7 @@
     "# heads are located by creating an integer array\n",
     "# that starts with ones everywhere, but is assigned\n",
     "# a -1 where chds are located\n",
-    "ibd = np.ones((Nlay, N, N), dtype=np.int)\n",
+    "ibd = np.ones((Nlay, N, N), dtype=int)\n",
     "for k, i, j in ra['cellid']:\n",
     "    ibd[k, i, j] = -1\n",
     "\n",
diff --git a/examples/Notebooks/flopy3_modpath6_example.ipynb b/examples/Notebooks/flopy3_modpath6_example.ipynb
index 18fb49af5f..b430987874 100644
--- a/examples/Notebooks/flopy3_modpath6_example.ipynb
+++ b/examples/Notebooks/flopy3_modpath6_example.ipynb
@@ -623,15 +623,15 @@
    "source": [
     "node_data = np.array([(3, 12, 12, 'well1', 'skin', -1, 0, 0, 0, 1., 2., 5., 6.2),\n",
     "                      (4, 12, 12, 'well1', 'skin', -1, 0, 0, 0, 0.5, 2., 5., 6.2)], \n",
-    "                     dtype=[('k', np.int), ('i', np.int), ('j', np.int), \n",
-    "                            ('wellid', np.object), ('losstype', np.object), \n",
-    "                            ('pumploc', np.int), ('qlimit', np.int), \n",
-    "                            ('ppflag', np.int), ('pumpcap', np.int), \n",
-    "                            ('rw', np.float), ('rskin', np.float), \n",
-    "                            ('kskin', np.float), ('zpump', np.float)]).view(np.recarray)\n",
+    "                     dtype=[('k', int), ('i', int), ('j', int), \n",
+    "                            ('wellid', object), ('losstype', object), \n",
+    "                            ('pumploc', int), ('qlimit', int), \n",
+    "                            ('ppflag', int), ('pumpcap', int), \n",
+    "                            ('rw', float), ('rskin', float), \n",
+    "                            ('kskin', float), ('zpump', float)]).view(np.recarray)\n",
     "\n",
-    "stress_period_data = {0: np.array([(0, 'well1', -150000.0)], dtype=[('per', np.int), ('wellid', np.object), \n",
-    "                                                            ('qdes', np.float)])}"
+    "stress_period_data = {0: np.array([(0, 'well1', -150000.0)], dtype=[('per', int), ('wellid', object), \n",
+    "                                                            ('qdes', float)])}"
    ]
   },
   {
diff --git a/examples/Notebooks/flopy3_modpath7_unstructured_example.ipynb b/examples/Notebooks/flopy3_modpath7_unstructured_example.ipynb
index 4d75685487..6f5f2040c9 100644
--- a/examples/Notebooks/flopy3_modpath7_unstructured_example.ipynb
+++ b/examples/Notebooks/flopy3_modpath7_unstructured_example.ipynb
@@ -361,7 +361,7 @@
     "                              stress_period_data=welspd)\n",
     "\n",
     "# rch\n",
-    "aux = [np.ones(ncpl, dtype=np.int)*6]\n",
+    "aux = [np.ones(ncpl, dtype=int)*6]\n",
     "rch = flopy.mf6.ModflowGwfrcha(gwf, recharge=0.005, \n",
     "                               auxiliary=[('iface',)], \n",
     "                               aux={0: [6]})\n",
@@ -506,7 +506,7 @@
     }
    ],
    "source": [
-    "ibd = np.zeros((ncpl), dtype=np.int)\n",
+    "ibd = np.zeros((ncpl), dtype=int)\n",
     "ibd[welcells['nodenumber']] = 1\n",
     "ibd[rivcells['nodenumber']] = 2\n",
     "ibd = np.ma.masked_equal(ibd, 0)\n",
diff --git a/examples/Notebooks/flopy3_multi-component_SSM.ipynb b/examples/Notebooks/flopy3_multi-component_SSM.ipynb
index df369ab8d5..5d6c5b411a 100644
--- a/examples/Notebooks/flopy3_multi-component_SSM.ipynb
+++ b/examples/Notebooks/flopy3_multi-component_SSM.ipynb
@@ -59,10 +59,10 @@
    "outputs": [],
    "source": [
     "nlay, nrow, ncol = 10, 10, 10\n",
-    "perlen = np.zeros((10), dtype=np.float) + 10\n",
+    "perlen = np.zeros((10), dtype=float) + 10\n",
     "nper = len(perlen)\n",
     "\n",
-    "ibound = np.ones((nlay,nrow,ncol), dtype=np.int)\n",
+    "ibound = np.ones((nlay,nrow,ncol), dtype=int)\n",
     "\n",
     "botm = np.arange(-1,-11,-1)\n",
     "top = 0."
diff --git a/examples/Notebooks/flopy3_swi2package_ex1.ipynb b/examples/Notebooks/flopy3_swi2package_ex1.ipynb
index dc12c6e6ca..435cf525cf 100644
--- a/examples/Notebooks/flopy3_swi2package_ex1.ipynb
+++ b/examples/Notebooks/flopy3_swi2package_ex1.ipynb
@@ -201,11 +201,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "z = np.zeros((nrow, ncol), np.float)\n",
+    "z = np.zeros((nrow, ncol), float)\n",
     "z[0, 16:24] = np.arange(-2.5, -40, -5)\n",
     "z[0, 24:] = -40\n",
     "z = [z]  # zeta needs to be \n",
-    "isource = np.ones((nrow, ncol), np.int)\n",
+    "isource = np.ones((nrow, ncol), int)\n",
     "isource[0, 0] = 2\n",
     "#\n",
     "swi = flopy.modflow.ModflowSwi2(ml, nsrf=1, istrat=1, \n",
diff --git a/examples/Notebooks/flopy3_swi2package_ex4.ipynb b/examples/Notebooks/flopy3_swi2package_ex4.ipynb
index ea1ea8ce57..46c2050727 100644
--- a/examples/Notebooks/flopy3_swi2package_ex4.ipynb
+++ b/examples/Notebooks/flopy3_swi2package_ex4.ipynb
@@ -149,13 +149,13 @@
    "source": [
     "# bas data\n",
     "# ibound - active except for the corners\n",
-    "ibound = np.ones((nlay, nrow, ncol), dtype= np.int)\n",
+    "ibound = np.ones((nlay, nrow, ncol), dtype= int)\n",
     "ibound[:, 0, 0] = 0\n",
     "ibound[:, 0, -1] = 0\n",
     "ibound[:, -1, 0] = 0\n",
     "ibound[:, -1, -1] = 0\n",
     "# initial head data\n",
-    "ihead = np.zeros((nlay, nrow, ncol), dtype=np.float)"
+    "ihead = np.zeros((nlay, nrow, ncol), dtype=float)"
    ]
   },
   {
@@ -193,7 +193,7 @@
     "# boundary condition data\n",
     "# ghb data\n",
     "colcell, rowcell = np.meshgrid(np.arange(0, ncol), np.arange(0, nrow))\n",
-    "index = np.zeros((nrow, ncol), dtype=np.int)\n",
+    "index = np.zeros((nrow, ncol), dtype=int)\n",
     "index[:, :10] = 1\n",
     "index[:, -10:] = 1\n",
     "index[:10, :] = 1\n",
@@ -209,7 +209,7 @@
     "ghb_data = {0:lrchc}\n",
     "\n",
     "# recharge data\n",
-    "rch = np.zeros((nrow, ncol), dtype=np.float)\n",
+    "rch = np.zeros((nrow, ncol), dtype=float)\n",
     "rch[index == 0] = 0.0004\n",
     "# create recharge dictionary\n",
     "rch_data = {0: rch}\n",
@@ -245,7 +245,7 @@
     "z1 = -10.0 * np.ones((nrow, ncol))\n",
     "z1[index == 0] = -11.0\n",
     "z = np.array([[z1, z1]])\n",
-    "iso = np.zeros((nlay, nrow, ncol), dtype=np.int)\n",
+    "iso = np.zeros((nlay, nrow, ncol), dtype=int)\n",
     "iso[0, :, :][index == 0] = 1\n",
     "iso[0, :, :][index == 1] = -2\n",
     "iso[1, 30, 35] = 2\n",
diff --git a/examples/Notebooks/flopy3_uzf_example.ipynb b/examples/Notebooks/flopy3_uzf_example.ipynb
index b9cc6436e1..8abe35a8c4 100644
--- a/examples/Notebooks/flopy3_uzf_example.ipynb
+++ b/examples/Notebooks/flopy3_uzf_example.ipynb
@@ -1073,10 +1073,10 @@
     "fpth = os.path.join(path, 'UZFtest2.uzf68.out')\n",
     "avail = os.path.isfile(fpth)\n",
     "if avail:\n",
-    "    dtype = [('TIME', np.float), ('APPLIED-INFIL', np.float), ('RUNOFF', np.float), \n",
-    "             ('ACTUAL-INFIL', np.float), ('SURFACE-LEAK', np.float), \n",
-    "             ('UZ-ET', np.float), ('GW-ET', np.float), ('UZSTOR-CHANGE', np.float), \n",
-    "             ('RECHARGE', np.float)]\n",
+    "    dtype = [('TIME', float), ('APPLIED-INFIL', float), ('RUNOFF', float), \n",
+    "             ('ACTUAL-INFIL', float), ('SURFACE-LEAK', float), \n",
+    "             ('UZ-ET', float), ('GW-ET', float), ('UZSTOR-CHANGE', float), \n",
+    "             ('RECHARGE', float)]\n",
     "    # read data from file\n",
     "    df = np.genfromtxt(fpth, skip_header=3, dtype=dtype)\n",
     "    # convert numpy recarray to pandas dataframe\n",
diff --git a/examples/Notebooks/flopy3_voronoi.ipynb b/examples/Notebooks/flopy3_voronoi.ipynb
index 5fd8d1c887..0f725bf73f 100644
--- a/examples/Notebooks/flopy3_voronoi.ipynb
+++ b/examples/Notebooks/flopy3_voronoi.ipynb
@@ -201,7 +201,7 @@
    "source": [
     "gridprops = voronoi_grid.get_gridprops_vertexgrid()\n",
     "vgrid = flopy.discretization.VertexGrid(**gridprops, nlay=1)\n",
-    "ibd = np.zeros(vgrid.ncpl, dtype=np.int)\n",
+    "ibd = np.zeros(vgrid.ncpl, dtype=int)\n",
     "gi = flopy.utils.GridIntersect(vgrid)\n",
     "\n",
     "# identify cells on left edge\n",
diff --git a/examples/Testing/util_list_tester.py b/examples/Testing/util_list_tester.py
index a110b32ae1..1922ec5785 100644
--- a/examples/Testing/util_list_tester.py
+++ b/examples/Testing/util_list_tester.py
@@ -7,7 +7,7 @@
 
 # instance testing
 #override the default data type
-dtype = np.dtype([("k",int),("i",np.int),("j",np.int),\
+dtype = np.dtype([("k", int), ("i", int), ("j", int),\
                   ("stage",np.float32),("cond",np.float32),("rbot",np.float32),\
     ("aux1",np.float32),("aux2",np.float32),])#("comments",object)])
 
diff --git a/examples/Tutorials/modflow6/tutorial01_mf6.py b/examples/Tutorials/modflow6/tutorial01_mf6.py
index 8bd541d7fd..e069cb1e0b 100644
--- a/examples/Tutorials/modflow6/tutorial01_mf6.py
+++ b/examples/Tutorials/modflow6/tutorial01_mf6.py
@@ -224,7 +224,7 @@
 # Before we start we will create a MODFLOW-2005 ibound array to use to plot
 # the locations of the constant heads.
 
-ibd = np.ones((Nlay, N, N), dtype=np.int)
+ibd = np.ones((Nlay, N, N), dtype=int)
 for k, i, j in ra["cellid"]:
     ibd[k, i, j] = -1
 
diff --git a/examples/Tutorials/seawat/tutorial01_seawat.py b/examples/Tutorials/seawat/tutorial01_seawat.py
index dcde4a75bb..a2969d540e 100644
--- a/examples/Tutorials/seawat/tutorial01_seawat.py
+++ b/examples/Tutorials/seawat/tutorial01_seawat.py
@@ -161,7 +161,7 @@
 cbbobj = bf.CellBudgetFile("henry.cbc")
 times = cbbobj.get_times()
 qx = cbbobj.get_data(text="flow right face", totim=times[-1])[0]
-qy = np.zeros((nlay, nrow, ncol), dtype=np.float)
+qy = np.zeros((nlay, nrow, ncol), dtype=float)
 qz = cbbobj.get_data(text="flow lower face", totim=times[-1])[0]
 
 # ### Create a plot with concentrations and flow vectors
diff --git a/examples/common/setup_pmv_demo.py b/examples/common/setup_pmv_demo.py
index 52b70680ba..ec5c6e3e4c 100644
--- a/examples/common/setup_pmv_demo.py
+++ b/examples/common/setup_pmv_demo.py
@@ -191,7 +191,7 @@ def run():
     )
 
     # rch
-    aux = [np.ones(ncpl, dtype=np.int) * 6]
+    aux = [np.ones(ncpl, dtype=int) * 6]
     rch = flopy.mf6.ModflowGwfrcha(
         gwf, recharge=0.005, auxiliary=[("iface",)], aux={0: [6]}
     )
diff --git a/examples/groundwater_paper/scripts/uspb_capture.py b/examples/groundwater_paper/scripts/uspb_capture.py
index 046ce380d9..a0bdfcbfdf 100644
--- a/examples/groundwater_paper/scripts/uspb_capture.py
+++ b/examples/groundwater_paper/scripts/uspb_capture.py
@@ -24,7 +24,7 @@ def cf_model(model, k, i, j, base, Q=-100):
         os.path.join(cf_pth, "DG.cbc"), precision="double"
     )
     kk = hedObj.get_kstpkper()
-    v = np.zeros((len(kk)), dtype=np.float)
+    v = np.zeros((len(kk)), dtype=float)
     h = hedObj.get_ts((k, i, j))
     for idx, kon in enumerate(kk):
         if h[idx, 1] == model.lpf.hdry:
@@ -116,7 +116,7 @@ def cf_model(model, k, i, j, base, Q=-100):
 )
 
 # create array to store capture fraction data (subset of model)
-cf_array = np.empty((10, nrow2, ncol2), dtype=np.float)
+cf_array = np.empty((10, nrow2, ncol2), dtype=float)
 cf_array.fill(np.nan)
 
 # timer for capture fraction analysis
diff --git a/examples/groundwater_paper/scripts/uspb_capture_par.py b/examples/groundwater_paper/scripts/uspb_capture_par.py
index 5bb06e6879..8bb4e36f1b 100644
--- a/examples/groundwater_paper/scripts/uspb_capture_par.py
+++ b/examples/groundwater_paper/scripts/uspb_capture_par.py
@@ -188,7 +188,7 @@ def cf_model(imod, ion, nmax, k, i, j, Qt, base, hdry):
     line += "  row {} - col {}\n".format(i + 1, j + 1)
     line += "  {}\n".format(elt)
     # get the results
-    v = np.zeros((10), dtype=np.float)
+    v = np.zeros((10), dtype=float)
     if success:
         try:
             hedObj = flopy.utils.HeadFile(
@@ -289,7 +289,7 @@ def doit():
     )
 
     # create array to store capture fraction data (subset of model)
-    cf_array = np.empty((10, nrow2, ncol2), dtype=np.float)
+    cf_array = np.empty((10, nrow2, ncol2), dtype=float)
     cf_array.fill(np.nan)
 
     # timer for capture fraction analysis
diff --git a/examples/scripts/flopy_lake_example.py b/examples/scripts/flopy_lake_example.py
index cb6720aadf..799d536a35 100644
--- a/examples/scripts/flopy_lake_example.py
+++ b/examples/scripts/flopy_lake_example.py
@@ -75,7 +75,7 @@ def run():
     # the other starting values matter). So we set the starting heads to `h1` everywhere, except for
     # the head at the center of the model in the top layer.
     Nhalf = int((N - 1) / 2)
-    ibound = np.ones((Nlay, N, N), dtype=np.int)
+    ibound = np.ones((Nlay, N, N), dtype=int)
     ibound[:, 0, :] = -1
     ibound[:, -1, :] = -1
     ibound[:, :, 0] = -1
diff --git a/examples/scripts/flopy_swi2_ex1.py b/examples/scripts/flopy_swi2_ex1.py
index 5fb36d003a..38f8368a82 100755
--- a/examples/scripts/flopy_swi2_ex1.py
+++ b/examples/scripts/flopy_swi2_ex1.py
@@ -62,16 +62,16 @@ def run():
     delr = 5.0
     delc = 1.0
 
-    ibound = np.ones((nrow, ncol), np.int)
+    ibound = np.ones((nrow, ncol), int)
     ibound[0, -1] = -1
 
     # create initial zeta surface
-    z = np.zeros((nrow, ncol), np.float)
+    z = np.zeros((nrow, ncol), float)
     z[0, 16:24] = np.arange(-2.5, -40, -5)
     z[0, 24:] = -40
     z = [z]
     # create isource for SWI2
-    isource = np.ones((nrow, ncol), np.int)
+    isource = np.ones((nrow, ncol), int)
     isource[0, 0] = 2
 
     ocdict = {}
diff --git a/examples/scripts/flopy_swi2_ex2.py b/examples/scripts/flopy_swi2_ex2.py
index 3c6268238a..b43e3bb8b3 100755
--- a/examples/scripts/flopy_swi2_ex2.py
+++ b/examples/scripts/flopy_swi2_ex2.py
@@ -90,11 +90,11 @@ def run():
     nsurf = 2
     x = np.arange(0.5 * delr, ncol * delr, delr)
     xedge = np.linspace(0, float(ncol) * delr, len(x) + 1)
-    ibound = np.ones((nrow, ncol), np.int)
+    ibound = np.ones((nrow, ncol), int)
     ibound[0, 0] = -1
     # swi2 data
-    z0 = np.zeros((nlay, nrow, ncol), np.float)
-    z1 = np.zeros((nlay, nrow, ncol), np.float)
+    z0 = np.zeros((nlay, nrow, ncol), float)
+    z1 = np.zeros((nlay, nrow, ncol), float)
     z0[0, 0, 30:38] = np.arange(-2.5, -40, -5)
     z0[0, 0, 38:] = -40
     z1[0, 0, 22:30] = np.arange(-2.5, -40, -5)
@@ -217,15 +217,15 @@ def run():
     swt_nrow = 1
     swt_nlay = int(swt_zmax / swt_delz)  # 80
     print(swt_nlay, swt_nrow, swt_ncol)
-    swt_ibound = np.ones((swt_nlay, swt_nrow, swt_ncol), np.int)
+    swt_ibound = np.ones((swt_nlay, swt_nrow, swt_ncol), int)
     # swt_ibound[0, swt_ncol-1, 0] = -1
     swt_ibound[0, 0, 0] = -1
     swt_x = np.arange(0.5 * swt_delr, swt_ncol * swt_delr, swt_delr)
     swt_xedge = np.linspace(0, float(ncol) * delr, len(swt_x) + 1)
     swt_top = 0.0
     z0 = swt_top
-    swt_botm = np.zeros((swt_nlay), np.float)
-    swt_z = np.zeros((swt_nlay), np.float)
+    swt_botm = np.zeros((swt_nlay), float)
+    swt_z = np.zeros((swt_nlay), float)
     zcell = -swt_delz / 2.0
     for ilay in range(0, swt_nlay):
         z0 -= swt_delz
@@ -236,10 +236,10 @@ def run():
     swt_X, swt_Z = np.meshgrid(swt_x, swt_z)
     # mt3d
     # mt3d boundary array set to all active
-    icbund = np.ones((swt_nlay, swt_nrow, swt_ncol), np.int)
+    icbund = np.ones((swt_nlay, swt_nrow, swt_ncol), int)
     # create initial concentrations for MT3D
-    sconc = np.ones((swt_nlay, swt_nrow, swt_ncol), np.float)
-    sconcp = np.zeros((swt_nlay, swt_ncol), np.float)
+    sconc = np.ones((swt_nlay, swt_nrow, swt_ncol), float)
+    sconcp = np.zeros((swt_nlay, swt_ncol), float)
     xsb = 110
     xbf = 150
     for ilay in range(0, swt_nlay):
@@ -346,7 +346,7 @@ def run():
     ukstpkper = uobj.get_kstpkper()
     print(ukstpkper)
     c = uobj.get_data(totim=times[-1])
-    conc = np.zeros((swt_nlay, swt_ncol), np.float)
+    conc = np.zeros((swt_nlay, swt_ncol), float)
     for icol in range(0, swt_ncol):
         for ilay in range(0, swt_nlay):
             conc[ilay, icol] = c[ilay, 0, icol]
diff --git a/examples/scripts/flopy_swi2_ex3.py b/examples/scripts/flopy_swi2_ex3.py
index ea28b3bf21..ef061d0f70 100755
--- a/examples/scripts/flopy_swi2_ex3.py
+++ b/examples/scripts/flopy_swi2_ex3.py
@@ -25,9 +25,9 @@
 
 def MergeData(ndim, zdata, tb):
     sv = 0.05
-    md = np.empty((ndim), np.float)
+    md = np.empty((ndim), float)
     md.fill(np.nan)
-    found = np.empty((ndim), np.bool)
+    found = np.empty((ndim), bool)
     found.fill(False)
     for idx, layer in enumerate(zdata):
         for jdx, z in enumerate(layer):
@@ -102,7 +102,7 @@ def run():
     zini = np.hstack(
         (-9 * np.ones(24), np.arange(-9, -50, -0.5), -50 * np.ones(94))
     )[np.newaxis, :]
-    iso = np.zeros((1, 200), dtype=np.int)
+    iso = np.zeros((1, 200), dtype=int)
     iso[:, :30] = -2
     # model objects
     ml = flopy.modflow.Modflow(
@@ -279,8 +279,8 @@ def run():
     )
     # ghyben herzberg
     zeta1 = -9 - 40.0 * (head[0, 0, :])
-    gbh = np.empty(len(zeta1), np.float)
-    gbho = np.empty(len(zeta1), np.float)
+    gbh = np.empty(len(zeta1), float)
+    gbho = np.empty(len(zeta1), float)
     for idx, z1 in enumerate(zeta1):
         if z1 >= -9.0 or z1 <= -50.0:
             gbh[idx] = np.nan
diff --git a/examples/scripts/flopy_swi2_ex4.py b/examples/scripts/flopy_swi2_ex4.py
index f1a8842dc1..d8125508cf 100644
--- a/examples/scripts/flopy_swi2_ex4.py
+++ b/examples/scripts/flopy_swi2_ex4.py
@@ -112,13 +112,13 @@ def run():
 
     # bas data
     # ibound - active except for the corners
-    ibound = np.ones((nlay, nrow, ncol), dtype=np.int)
+    ibound = np.ones((nlay, nrow, ncol), dtype=int)
     ibound[:, 0, 0] = 0
     ibound[:, 0, -1] = 0
     ibound[:, -1, 0] = 0
     ibound[:, -1, -1] = 0
     # initial head data
-    ihead = np.zeros((nlay, nrow, ncol), dtype=np.float)
+    ihead = np.zeros((nlay, nrow, ncol), dtype=float)
 
     # lpf data
     laytyp = 0
@@ -128,7 +128,7 @@ def run():
     # boundary condition data
     # ghb data
     colcell, rowcell = np.meshgrid(np.arange(0, ncol), np.arange(0, nrow))
-    index = np.zeros((nrow, ncol), dtype=np.int)
+    index = np.zeros((nrow, ncol), dtype=int)
     index[:, :10] = 1
     index[:, -10:] = 1
     index[:10, :] = 1
@@ -144,7 +144,7 @@ def run():
     ghb_data = {0: lrchc}
 
     # recharge data
-    rch = np.zeros((nrow, ncol), dtype=np.float)
+    rch = np.zeros((nrow, ncol), dtype=float)
     rch[index == 0] = 0.0004
     # create recharge dictionary
     rch_data = {0: rch}
@@ -174,7 +174,7 @@ def run():
     z1 = -10.0 * np.ones((nrow, ncol))
     z1[index == 0] = -11.0
     z = np.array([[z1, z1]])
-    iso = np.zeros((nlay, nrow, ncol), np.int)
+    iso = np.zeros((nlay, nrow, ncol), int)
     iso[0, :, :][index == 0] = 1
     iso[0, :, :][index == 1] = -2
     iso[1, 30, 35] = 2
diff --git a/examples/scripts/flopy_swi2_ex5.py b/examples/scripts/flopy_swi2_ex5.py
index ce9c6ae00f..937a0ab738 100755
--- a/examples/scripts/flopy_swi2_ex5.py
+++ b/examples/scripts/flopy_swi2_ex5.py
@@ -83,11 +83,11 @@ def run(silent=False):
     nlay = 6
     nrow = 1
     ncol = 113
-    delr = np.zeros((ncol), np.float)
+    delr = np.zeros((ncol), float)
     delc = 1.0
-    r = np.zeros((ncol), np.float)
-    x = np.zeros((ncol), np.float)
-    edge = np.zeros((ncol), np.float)
+    r = np.zeros((ncol), float)
+    x = np.zeros((ncol), float)
+    edge = np.zeros((ncol), float)
     dx = 25.0
     for i in range(0, ncol):
         delr[i] = dx
@@ -109,10 +109,10 @@ def run(silent=False):
 
     nsave_zeta = 8
     ndecay = 4
-    ibound = np.ones((nlay, nrow, ncol), np.int)
+    ibound = np.ones((nlay, nrow, ncol), int)
     for k in range(0, nlay):
         ibound[k, 0, ncol - 1] = -1
-    bot = np.zeros((nlay, nrow, ncol), np.float)
+    bot = np.zeros((nlay, nrow, ncol), float)
     dz = 100.0 / float(nlay - 1)
     zall = -np.arange(0, 100 + dz, dz)
     zall = np.append(zall, -120.0)
@@ -121,7 +121,7 @@ def run(silent=False):
     for k in range(0, nlay):
         for i in range(0, ncol):
             bot[k, 0, i] = tb[k]
-    isource = np.zeros((nlay, nrow, ncol), np.int)
+    isource = np.zeros((nlay, nrow, ncol), int)
     isource[:, 0, ncol - 1] = 1
     isource[nlay - 1, 0, ncol - 1] = 2
 
@@ -129,10 +129,10 @@ def run(silent=False):
     kvb = (0.0000000000100 * 1000.0 * 9.81 / 0.001) * 60 * 60 * 24
     ssb = 1e-5
     sszb = 0.2
-    kh = np.zeros((nlay, nrow, ncol), np.float)
-    kv = np.zeros((nlay, nrow, ncol), np.float)
-    ss = np.zeros((nlay, nrow, ncol), np.float)
-    ssz = np.zeros((nlay, nrow, ncol), np.float)
+    kh = np.zeros((nlay, nrow, ncol), float)
+    kv = np.zeros((nlay, nrow, ncol), float)
+    ss = np.zeros((nlay, nrow, ncol), float)
+    ssz = np.zeros((nlay, nrow, ncol), float)
     for k in range(0, nlay):
         for i in range(0, ncol):
             f = r[i] * 2.0 * math.pi
@@ -140,7 +140,7 @@ def run(silent=False):
             kv[k, 0, i] = kvb * f
             ss[k, 0, i] = ssb * f
             ssz[k, 0, i] = sszb * f
-    z = np.ones((nlay), np.float)
+    z = np.ones((nlay), float)
     z = -100.0 * z
 
     nwell = 1
@@ -153,7 +153,7 @@ def run(silent=False):
     wellQpm = wellQ / abs(wellbtm)
     well_data = {}
     for ip in range(0, nper):
-        welllist = np.zeros((nwell, 4), np.float)
+        welllist = np.zeros((nwell, 4), float)
         for iw in range(0, nwell):
             if ip == 0:
                 b = zall[iw] - zall[iw + 1]
@@ -168,7 +168,7 @@ def run(silent=False):
             welllist[iw, 3] = q
         well_data[ip] = welllist.copy()
 
-    ihead = np.zeros((nlay), np.float)
+    ihead = np.zeros((nlay), float)
 
     ocspd = {}
     for i in range(0, nper):
@@ -263,13 +263,13 @@ def run(silent=False):
     ibound = np.ones((nlay_swt, nrow, ncol), "int")
     for k in range(0, nlay_swt):
         ibound[k, 0, ncol - 1] = -1
-    bot = np.zeros((nlay_swt, nrow, ncol), np.float)
+    bot = np.zeros((nlay_swt, nrow, ncol), float)
     zall = [0, -20.0, -40.0, -60.0, -80.0, -100.0, -120.0]
     dz = 120.0 / nlay_swt
     tb = np.arange(nlay_swt) * -dz - dz
-    sconc = np.zeros((nlay_swt, nrow, ncol), np.float)
-    icbund = np.ones((nlay_swt, nrow, ncol), np.int)
-    strt = np.zeros((nlay_swt, nrow, ncol), np.float)
+    sconc = np.zeros((nlay_swt, nrow, ncol), float)
+    icbund = np.ones((nlay_swt, nrow, ncol), int)
+    strt = np.zeros((nlay_swt, nrow, ncol), float)
     pressure = 0.0
     g = 9.81
     z = -dz / 2.0  # cell center
@@ -295,10 +295,10 @@ def run(silent=False):
     kvb = (0.0000000000100 * 1000.0 * 9.81 / 0.001) * 60 * 60 * 24
     ssb = 1e-5
     sszb = 0.2
-    kh = np.zeros((nlay_swt, nrow, ncol), np.float)
-    kv = np.zeros((nlay_swt, nrow, ncol), np.float)
-    ss = np.zeros((nlay_swt, nrow, ncol), np.float)
-    ssz = np.zeros((nlay_swt, nrow, ncol), np.float)
+    kh = np.zeros((nlay_swt, nrow, ncol), float)
+    kv = np.zeros((nlay_swt, nrow, ncol), float)
+    ss = np.zeros((nlay_swt, nrow, ncol), float)
+    ssz = np.zeros((nlay_swt, nrow, ncol), float)
     for k in range(0, nlay_swt):
         for i in range(0, ncol):
             f = r[i] * 2.0 * math.pi
@@ -319,7 +319,7 @@ def run(silent=False):
     wellbtm = -20.0
     wellQpm = wellQ / abs(wellbtm)
     for ip in range(0, nper):
-        welllist = np.zeros((nwell, 4), np.float)
+        welllist = np.zeros((nwell, 4), float)
         ssmlist = []
         for iw in range(0, nwell):
             if ip == 0:
@@ -426,7 +426,7 @@ def run(silent=False):
     uobj = flopy.utils.UcnFile(ucnfile)
     times = uobj.get_times()
     print(times)
-    conc = np.zeros((len(times), nlay_swt, ncol), np.float)
+    conc = np.zeros((len(times), nlay_swt, ncol), float)
     for idx, tt in enumerate(times):
         c = uobj.get_data(totim=tt)
         for ilay in range(0, nlay_swt):
@@ -435,7 +435,7 @@ def run(silent=False):
 
     # spatial data
     # swi2
-    bot = np.zeros((1, ncol, nlay), np.float)
+    bot = np.zeros((1, ncol, nlay), float)
     dz = 100.0 / float(nlay - 1)
     zall = -np.arange(0, 100 + dz, dz)
     zall = np.append(zall, -120.0)
@@ -446,7 +446,7 @@ def run(silent=False):
             bot[0, i, k] = tb[k]
     # seawat
     swt_dz = 120.0 / nlay_swt
-    swt_tb = np.zeros((nlay_swt), np.float)
+    swt_tb = np.zeros((nlay_swt), float)
     zc = -swt_dz / 2.0
     for klay in range(0, nlay_swt):
         swt_tb[klay] = zc
@@ -529,8 +529,8 @@ def run(silent=False):
     leg._drawFrame = False
     # data items
     for itime in range(0, nswi_times):
-        zb = np.zeros((ncol), np.float)
-        zs = np.zeros((ncol), np.float)
+        zb = np.zeros((ncol), float)
+        zs = np.zeros((ncol), float)
         for icol in range(0, ncol):
             for klay in range(0, nlay):
                 # top and bottom of layer
diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py
index 8e5319396d..e2f2ba5310 100644
--- a/flopy/discretization/grid.py
+++ b/flopy/discretization/grid.py
@@ -30,13 +30,13 @@ class Grid(object):
     ----------
     grid_type : enumeration
         type of model grid ('structured', 'vertex', 'unstructured')
-    top : ndarray(np.float)
+    top : ndarray(float)
         top elevations of cells in topmost layer
-    botm : ndarray(np.float)
+    botm : ndarray(float)
         bottom elevations of all cells
-    idomain : ndarray(np.int)
+    idomain : ndarray(int)
         ibound/idomain value for each cell
-    lenuni : ndarray(np.int)
+    lenuni : ndarray(int)
         model length units
     origin_loc : str
         Corner of the model grid that is the model origin
@@ -54,11 +54,11 @@ class Grid(object):
     ----------
     grid_type : enumeration
         type of model grid ('structured', 'vertex', 'unstructured')
-    top : ndarray(np.float)
+    top : ndarray(float)
         top elevations of cells in topmost layer
-    botm : ndarray(np.float)
+    botm : ndarray(float)
         bottom elevations of all cells
-    idomain : ndarray(np.int)
+    idomain : ndarray(int)
         ibound/idomain value for each cell
     proj4 : proj4 SpatialReference
         spatial reference locates the grid in a coordinate system
diff --git a/flopy/discretization/unstructuredgrid.py b/flopy/discretization/unstructuredgrid.py
index 7560f7deab..9521d26fae 100644
--- a/flopy/discretization/unstructuredgrid.py
+++ b/flopy/discretization/unstructuredgrid.py
@@ -27,10 +27,10 @@ class UnstructuredGrid(Grid):
     ncpl : ndarray
         one dimensional array of size nlay with the number of cells in each
         layer.  This can also be passed in as a tuple or list as long as it
-        can be set using ncpl = np.array(ncpl, dtype=np.int).  The sum of ncpl
+        can be set using ncpl = np.array(ncpl, dtype=int).  The sum of ncpl
         must be equal to the number of cells in the grid.  ncpl is optional
         and if it is not passed in, then it is is set using
-        ncpl = np.array([len(iverts)], dtype=np.int), which means that all
+        ncpl = np.array([len(iverts)], dtype=int), which means that all
         cells in the grid are contained in a single plottable layer.
         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
@@ -140,9 +140,9 @@ def __init__(
 
     def set_ncpl(self, ncpl):
         if isinstance(ncpl, int):
-            ncpl = np.array([ncpl], dtype=np.int)
+            ncpl = np.array([ncpl], dtype=int)
         if isinstance(ncpl, (list, tuple, np.ndarray)):
-            ncpl = np.array(ncpl, dtype=np.int)
+            ncpl = np.array(ncpl, dtype=int)
         else:
             raise TypeError("ncpl must be a list, tuple or ndarray")
         assert ncpl.ndim == 1, "ncpl must be 1d"
@@ -499,9 +499,9 @@ def from_argus_export(cls, fname, nlay=1):
         ncells, nverts = ll[0:2]
         ncells = int(ncells)
         nverts = int(nverts)
-        verts = np.empty((nverts, 3), dtype=np.float)
-        xc = np.empty((ncells), dtype=np.float)
-        yc = np.empty((ncells), dtype=np.float)
+        verts = np.empty((nverts, 3), dtype=float)
+        xc = np.empty((ncells), dtype=float)
+        yc = np.empty((ncells), dtype=float)
 
         # read the vertices
         f.readline()
diff --git a/flopy/export/shapefile_utils.py b/flopy/export/shapefile_utils.py
index 36b5410192..6dcfe52d85 100755
--- a/flopy/export/shapefile_utils.py
+++ b/flopy/export/shapefile_utils.py
@@ -192,7 +192,7 @@ def write_grid_shapefile(
         names = enforce_10ch_limit(names)
 
     # flag nan values and explicitly set the dtypes
-    if at.dtype in [np.float, np.float32, np.float64]:
+    if at.dtype in [float, np.float32, np.float64]:
         at[np.isnan(at)] = nan_val
     at = np.array([tuple(i) for i in at], dtype=dtypes)
 
@@ -490,12 +490,12 @@ def get_pyshp_field_info(dtypename):
 def get_pyshp_field_dtypes(code):
     """Returns a numpy dtype for a pyshp field type."""
     dtypes = {
-        "N": np.int,
-        "F": np.float,
-        "L": np.bool,
-        "C": np.object,
+        "N": int,
+        "F": float,
+        "L": bool,
+        "C": object,
     }
-    return dtypes.get(code, np.object)
+    return dtypes.get(code, object)
 
 
 def shp2recarray(shpname):
@@ -524,7 +524,7 @@ def shp2recarray(shpname):
     records = [
         tuple(r) + (geoms[i],) for i, r in enumerate(sfobj.iterRecords())
     ]
-    dtype += [("geometry", np.object)]
+    dtype += [("geometry", object)]
 
     recarray = np.array(records, dtype=dtype).view(np.recarray)
     return recarray
diff --git a/flopy/export/utils.py b/flopy/export/utils.py
index a9d0dd0c0f..e3f5793939 100644
--- a/flopy/export/utils.py
+++ b/flopy/export/utils.py
@@ -20,7 +20,7 @@
 NC_PRECISION_TYPE = {
     np.float64: "f8",
     np.float32: "f4",
-    np.int: "i4",
+    int: "i4",
     np.int64: "i4",
     np.int32: "i4",
 }
@@ -1040,7 +1040,7 @@ def transient2d_export(f, t2d, fmt=None, **kwargs):
         array = t2d.array
         # f.log("getting 4D array for {0}".format(t2d.name_base))
         with np.errstate(invalid="ignore"):
-            if array.dtype not in [int, np.int, np.int32, np.int64]:
+            if array.dtype not in [int, np.int32, np.int64]:
                 if mask is not None:
                     array[:, 0, mask] = np.NaN
                 array[array <= min_valid] = np.NaN
@@ -1216,7 +1216,7 @@ def array3d_export(f, u3d, fmt=None, **kwargs):
         # runtime warning issued in some cases - need to track down cause
         # happens when NaN is already in array
         with np.errstate(invalid="ignore"):
-            if array.dtype not in [int, np.int, np.int32, np.int64]:
+            if array.dtype not in [int, np.int32, np.int64]:
                 # if u3d.model.modelgrid.bas6 is not None and "ibound" not
                 # in var_name:
                 #    array[u3d.model.modelgrid.bas6.ibound.array == 0] =
@@ -1359,7 +1359,7 @@ def array2d_export(f, u2d, fmt=None, **kwargs):
         # f.log("getting 2D array for {0}".format(u2d.name))
 
         with np.errstate(invalid="ignore"):
-            if array.dtype not in [int, np.int, np.int32, np.int64]:
+            if array.dtype not in [int, np.int32, np.int64]:
                 if (
                     modelgrid.idomain is not None
                     and "ibound" not in u2d.name.lower()
diff --git a/flopy/export/vtk.py b/flopy/export/vtk.py
index d6dc2cf375..6328da2d10 100644
--- a/flopy/export/vtk.py
+++ b/flopy/export/vtk.py
@@ -1070,12 +1070,12 @@ def _configure_data_arrays(self):
         """
         # build 1d index array
         shape1d = self.shape[0] * self.shape[1] * self.shape[2]
-        actwcells1d = np.zeros(shape1d, dtype=np.int)
+        actwcells1d = np.zeros(shape1d, dtype=int)
         if self.has_point_data:
             shape1d_verts = (
                 self.shape_verts[0] * self.shape_verts[1] * self.shape_verts[2]
             )
-            actwcells1d_verts = np.zeros(shape1d_verts, dtype=np.int)
+            actwcells1d_verts = np.zeros(shape1d_verts, dtype=int)
 
         # loop through arrays
         for a in self.arrays.values():
diff --git a/flopy/mf6/coordinates/simulationtime.py b/flopy/mf6/coordinates/simulationtime.py
index 45e0b172b2..a3dc6b15ce 100644
--- a/flopy/mf6/coordinates/simulationtime.py
+++ b/flopy/mf6/coordinates/simulationtime.py
@@ -98,7 +98,7 @@ def __init__(self, simdata):
     #        self.calendar_start_time = calendar_start_time
 
     #    def get_stress_period_array(self):
-    #        return np.arange(1, self.get_num_stress_periods(), 1, np.int)
+    #        return np.arange(1, self.get_num_stress_periods(), 1, int)
 
     def get_time_units(self):
         time_units = self.simdata.mfdata[
diff --git a/flopy/mf6/data/mffileaccess.py b/flopy/mf6/data/mffileaccess.py
index 0499d5cd7f..336abc0b61 100644
--- a/flopy/mf6/data/mffileaccess.py
+++ b/flopy/mf6/data/mffileaccess.py
@@ -190,7 +190,7 @@ def datum_to_numpy_type(datum_type):
         elif datum_type == DatumType.double_precision:
             return np.float64, "double"
         elif datum_type == DatumType.string or datum_type == DatumType.keyword:
-            return np.str, "str"
+            return str, "str"
         else:
             return None, None
 
diff --git a/flopy/mf6/data/mfstructure.py b/flopy/mf6/data/mfstructure.py
index 2428302b0f..22d46aeb51 100644
--- a/flopy/mf6/data/mfstructure.py
+++ b/flopy/mf6/data/mfstructure.py
@@ -1862,7 +1862,7 @@ def get_datum_type(self, numpy_type=False, return_enum_type=False):
                         elif var_type[0] == DatumType.integer:
                             return np.int32
                         else:
-                            return np.object
+                            return object
                     else:
                         return var_type[2]
         return None
diff --git a/flopy/mf6/mfmodel.py b/flopy/mf6/mfmodel.py
index 1340cd0b44..18cf5678fb 100644
--- a/flopy/mf6/mfmodel.py
+++ b/flopy/mf6/mfmodel.py
@@ -321,7 +321,7 @@ def modeldiscrit(self):
         elif self.get_grid_type() == DiscretizationType.DISU:
             dis = self.get_package("disu")
             nodes = dis.nodes.get_data()
-            ncpl = np.array([nodes], dtype=np.int)
+            ncpl = np.array([nodes], dtype=int)
             return UnstructuredGrid(ncpl=ncpl)
 
     @property
@@ -410,7 +410,7 @@ def modelgrid(self):
             iac = dis.iac.array
             ncpl = UnstructuredGrid.ncpl_from_ihc(ihc, iac)
             if ncpl is None:
-                ncpl = np.array([dis.nodes.get_data()], dtype=np.int)
+                ncpl = np.array([dis.nodes.get_data()], dtype=int)
             cell2d = dis.cell2d.array
             idomain = np.ones(dis.nodes.array, np.int32)
             if cell2d is None:
diff --git a/flopy/mf6/utils/mfobservation.py b/flopy/mf6/utils/mfobservation.py
index ecec8b3e0c..04a1c6c6f4 100644
--- a/flopy/mf6/utils/mfobservation.py
+++ b/flopy/mf6/utils/mfobservation.py
@@ -373,11 +373,11 @@ def get_obs_data(self, key=None, idx=None, totim=None):
         """
         data = self.get_data(key=key, idx=idx, totim=totim)
         # create x-array coordinates from time and header
-        totim = data.T[0][1:].astype(np.float)
-        header = data[0][1:].astype(np.str)
+        totim = data.T[0][1:].astype(float)
+        header = data[0][1:].astype(str)
 
         # strip time and header off of data
-        data = data[1:, 1:].astype(np.float)
+        data = data[1:, 1:].astype(float)
 
         return data
 
diff --git a/flopy/modflow/mfag.py b/flopy/modflow/mfag.py
index bef11307d4..3478ad9304 100644
--- a/flopy/modflow/mfag.py
+++ b/flopy/modflow/mfag.py
@@ -618,64 +618,64 @@ def get_default_dtype(maxells=0, block="well"):
         """
         if block == "well":
             dtype = [
-                ("k", np.int),
-                ("i", np.int),
-                ("j", np.int),
-                ("flux", np.float),
+                ("k", int),
+                ("i", int),
+                ("j", int),
+                ("flux", float),
             ]
 
         elif block == "tabfile_well":
             dtype = [
-                ("unit", np.int),
-                ("tabval", np.int),
-                ("k", np.int),
-                ("i", np.int),
-                ("j", np.int),
+                ("unit", int),
+                ("tabval", int),
+                ("k", int),
+                ("i", int),
+                ("j", int),
             ]
 
         elif block == "time series":
-            dtype = [("keyword", np.object), ("id", np.int), ("unit", np.int)]
+            dtype = [("keyword", object), ("id", int), ("unit", int)]
 
         elif block == "irrdiversion":
             dtype = [
-                ("segid", np.int),
-                ("numcell", np.int),
-                ("period", np.float),
-                ("triggerfact", np.float),
+                ("segid", int),
+                ("numcell", int),
+                ("period", float),
+                ("triggerfact", float),
             ]
 
             for i in range(maxells):
                 dtype += [
-                    ("i{}".format(i), np.int),
-                    ("j{}".format(i), np.int),
-                    ("eff_fact{}".format(i), np.float),
-                    ("field_fact{}".format(i), np.float),
+                    ("i{}".format(i), int),
+                    ("j{}".format(i), int),
+                    ("eff_fact{}".format(i), float),
+                    ("field_fact{}".format(i), float),
                 ]
 
         elif block == "irrwell":
             dtype = [
-                ("wellid", np.int),
-                ("numcell", np.int),
-                ("period", np.float),
-                ("triggerfact", np.float),
+                ("wellid", int),
+                ("numcell", int),
+                ("period", float),
+                ("triggerfact", float),
             ]
 
             for i in range(maxells):
                 dtype += [
-                    ("i{}".format(i), np.int),
-                    ("j{}".format(i), np.int),
-                    ("eff_fact{}".format(i), np.float),
-                    ("field_fact{}".format(i), np.float),
+                    ("i{}".format(i), int),
+                    ("j{}".format(i), int),
+                    ("eff_fact{}".format(i), float),
+                    ("field_fact{}".format(i), float),
                 ]
 
         elif block == "supwell":
-            dtype = [("wellid", np.int), ("numcell", np.int)]
+            dtype = [("wellid", int), ("numcell", int)]
 
             for i in range(maxells):
                 dtype += [
-                    ("segid{}".format(i), np.int),
-                    ("fracsup{}".format(i), np.float),
-                    ("fracsupmax{}".format(i), np.float),
+                    ("segid{}".format(i), int),
+                    ("fracsup{}".format(i), float),
+                    ("fracsupmax{}".format(i), float),
                 ]
 
         else:
diff --git a/flopy/modflow/mfchd.py b/flopy/modflow/mfchd.py
index 3797e1a84d..ae9b1b07ac 100644
--- a/flopy/modflow/mfchd.py
+++ b/flopy/modflow/mfchd.py
@@ -212,9 +212,9 @@ def get_default_dtype(structured=True):
         if structured:
             dtype = np.dtype(
                 [
-                    ("k", np.int),
-                    ("i", np.int),
-                    ("j", np.int),
+                    ("k", int),
+                    ("i", int),
+                    ("j", int),
                     ("shead", np.float32),
                     ("ehead", np.float32),
                 ]
@@ -222,7 +222,7 @@ def get_default_dtype(structured=True):
         else:
             dtype = np.dtype(
                 [
-                    ("node", np.int),
+                    ("node", int),
                     ("shead", np.float32),
                     ("ehead", np.float32),
                 ]
diff --git a/flopy/modflow/mfdis.py b/flopy/modflow/mfdis.py
index 83f98a6fcc..5dc7c706ec 100644
--- a/flopy/modflow/mfdis.py
+++ b/flopy/modflow/mfdis.py
@@ -62,7 +62,7 @@ class ModflowDis(Package):
         Number of time steps in each stress period (default is 1).
     tsmult : float or array of floats (nper)
         Time step multiplier (default is 1.0).
-    steady : boolean or array of boolean (nper)
+    steady : bool or array of bool (nper)
         true or False indicating whether or not stress period is steady state
         (default is True).
     itmuni : int
@@ -236,9 +236,7 @@ def __init__(
         self.tsmult = Util2d(
             model, (self.nper,), np.float32, tsmult, name="tsmult"
         )
-        self.steady = Util2d(
-            model, (self.nper,), np.bool, steady, name="steady"
-        )
+        self.steady = Util2d(model, (self.nper,), bool, steady, name="steady")
 
         try:
             self.itmuni = int(itmuni)
@@ -356,7 +354,7 @@ def get_totim(self):
                 totim.append(t)
                 if m > 1:
                     dt *= m
-        return np.array(totim, dtype=np.float)
+        return np.array(totim, dtype=float)
 
     def get_final_totim(self):
         """
@@ -674,7 +672,7 @@ def write_file(self, check=True):
 
         Parameters
         ----------
-        check : boolean
+        check : bool
             Check package data for common errors. (default True)
 
         Returns
@@ -823,7 +821,7 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
             handle.  In this case ext_unit_dict is required, which can be
             constructed using the function
             :class:`flopy.utils.mfreadnam.parsenamefile`.
-        check : boolean
+        check : bool
             Check package data for common errors. (default True)
 
         Returns
@@ -937,13 +935,13 @@ def load(cls, f, model, ext_unit_dict=None, check=True):
                 )
             )
             print("   loading laycbd...")
-        laycbd = np.zeros(nlay, dtype=np.int)
+        laycbd = np.zeros(nlay, dtype=int)
         d = 0
         while True:
             line = f.readline()
             raw = line.strip("\n").split()
             for val in raw:
-                if (np.int(val)) != 0:
+                if int(val) != 0:
                     laycbd[d] = 1
                 d += 1
                 if d == nlay:
diff --git a/flopy/modflow/mfdisu.py b/flopy/modflow/mfdisu.py
index 27ca90a84a..d960e45c1c 100644
--- a/flopy/modflow/mfdisu.py
+++ b/flopy/modflow/mfdisu.py
@@ -157,8 +157,8 @@ class ModflowDisU(Package):
         Number of time steps in each stress period (default is 1).
     tsmult : float or array of floats (nper)
         Time step multiplier (default is 1.0).
-    steady : boolean or array of boolean (nper)
-        true or False indicating whether or not stress period is steady state
+    steady : bool or array of bool (nper)
+        True or False indicating whether or not stress period is steady state
         (default is True).
     extension : string
         Filename extension (default is 'dis')
@@ -442,9 +442,7 @@ def __init__(
         self.tsmult = Util2d(
             model, (self.nper,), np.float32, tsmult, name="tsmult"
         )
-        self.steady = Util2d(
-            model, (self.nper,), np.bool, steady, name="steady"
-        )
+        self.steady = Util2d(model, (self.nper,), bool, steady, name="steady")
 
         self.itmuni_dict = {
             0: "undefined",
@@ -550,7 +548,7 @@ def load(cls, f, model, ext_unit_dict=None, check=False):
             handle.  In this case ext_unit_dict is required, which can be
             constructed using the function
             :class:`flopy.utils.mfreadnam.parsenamefile`.
-        check : boolean
+        check : bool
             Check package data for common errors. (default False; not setup yet)
 
         Returns
diff --git a/flopy/modflow/mfdrn.py b/flopy/modflow/mfdrn.py
index 26dd2df390..3de298664b 100644
--- a/flopy/modflow/mfdrn.py
+++ b/flopy/modflow/mfdrn.py
@@ -197,9 +197,9 @@ def get_default_dtype(structured=True, is_drt=False):
             if not is_drt:
                 dtype = np.dtype(
                     [
-                        ("k", np.int),
-                        ("i", np.int),
-                        ("j", np.int),
+                        ("k", int),
+                        ("i", int),
+                        ("j", int),
                         ("elev", np.float32),
                         ("cond", np.float32),
                     ]
@@ -207,20 +207,20 @@ def get_default_dtype(structured=True, is_drt=False):
             else:
                 dtype = np.dtype(
                     [
-                        ("k", np.int),
-                        ("i", np.int),
-                        ("j", np.int),
+                        ("k", int),
+                        ("i", int),
+                        ("j", int),
                         ("elev", np.float32),
                         ("cond", np.float32),
-                        ("layr", np.int),
-                        ("rowr", np.int),
-                        ("colr", np.int),
+                        ("layr", int),
+                        ("rowr", int),
+                        ("colr", int),
                         ("rfprop", np.float32),
                     ]
                 )
         else:
             dtype = np.dtype(
-                [("node", np.int), ("elev", np.float32), ("cond", np.float32)]
+                [("node", int), ("elev", np.float32), ("cond", np.float32)]
             )
         return dtype
 
diff --git a/flopy/modflow/mfdrt.py b/flopy/modflow/mfdrt.py
index 8664a07a37..294ccb637c 100644
--- a/flopy/modflow/mfdrt.py
+++ b/flopy/modflow/mfdrt.py
@@ -194,26 +194,26 @@ def get_default_dtype(structured=True):
         if structured:
             dtype = np.dtype(
                 [
-                    ("k", np.int),
-                    ("i", np.int),
-                    ("j", np.int),
+                    ("k", int),
+                    ("i", int),
+                    ("j", int),
                     ("elev", np.float32),
                     ("cond", np.float32),
-                    ("layr", np.int),
-                    ("rowr", np.int),
-                    ("colr", np.int),
+                    ("layr", int),
+                    ("rowr", int),
+                    ("colr", int),
                     ("rfprop", np.float32),
                 ]
             )
         else:
             dtype = np.dtype(
                 [
-                    ("inode", np.int),
+                    ("inode", int),
                     ("elev", np.float32),
                     ("cond", np.float32),
-                    ("layr", np.int),
-                    ("rowr", np.int),
-                    ("colr", np.int),
+                    ("layr", int),
+                    ("rowr", int),
+                    ("colr", int),
                     ("rfprop", np.float32),
                 ]
             )
diff --git a/flopy/modflow/mffhb.py b/flopy/modflow/mffhb.py
index 2c02cc9f0a..49d979f30f 100644
--- a/flopy/modflow/mffhb.py
+++ b/flopy/modflow/mffhb.py
@@ -325,10 +325,10 @@ def get_empty(ncells=0, nbdtim=1, structured=True, head=False):
     @staticmethod
     def get_default_dtype(nbdtim=1, structured=True, head=False):
         if structured:
-            dtype = [("k", np.int), ("i", np.int), ("j", np.int)]
+            dtype = [("k", int), ("i", int), ("j", int)]
         else:
-            dtype = [("node", np.int)]
-        dtype.append(("iaux", np.int))
+            dtype = [("node", int)]
+        dtype.append(("iaux", int))
         for n in range(nbdtim):
             if head:
                 name = "sbhed{}".format(n + 1)
diff --git a/flopy/modflow/mfgage.py b/flopy/modflow/mfgage.py
index 0059e4a412..36e45250e6 100644
--- a/flopy/modflow/mfgage.py
+++ b/flopy/modflow/mfgage.py
@@ -230,10 +230,10 @@ def __init__(
     def get_default_dtype():
         dtype = np.dtype(
             [
-                ("gageloc", np.int),
-                ("gagerch", np.int),
-                ("unit", np.int),
-                ("outtype", np.int),
+                ("gageloc", int),
+                ("gagerch", int),
+                ("unit", int),
+                ("outtype", int),
             ]
         )
         return dtype
diff --git a/flopy/modflow/mfghb.py b/flopy/modflow/mfghb.py
index ea41d98931..e6bed6d656 100644
--- a/flopy/modflow/mfghb.py
+++ b/flopy/modflow/mfghb.py
@@ -251,16 +251,16 @@ def get_default_dtype(structured=True):
         if structured:
             dtype = np.dtype(
                 [
-                    ("k", np.int),
-                    ("i", np.int),
-                    ("j", np.int),
+                    ("k", int),
+                    ("i", int),
+                    ("j", int),
                     ("bhead", np.float32),
                     ("cond", np.float32),
                 ]
             )
         else:
             dtype = np.dtype(
-                [("node", np.int), ("bhead", np.float32), ("cond", np.float32)]
+                [("node", int), ("bhead", np.float32), ("cond", np.float32)]
             )
         return dtype
 
diff --git a/flopy/modflow/mfhfb.py b/flopy/modflow/mfhfb.py
index f6adb3099f..30ad445358 100644
--- a/flopy/modflow/mfhfb.py
+++ b/flopy/modflow/mfhfb.py
@@ -241,11 +241,11 @@ def get_default_dtype(structured=True):
         if structured:
             dtype = np.dtype(
                 [
-                    ("k", np.int),
-                    ("irow1", np.int),
-                    ("icol1", np.int),
-                    ("irow2", np.int),
-                    ("icol2", np.int),
+                    ("k", int),
+                    ("irow1", int),
+                    ("icol1", int),
+                    ("irow2", int),
+                    ("icol2", int),
                     ("hydchr", np.float32),
                 ]
             )
@@ -372,12 +372,12 @@ def load(cls, f, model, ext_unit_dict=None):
 
                 #
                 if model.mfpar.pval is None:
-                    parval = np.float(par_dict["parval"])
+                    parval = float(par_dict["parval"])
                 else:
                     try:
-                        parval = np.float(model.mfpar.pval.pval_dict[pname])
+                        parval = float(model.mfpar.pval.pval_dict[pname])
                     except:
-                        parval = np.float(par_dict["parval"])
+                        parval = float(par_dict["parval"])
 
                 # fill current parameter data (par_current)
                 for ibnd, t in enumerate(data_dict):
diff --git a/flopy/modflow/mfhob.py b/flopy/modflow/mfhob.py
index 12c6ae5e30..d5b71f822c 100755
--- a/flopy/modflow/mfhob.py
+++ b/flopy/modflow/mfhob.py
@@ -637,7 +637,7 @@ def __init__(
 
         # convert passed time_series_data to a numpy array
         if isinstance(time_series_data, list):
-            time_series_data = np.array(time_series_data, dtype=np.float)
+            time_series_data = np.array(time_series_data, dtype=float)
 
         # if a single observation is passed as a list reshape to a
         # two-dimensional numpy array
@@ -730,7 +730,7 @@ def _get_dtype(self):
         dtype = np.dtype(
             [
                 ("totim", np.float32),
-                ("irefsp", np.int),
+                ("irefsp", int),
                 ("toffset", np.float32),
                 ("hobs", np.float32),
                 ("obsname", "|S12"),
diff --git a/flopy/modflow/mfhyd.py b/flopy/modflow/mfhyd.py
index d3341d81cc..53792707b0 100644
--- a/flopy/modflow/mfhyd.py
+++ b/flopy/modflow/mfhyd.py
@@ -202,7 +202,7 @@ def __init__(
                 obs["hydlbl"][idx] = obsdata[idx][6]
             obsdata = obs
         elif isinstance(obsdata, np.ndarray):
-            if obsdata.dtype == np.object:
+            if obsdata.dtype == object:
                 if obsdata.shape[1] != len(dtype):
                     raise IndexError("Incorrect number of fields for obsdata")
                 obsdata = obsdata.transpose()
@@ -278,7 +278,7 @@ def get_default_dtype():
                 ("pckg", "|S3"),
                 ("arr", "|S2"),
                 ("intyp", "|S1"),
-                ("klay", np.int),
+                ("klay", int),
                 ("xl", np.float32),
                 ("yl", np.float32),
                 ("hydlbl", "|S14"),
diff --git a/flopy/modflow/mflak.py b/flopy/modflow/mflak.py
index d6b5002097..acfeab3066 100644
--- a/flopy/modflow/mflak.py
+++ b/flopy/modflow/mflak.py
@@ -385,7 +385,7 @@ def __init__(
 
         if isinstance(stages, float):
             if self.nlakes == 1:
-                stages = np.array([self.nlakes], dtype=np.float) * stages
+                stages = np.array([self.nlakes], dtype=float) * stages
             else:
                 stages = np.ones(self.nlakes, dtype=float) * stages
         elif isinstance(stages, list):
@@ -397,7 +397,7 @@ def __init__(
             raise Exception(err)
         self.stages = stages
         if stage_range is None:
-            stage_range = np.ones((nlakes, 2), dtype=np.float)
+            stage_range = np.ones((nlakes, 2), dtype=float)
             stage_range[:, 0] = -10000.0
             stage_range[:, 1] = 10000.0
         else:
@@ -466,7 +466,7 @@ def __init__(
                 elif isinstance(value, float) or isinstance(value, int):
                     td = {}
                     for k in range(self.nlakes):
-                        td[k] = (np.ones(6, dtype=np.float) * value).tolist()
+                        td[k] = (np.ones(6, dtype=float) * value).tolist()
                     flux_data[key] = td
                 elif isinstance(value, dict):
                     try:
diff --git a/flopy/modflow/mfmnw1.py b/flopy/modflow/mfmnw1.py
index 01788a5a62..0e2d324518 100644
--- a/flopy/modflow/mfmnw1.py
+++ b/flopy/modflow/mfmnw1.py
@@ -186,24 +186,24 @@ def get_default_dtype(structured=True):
         if structured:
             return np.dtype(
                 [
-                    ("mnw_no", np.int),
-                    ("k", np.int),
-                    ("i", np.int),
-                    ("j", np.int),
+                    ("mnw_no", int),
+                    ("k", int),
+                    ("i", int),
+                    ("j", int),
                     ("qdes", np.float32),
-                    ("mntxt", np.object),
+                    ("mntxt", object),
                     ("qwval", np.float32),
                     ("rw", np.float32),
                     ("skin", np.float32),
                     ("hlim", np.float32),
                     ("href", np.float32),
-                    ("dd", np.object),
-                    ("iqwgrp", np.object),
-                    ("cpc", np.object),
-                    ("qcut", np.object),
+                    ("dd", object),
+                    ("iqwgrp", object),
+                    ("cpc", object),
+                    ("qcut", object),
                     ("qfrcmn", np.float32),
                     ("qfrcmx", np.float32),
-                    ("label", np.object),
+                    ("label", object),
                 ]
             )
         else:
diff --git a/flopy/modflow/mfmnw2.py b/flopy/modflow/mfmnw2.py
index a9023a201b..3a74df17b4 100644
--- a/flopy/modflow/mfmnw2.py
+++ b/flopy/modflow/mfmnw2.py
@@ -545,15 +545,15 @@ def get_default_spd_dtype(structured=True):
         if structured:
             return np.dtype(
                 [
-                    ("k", np.int),
-                    ("i", np.int),
-                    ("j", np.int),
-                    ("per", np.int),
+                    ("k", int),
+                    ("i", int),
+                    ("j", int),
+                    ("per", int),
                     ("qdes", np.float32),
-                    ("capmult", np.int),
+                    ("capmult", int),
                     ("cprime", np.float32),
                     ("hlim", np.float32),
-                    ("qcut", np.int),
+                    ("qcut", int),
                     ("qfrcmn", np.float32),
                     ("qfrcmx", np.float32),
                 ]
@@ -1218,17 +1218,17 @@ def get_default_node_dtype(structured=True):
         if structured:
             return np.dtype(
                 [
-                    ("k", np.int),
-                    ("i", np.int),
-                    ("j", np.int),
+                    ("k", int),
+                    ("i", int),
+                    ("j", int),
                     ("ztop", np.float32),
                     ("zbotm", np.float32),
-                    ("wellid", np.object),
-                    ("losstype", np.object),
-                    ("pumploc", np.int),
-                    ("qlimit", np.int),
-                    ("ppflag", np.int),
-                    ("pumpcap", np.int),
+                    ("wellid", object),
+                    ("losstype", object),
+                    ("pumploc", int),
+                    ("qlimit", int),
+                    ("ppflag", int),
+                    ("pumpcap", int),
                     ("rw", np.float32),
                     ("rskin", np.float32),
                     ("kskin", np.float32),
@@ -1237,12 +1237,12 @@ def get_default_node_dtype(structured=True):
                     ("P", np.float32),
                     ("cwc", np.float32),
                     ("pp", np.float32),
-                    ("pumplay", np.int),
-                    ("pumprow", np.int),
-                    ("pumpcol", np.int),
+                    ("pumplay", int),
+                    ("pumprow", int),
+                    ("pumpcol", int),
                     ("zpump", np.float32),
                     ("hlim", np.float32),
-                    ("qcut", np.int),
+                    ("qcut", int),
                     ("qfrcmn", np.float32),
                     ("qfrcmx", np.float32),
                     ("hlift", np.float32),
@@ -1307,15 +1307,15 @@ def get_default_spd_dtype(structured=True):
         if structured:
             return np.dtype(
                 [
-                    ("k", np.int),
-                    ("i", np.int),
-                    ("j", np.int),
-                    ("wellid", np.object),
+                    ("k", int),
+                    ("i", int),
+                    ("j", int),
+                    ("wellid", object),
                     ("qdes", np.float32),
-                    ("capmult", np.int),
+                    ("capmult", int),
                     ("cprime", np.float32),
                     ("hlim", np.float32),
-                    ("qcut", np.int),
+                    ("qcut", int),
                     ("qfrcmn", np.float32),
                     ("qfrcmx", np.float32),
                 ]
diff --git a/flopy/modflow/mfpar.py b/flopy/modflow/mfpar.py
index cc5a1ce58d..ee4e5f813d 100644
--- a/flopy/modflow/mfpar.py
+++ b/flopy/modflow/mfpar.py
@@ -233,13 +233,13 @@ def load(f, npar, verbose=False):
                 partyp = t[1].lower()
                 if partyp not in par_types:
                     par_types.append(partyp)
-                parval = np.float(t[2])
-                nclu = np.int(t[3])
+                parval = float(t[2])
+                nclu = int(t[3])
                 clusters = []
                 for nc in range(nclu):
                     line = f.readline()
                     t = line.strip().split()
-                    lay = np.int(t[0])
+                    lay = int(t[0])
                     s = t[1]
                     if len(s) > 10:
                         s = s[0:10]
@@ -251,7 +251,7 @@ def load(f, npar, verbose=False):
                     iarr = []
                     for iv in t[3:]:
                         try:
-                            iz = int(np.int(iv))
+                            iz = int(iv)
                             if iz != 0:
                                 iarr.append(iz)
                         except:
@@ -313,12 +313,12 @@ def parameter_fill(model, shape, findkey, parm_dict, findlayer=None):
             partyp, parval = tdict["partyp"], tdict["parval"]
             nclu, clusters = tdict["nclu"], tdict["clusters"]
             if model.mfpar.pval is None:
-                pv = np.float(parval)
+                pv = float(parval)
             else:
                 try:
-                    pv = np.float(model.mfpar.pval.pval_dict[key.lower()])
+                    pv = float(model.mfpar.pval.pval_dict[key.lower()])
                 except:
-                    pv = np.float(parval)
+                    pv = float(parval)
             # print partyp, parval, nclu, clusters
             if partyp == findkey:
                 for [layer, mltarr, zonarr, izones] in clusters:
diff --git a/flopy/modflow/mfparbc.py b/flopy/modflow/mfparbc.py
index a12d9dd87a..68be7e1256 100644
--- a/flopy/modflow/mfparbc.py
+++ b/flopy/modflow/mfparbc.py
@@ -85,12 +85,12 @@ def load(cls, f, npar, dt, model, ext_unit_dict=None, verbose=False):
                     print('   loading parameter "{}"...'.format(parnam))
                 partyp = t[1].lower()
                 parval = t[2]
-                nlst = np.int(t[3])
+                nlst = int(t[3])
                 numinst = 1
                 timeVarying = False
                 if len(t) > 4:
                     if "instances" in t[4].lower():
-                        numinst = np.int(t[5])
+                        numinst = int(t[5])
                         timeVarying = True
                 pinst = {}
                 for inst in range(numinst):
@@ -156,12 +156,12 @@ def loadarray(f, npar, verbose=False):
                     print('   loading parameter "{}"...'.format(parnam))
                 partyp = t[1].lower()
                 parval = t[2]
-                nclu = np.int(t[3])
+                nclu = int(t[3])
                 numinst = 1
                 timeVarying = False
                 if len(t) > 4:
                     if "instances" in t[4].lower():
-                        numinst = np.int(t[5])
+                        numinst = int(t[5])
                         timeVarying = True
                 pinst = {}
                 for inst in range(numinst):
@@ -248,12 +248,12 @@ def parameter_bcfill(model, shape, parm_dict, pak_parms):
             pdict, idict = pak_parms.bc_parms[key]
             inst_data = idict[value]
             if model.mfpar.pval is None:
-                pv = np.float(pdict["parval"])
+                pv = float(pdict["parval"])
             else:
                 try:
-                    pv = np.float(model.mfpar.pval.pval_dict[key.lower()])
+                    pv = float(model.mfpar.pval.pval_dict[key.lower()])
                 except:
-                    pv = np.float(pdict["parval"])
+                    pv = float(pdict["parval"])
             for [mltarr, zonarr, izones] in inst_data:
                 model.parameter_load = True
                 # print mltarr, zonarr, izones
diff --git a/flopy/modflow/mfrch.py b/flopy/modflow/mfrch.py
index 8c1ed33849..063347b789 100644
--- a/flopy/modflow/mfrch.py
+++ b/flopy/modflow/mfrch.py
@@ -426,7 +426,7 @@ def load(cls, f, model, nper=None, ext_unit_dict=None, check=True):
         npar = 0
         if "parameter" in line.lower():
             raw = line.strip().split()
-            npar = np.int(raw[1])
+            npar = int(raw[1])
             if npar > 0:
                 if model.verbose:
                     txt = (
diff --git a/flopy/modflow/mfriv.py b/flopy/modflow/mfriv.py
index 61afd9d62a..b5d625fcbe 100644
--- a/flopy/modflow/mfriv.py
+++ b/flopy/modflow/mfriv.py
@@ -277,9 +277,9 @@ def get_default_dtype(structured=True):
         if structured:
             dtype = np.dtype(
                 [
-                    ("k", np.int),
-                    ("i", np.int),
-                    ("j", np.int),
+                    ("k", int),
+                    ("i", int),
+                    ("j", int),
                     ("stage", np.float32),
                     ("cond", np.float32),
                     ("rbot", np.float32),
@@ -288,7 +288,7 @@ def get_default_dtype(structured=True):
         else:
             dtype = np.dtype(
                 [
-                    ("node", np.int),
+                    ("node", int),
                     ("stage", np.float32),
                     ("cond", np.float32),
                     ("rbot", np.float32),
diff --git a/flopy/modflow/mfsfr2.py b/flopy/modflow/mfsfr2.py
index ba46411693..e26eab7669 100644
--- a/flopy/modflow/mfsfr2.py
+++ b/flopy/modflow/mfsfr2.py
@@ -754,12 +754,12 @@ def get_default_reach_dtype(structured=True):
             # include node column for structured grids (useful for indexing)
             return np.dtype(
                 [
-                    ("node", np.int),
-                    ("k", np.int),
-                    ("i", np.int),
-                    ("j", np.int),
-                    ("iseg", np.int),
-                    ("ireach", np.int),
+                    ("node", int),
+                    ("k", int),
+                    ("i", int),
+                    ("j", int),
+                    ("iseg", int),
+                    ("ireach", int),
                     ("rchlen", np.float32),
                     ("strtop", np.float32),
                     ("slope", np.float32),
@@ -769,16 +769,16 @@ def get_default_reach_dtype(structured=True):
                     ("thti", np.float32),
                     ("eps", np.float32),
                     ("uhc", np.float32),
-                    ("reachID", np.int),
-                    ("outreach", np.int),
+                    ("reachID", int),
+                    ("outreach", int),
                 ]
             )
         else:
             return np.dtype(
                 [
-                    ("node", np.int),
-                    ("iseg", np.int),
-                    ("ireach", np.int),
+                    ("node", int),
+                    ("iseg", int),
+                    ("ireach", int),
                     ("rchlen", np.float32),
                     ("strtop", np.float32),
                     ("slope", np.float32),
@@ -788,8 +788,8 @@ def get_default_reach_dtype(structured=True):
                     ("thti", np.float32),
                     ("eps", np.float32),
                     ("uhc", np.float32),
-                    ("reachID", np.int),
-                    ("outreach", np.int),
+                    ("reachID", int),
+                    ("outreach", int),
                 ]
             )
 
@@ -797,12 +797,12 @@ def get_default_reach_dtype(structured=True):
     def get_default_segment_dtype():
         return np.dtype(
             [
-                ("nseg", np.int),
-                ("icalc", np.int),
-                ("outseg", np.int),
-                ("iupseg", np.int),
-                ("iprior", np.int),
-                ("nstrpts", np.int),
+                ("nseg", int),
+                ("icalc", int),
+                ("outseg", int),
+                ("iupseg", int),
+                ("iprior", int),
+                ("nstrpts", int),
                 ("flow", np.float32),
                 ("runoff", np.float32),
                 ("etsw", np.float32),
diff --git a/flopy/modflow/mfstr.py b/flopy/modflow/mfstr.py
index 7fbc1fe412..6a9cebfbf5 100644
--- a/flopy/modflow/mfstr.py
+++ b/flopy/modflow/mfstr.py
@@ -508,11 +508,11 @@ def get_default_dtype(structured=True):
         if structured:
             dtype = np.dtype(
                 [
-                    ("k", np.int),
-                    ("i", np.int),
-                    ("j", np.int),
-                    ("segment", np.int),
-                    ("reach", np.int),
+                    ("k", int),
+                    ("i", int),
+                    ("j", int),
+                    ("segment", int),
+                    ("reach", int),
                     ("flow", np.float32),
                     ("stage", np.float32),
                     ("cond", np.float32),
@@ -526,9 +526,9 @@ def get_default_dtype(structured=True):
         else:
             dtype = np.dtype(
                 [
-                    ("node", np.int),
-                    ("segment", np.int),
-                    ("reach", np.int),
+                    ("node", int),
+                    ("segment", int),
+                    ("reach", int),
                     ("flow", np.float32),
                     ("stage", np.float32),
                     ("cond", np.float32),
@@ -542,17 +542,17 @@ def get_default_dtype(structured=True):
 
         dtype2 = np.dtype(
             [
-                ("itrib01", np.int),
-                ("itrib02", np.int),
-                ("itrib03", np.int),
-                ("itrib04", np.int),
-                ("itrib05", np.int),
-                ("itrib06", np.int),
-                ("itrib07", np.int),
-                ("itrib08", np.int),
-                ("itrib09", np.int),
-                ("itrib10", np.int),
-                ("iupseg", np.int),
+                ("itrib01", int),
+                ("itrib02", int),
+                ("itrib03", int),
+                ("itrib04", int),
+                ("itrib05", int),
+                ("itrib06", int),
+                ("itrib07", int),
+                ("itrib08", int),
+                ("itrib09", int),
+                ("itrib10", int),
+                ("iupseg", int),
             ]
         )
         return dtype, dtype2
@@ -890,14 +890,14 @@ def load(cls, f, model, nper=None, ext_unit_dict=None):
 
                         #  get appropriate parval
                         if model.mfpar.pval is None:
-                            parval = np.float(par_dict["parval"])
+                            parval = float(par_dict["parval"])
                         else:
                             try:
-                                parval = np.float(
+                                parval = float(
                                     model.mfpar.pval.pval_dict[pname]
                                 )
                             except:
-                                parval = np.float(par_dict["parval"])
+                                parval = float(par_dict["parval"])
 
                         # fill current parameter data (par_current)
                         for ibnd, t in enumerate(data_dict):
diff --git a/flopy/modflow/mfwel.py b/flopy/modflow/mfwel.py
index 30de6f114f..22fa8cc084 100644
--- a/flopy/modflow/mfwel.py
+++ b/flopy/modflow/mfwel.py
@@ -230,8 +230,8 @@ def __init__(
                 if "specify" in opt:
                     t = opt.strip().split()
                     self.specify = True
-                    self.phiramp = np.float(t[1])
-                    self.iunitramp = np.int(t[2])
+                    self.phiramp = float(t[1])
+                    self.iunitramp = int(t[2])
                     self.options.pop(idx)
                     break
 
@@ -384,14 +384,14 @@ def get_default_dtype(structured=True):
         if structured:
             dtype = np.dtype(
                 [
-                    ("k", np.int),
-                    ("i", np.int),
-                    ("j", np.int),
+                    ("k", int),
+                    ("i", int),
+                    ("j", int),
                     ("flux", np.float32),
                 ]
             )
         else:
-            dtype = np.dtype([("node", np.int), ("flux", np.float32)])
+            dtype = np.dtype([("node", int), ("flux", np.float32)])
         return dtype
 
     @staticmethod
diff --git a/flopy/modpath/mp6.py b/flopy/modpath/mp6.py
index 931156413e..9ac152f6bc 100644
--- a/flopy/modpath/mp6.py
+++ b/flopy/modpath/mp6.py
@@ -271,7 +271,7 @@ def create_mpsim(
             self.__mf.dis.nrow,
             self.__mf.dis.ncol,
         )
-        arr = np.zeros((nlay, nrow, ncol), dtype=np.int)
+        arr = np.zeros((nlay, nrow, ncol), dtype=int)
         group_name = []
         group_region = []
         group_placement = []
diff --git a/flopy/modpath/mp6sim.py b/flopy/modpath/mp6sim.py
index 9cade92e5b..5630896a78 100644
--- a/flopy/modpath/mp6sim.py
+++ b/flopy/modpath/mp6sim.py
@@ -448,12 +448,12 @@ def get_dtypes():
         """
         dtype = np.dtype(
             [
-                ("particleid", np.int),
-                ("particlegroup", np.int),
-                ("initialgrid", np.int),
-                ("k0", np.int),
-                ("i0", np.int),
-                ("j0", np.int),
+                ("particleid", int),
+                ("particlegroup", int),
+                ("initialgrid", int),
+                ("k0", int),
+                ("i0", int),
+                ("j0", int),
                 ("xloc0", np.float32),
                 ("yloc0", np.float32),
                 ("zloc0", np.float32),
diff --git a/flopy/modpath/mp7particledata.py b/flopy/modpath/mp7particledata.py
index 8e6dcda4b3..7ee9a0ab81 100644
--- a/flopy/modpath/mp7particledata.py
+++ b/flopy/modpath/mp7particledata.py
@@ -896,9 +896,9 @@ def __init__(self, subdivisiondata=None, nodes=None):
         if isinstance(subdivisiondata, (CellDataType, FaceDataType)):
             subdivisiondata = [subdivisiondata]
 
-        if isinstance(nodes, (int, np.int, np.int32, np.int64)):
+        if isinstance(nodes, (int, np.int32, np.int64)):
             nodes = [(nodes,)]
-        elif isinstance(nodes, (float, np.float, np.float32, np.float64)):
+        elif isinstance(nodes, (float, np.float32, np.float64)):
             msg = (
                 "{}: nodes is of type {} ".format(self.name, type(nodes))
                 + "but must be an int if a single value is passed"
diff --git a/flopy/mt3d/mtcts.py b/flopy/mt3d/mtcts.py
index ee599a8556..396d3bc088 100644
--- a/flopy/mt3d/mtcts.py
+++ b/flopy/mt3d/mtcts.py
@@ -274,8 +274,8 @@ def get_default_CTS_dtype(ncomp=1, iforce=0):
         raise NotImplementedError()
 
         # # Item 3
-        # type_list = [("icts", np.int), ("next", np.int), ("ninj", np.int),
-        #              ("itrtinj", np.int)]
+        # type_list = [("icts", int), ("next", int), ("ninj", int),
+        #              ("itrtinj", int)]
         #
         # # Create a list for storing items 5, 6, & 9
         # items_5_6_7_9_list = []
@@ -291,7 +291,7 @@ def get_default_CTS_dtype(ncomp=1, iforce=0):
         #     for comp in range(1, ncomp+1):
         #         ioptinj_name = "ioptinj{0:d}".format(comp)
         #         cmchginj_name = "cmchginj{0:d}".format(comp)
-        #         items_5_6_7_9_list.append((ioptinj_name, np.int))
+        #         items_5_6_7_9_list.append((ioptinj_name, int))
         #         items_5_6_7_9_list.append((cmchginj_name, np.float32))
         #
         #     if iforce == 0:
@@ -305,14 +305,14 @@ def get_default_CTS_dtype(ncomp=1, iforce=0):
         # type_list.append(items_5_6_7_9_list)
         #
         # # Now create a list for the records in Item 4
-        # ext_wels_list = [("kext", np.int), ("iext", np.int), ("jext", np.int),
-        #                  ("iwext", np.int)]
+        # ext_wels_list = [("kext", int), ("iext", int), ("jext", int),
+        #                  ("iwext", int)]
         #
         # type_list.append(ext_wels_list)
         #
         # # Now create a list for the records in Item 8
-        # inj_wels_list = [("kinj", np.int), ("iinj", np.int), ("jinj", np.int),
-        #                  ("iwinj", np.int)]
+        # inj_wels_list = [("kinj", int), ("iinj", int), ("jinj", int),
+        #                  ("iwinj", int)]
         # type_list.append(inj_wels_list)
         #
         # #
diff --git a/flopy/mt3d/mtlkt.py b/flopy/mt3d/mtlkt.py
index ffccb293a1..4b53bdfe17 100644
--- a/flopy/mt3d/mtlkt.py
+++ b/flopy/mt3d/mtlkt.py
@@ -532,8 +532,8 @@ def get_default_dtype(ncomp=1):
         conditions interacting with the lake (i.e., pumps, specified runoff...)
         """
         type_list = [
-            ("node", np.int),
-            ("ilkbctyp", np.int),
+            ("node", int),
+            ("ilkbctyp", int),
             ("cbclk0", np.float32),
         ]
         if ncomp > 1:
diff --git a/flopy/mt3d/mtsft.py b/flopy/mt3d/mtsft.py
index 9721f8ae6e..85273d4ece 100644
--- a/flopy/mt3d/mtsft.py
+++ b/flopy/mt3d/mtsft.py
@@ -367,8 +367,8 @@ def get_default_dtype(ncomp=1):
         water boundary conditions.
         """
         type_list = [
-            ("node", np.int),
-            ("isfbctyp", np.int),
+            ("node", int),
+            ("isfbctyp", int),
             ("cbcsf0", np.float32),
         ]
         if ncomp > 1:
diff --git a/flopy/mt3d/mtssm.py b/flopy/mt3d/mtssm.py
index 2e6cfff553..0b59d2050c 100644
--- a/flopy/mt3d/mtssm.py
+++ b/flopy/mt3d/mtssm.py
@@ -43,7 +43,7 @@ class Mt3dSsm(Package):
         can be specified for additional species by passing additional
         arguments to the Mt3dSsm constructor.  For example, to specify the
         recharge concentration for species two one could use
-        crch2={0: 0., 1: 10*np.ones((nrow, ncol), dtype=np.float)} as
+        crch2={0: 0., 1: 10*np.ones((nrow, ncol), dtype=float)} as
         and additional keyword argument that is passed to Mt3dSsm when making
         the ssm object.
     cevt : Transient2d, scalar, array of floats, or dictionary
@@ -63,9 +63,9 @@ class Mt3dSsm(Package):
         values in the dictionary are recarrays of SSM boundaries.  The
         dtype for the recarray can be obtained using ssm.dtype (after the
         ssm package has been created).  The default dtype for the recarray is
-        np.dtype([('k', np.int), ("i", np.int), ("j", np.int),
-        ("css", np.float32), ("itype", np.int),
-        ((cssms(n), np.float), n=1, ncomp)])
+        np.dtype([('k', int), ("i", int), ("j", int),
+        ("css", np.float32), ("itype", int),
+        ((cssms(n), float), n=1, ncomp)])
         If there are more than one component species, then additional entries
         will be added to the dtype as indicated by cssm(n).
         Note that if the number of dictionary entries is less than the number
@@ -445,11 +445,11 @@ def get_default_dtype(ncomp=1):
         and sinks
         """
         type_list = [
-            ("k", np.int),
-            ("i", np.int),
-            ("j", np.int),
+            ("k", int),
+            ("i", int),
+            ("j", int),
             ("css", np.float32),
-            ("itype", np.int),
+            ("itype", int),
         ]
         if ncomp > 1:
             for comp in range(1, ncomp + 1):
diff --git a/flopy/pakbase.py b/flopy/pakbase.py
index db14f400fe..94804733dd 100644
--- a/flopy/pakbase.py
+++ b/flopy/pakbase.py
@@ -899,10 +899,10 @@ def load(f, model, pak_type, ext_unit_dict=None, **kwargs):
         nppak = 0
         if "parameter" in line.lower():
             t = line.strip().split()
-            nppak = np.int(t[1])
+            nppak = int(t[1])
             mxl = 0
             if nppak > 0:
-                mxl = np.int(t[2])
+                mxl = int(t[2])
                 if model.verbose:
                     msg = (
                         3 * " "
@@ -931,7 +931,7 @@ def load(f, model, pak_type, ext_unit_dict=None, **kwargs):
                     msg = 3 * " " + "implicit nppak in {}".format(filename)
                     print(msg)
             if nppak > 0:
-                mxl = np.int(t[3])
+                mxl = int(t[3])
                 imax += 1
                 if model.verbose:
                     msg = (
@@ -1089,12 +1089,12 @@ def load(f, model, pak_type, ext_unit_dict=None, **kwargs):
 
                 #  get appropriate parval
                 if model.mfpar.pval is None:
-                    parval = np.float(par_dict["parval"])
+                    parval = float(par_dict["parval"])
                 else:
                     try:
-                        parval = np.float(model.mfpar.pval.pval_dict[pname])
+                        parval = float(model.mfpar.pval.pval_dict[pname])
                     except:
-                        parval = np.float(par_dict["parval"])
+                        parval = float(par_dict["parval"])
 
                 # fill current parameter data (par_current)
                 for ibnd, t in enumerate(data_dict):
diff --git a/flopy/plot/crosssection.py b/flopy/plot/crosssection.py
index 6712fa5938..ed5b08b0f5 100644
--- a/flopy/plot/crosssection.py
+++ b/flopy/plot/crosssection.py
@@ -207,7 +207,7 @@ def __init__(
             if l > 0:
                 self.ncb += 1
 
-        self.active = np.ones((self.mg.nlay + self.ncb), dtype=np.int)
+        self.active = np.ones((self.mg.nlay + self.ncb), dtype=int)
         kon = 0
 
         if len(self.laycbd) > 0:
@@ -379,7 +379,7 @@ def plot_array(self, a, masked_values=None, head=None, **kwargs):
             )
             if len(self.laycbd) > 0:
                 if self.laycbd[k] > 0:
-                    ta = np.empty((self.mg.nrow, self.mg.ncol), dtype=np.float)
+                    ta = np.empty((self.mg.nrow, self.mg.ncol), dtype=float)
                     ta[:, :] = -1e9
                     vpts.append(
                         plotutil.cell_value_points(self.xpts, xedge, yedge, ta)
@@ -513,7 +513,7 @@ def plot_fill_between(
             )
             if len(self.laycbd) > 0:
                 if self.laycbd[k] > 0:
-                    ta = np.empty((self.mg.nrow, self.mg.ncol), dtype=np.float)
+                    ta = np.empty((self.mg.nrow, self.mg.ncol), dtype=float)
                     ta[:, :] = self.mg.botm.array[k, :, :]
                     vpts.append(
                         plotutil.cell_value_points(
diff --git a/flopy/plot/map.py b/flopy/plot/map.py
index eb6b4d47c8..d535c89677 100644
--- a/flopy/plot/map.py
+++ b/flopy/plot/map.py
@@ -321,7 +321,7 @@ def plot_inactive(self, ibound=None, color_noflow="black", **kwargs):
 
             ibound = self.mg.idomain
 
-        plotarray = np.zeros(ibound.shape, dtype=np.int)
+        plotarray = np.zeros(ibound.shape, dtype=int)
         idx1 = ibound == 0
         plotarray[idx1] = 1
         plotarray = np.ma.masked_equal(plotarray, 0)
@@ -373,7 +373,7 @@ def plot_ibound(
 
             ibound = self.mg.idomain
 
-        plotarray = np.zeros(ibound.shape, dtype=np.int)
+        plotarray = np.zeros(ibound.shape, dtype=int)
         idx1 = ibound == 0
         idx2 = ibound < 0
         plotarray[idx1] = 1
@@ -528,9 +528,9 @@ def plot_bc(
             nlay = self.mg.nlay
 
         # Plot the list locations
-        plotarray = np.zeros(self.mg.shape, dtype=np.int)
+        plotarray = np.zeros(self.mg.shape, dtype=int)
         if plotAll:
-            pa = np.zeros(self.mg.shape[1:], dtype=np.int)
+            pa = np.zeros(self.mg.shape[1:], dtype=int)
             pa[tuple(idx[1:])] = 1
             for k in range(nlay):
                 plotarray[k] = pa.copy()
@@ -640,7 +640,7 @@ def contour_array_cvfd(self, vertc, a, masked_values=None, **kwargs):
             ncpl = kwargs.pop("ncpl")
             if isinstance(ncpl, int):
                 i = int(ncpl)
-                ncpl = np.ones((nlay,), dtype=np.int) * i
+                ncpl = np.ones((nlay,), dtype=int) * i
             elif isinstance(ncpl, list) or isinstance(ncpl, tuple):
                 ncpl = np.array(ncpl)
             i0 = 0
@@ -975,7 +975,7 @@ def plot_discharge(
                     laycbd = self.model.laycbd
 
             if laycbd is not None and 1 in laycbd:
-                active = np.ones((botm.shape[0],), dtype=np.int)
+                active = np.ones((botm.shape[0],), dtype=int)
                 kon = 0
                 for cbd in laycbd:
                     if cbd > 0:
@@ -987,7 +987,7 @@ def plot_discharge(
             # thickness by setting laytyp to zeros
             if head is None or laytyp is None:
                 head = np.zeros(botm.shape, np.float32)
-                laytyp = np.zeros((botm.shape[0],), dtype=np.int)
+                laytyp = np.zeros((botm.shape[0],), dtype=int)
 
             # calculate the saturated thickness
             sat_thk = plotutil.PlotUtilities.saturated_thickness(
@@ -1010,7 +1010,7 @@ def plot_discharge(
 
             spdis = np.recarray(
                 (len(temp),),
-                dtype=[("node", np.int), ("qx", np.float), ("qy", np.float)],
+                dtype=[("node", int), ("qx", float), ("qy", float)],
             )
             for ix, tup in enumerate(temp):
                 spdis[ix] = tup
diff --git a/flopy/plot/plotbase.py b/flopy/plot/plotbase.py
index c00c0b3d4b..cdc3ece80b 100644
--- a/flopy/plot/plotbase.py
+++ b/flopy/plot/plotbase.py
@@ -243,7 +243,7 @@ def plot_inactive(self, ibound=None, color_noflow="black", **kwargs):
             else:
                 ibound = self.mg.idomain
 
-        plotarray = np.zeros(ibound.shape, dtype=np.int)
+        plotarray = np.zeros(ibound.shape, dtype=int)
         idx1 = ibound == 0
         plotarray[idx1] = 1
         plotarray = np.ma.masked_equal(plotarray, 0)
@@ -297,7 +297,7 @@ def plot_ibound(
 
             ibound = self.mg.idomain
 
-        plotarray = np.zeros(ibound.shape, dtype=np.int)
+        plotarray = np.zeros(ibound.shape, dtype=int)
         idx1 = ibound == 0
         idx2 = ibound < 0
         plotarray[idx1] = 1
@@ -439,11 +439,11 @@ def plot_bc(
 
         # Plot the list locations, change this to self.mg.shape
         if len(self.mg.shape) != 3:
-            plotarray = np.zeros((self.mg.nlay, self.mg.ncpl), dtype=np.int)
+            plotarray = np.zeros((self.mg.nlay, self.mg.ncpl), dtype=int)
             plotarray[tuple(idx)] = 1
         else:
             plotarray = np.zeros(
-                (self.mg.nlay, self.mg.nrow, self.mg.ncol), dtype=np.int
+                (self.mg.nlay, self.mg.nrow, self.mg.ncol), dtype=int
             )
             plotarray[idx[0], idx[1], idx[2]] = 1
 
@@ -976,7 +976,7 @@ def plot_discharge(
             # thickness by setting laytyp to zeros
             if head is None or laytyp is None:
                 head = np.zeros(botm.shape, np.float32)
-                laytyp = np.zeros((nlay), dtype=np.int)
+                laytyp = np.zeros((nlay), dtype=int)
                 head[0, :, :] = top
                 if nlay > 1:
                     head[1:, :, :] = botm[:-1, :, :]
@@ -991,7 +991,7 @@ def plot_discharge(
             )
 
             if qz is None:
-                qz = np.zeros((qx.shape), dtype=np.float)
+                qz = np.zeros((qx.shape), dtype=float)
 
             ib = ib.ravel()
             qx = qx.ravel()
@@ -1006,10 +1006,10 @@ def plot_discharge(
             spdis = np.recarray(
                 (len(temp),),
                 dtype=[
-                    ("node", np.int),
-                    ("qx", np.float),
-                    ("qy", np.float),
-                    ("qz", np.float),
+                    ("node", int),
+                    ("qx", float),
+                    ("qy", float),
+                    ("qz", float),
                 ],
             )
             for ix, tup in enumerate(temp):
diff --git a/flopy/plot/plotutil.py b/flopy/plot/plotutil.py
index ca6e80164e..32949458f7 100644
--- a/flopy/plot/plotutil.py
+++ b/flopy/plot/plotutil.py
@@ -1915,9 +1915,7 @@ def __init__(self, model=None, botm=None, istrat=1, nu=None):
                 dis = model.get_package("DIS")
             except:
                 sys.stdout.write("Error: DIS package not available.\n")
-            self.__botm = np.zeros(
-                (dis.nlay + 1, dis.nrow, dis.ncol), np.float
-            )
+            self.__botm = np.zeros((dis.nlay + 1, dis.nrow, dis.ncol), float)
             self.__botm[0, :, :] = dis.top.array
             self.__botm[1:, :, :] = dis.botm.array
             try:
@@ -1959,7 +1957,7 @@ def calc_conc(self, zeta, layer=None):
         >>> conc = c.calc_conc(z, layer=0)
 
         """
-        conc = np.zeros((self.__nlay, self.__nrow, self.__ncol), np.float)
+        conc = np.zeros((self.__nlay, self.__nrow, self.__ncol), float)
 
         pct = {}
         for isrf in range(self.__nsrf):
@@ -2340,7 +2338,7 @@ def plot_cvfd(
         ncpl = kwargs.pop("ncpl")
         if isinstance(ncpl, int):
             i = int(ncpl)
-            ncpl = np.ones((nlay), dtype=np.int) * i
+            ncpl = np.ones((nlay), dtype=int) * i
         elif isinstance(ncpl, list) or isinstance(ncpl, tuple):
             ncpl = np.array(ncpl)
         i0 = 0
diff --git a/flopy/utils/binaryfile.py b/flopy/utils/binaryfile.py
index d4e4692b7e..f06cf5dc7b 100755
--- a/flopy/utils/binaryfile.py
+++ b/flopy/utils/binaryfile.py
@@ -413,10 +413,10 @@ def get_ts(self, idx):
                 )  # change ilay from header to zero-based
                 if ilay != k:
                     continue
-                ipos = np.long(self.iposarray[irec])
+                ipos = int(self.iposarray[irec])
 
                 # Calculate offset necessary to reach intended cell
-                self.file.seek(ipos + np.long(ioffset), 0)
+                self.file.seek(ipos + int(ioffset), 0)
 
                 # Find the time index and then put value into result in the
                 # correct location.
@@ -1526,7 +1526,7 @@ def get_record(self, idx, full3D=False):
             idx = np.array([idx])
 
         header = self.recordarray[idx]
-        ipos = np.long(self.iposarray[idx])
+        ipos = int(self.iposarray[idx])
         self.file.seek(ipos, 0)
         imeth = header["imeth"][0]
 
@@ -1755,9 +1755,9 @@ def get_residual(self, totim, scaled=False):
         nlay = self.nlay
         nrow = self.nrow
         ncol = self.ncol
-        residual = np.zeros((nlay, nrow, ncol), dtype=np.float)
+        residual = np.zeros((nlay, nrow, ncol), dtype=float)
         if scaled:
-            inflow = np.zeros((nlay, nrow, ncol), dtype=np.float)
+            inflow = np.zeros((nlay, nrow, ncol), dtype=float)
         select_indices = np.where((self.recordarray["totim"] == totim))[0]
 
         for i in select_indices:
@@ -1802,7 +1802,7 @@ def get_residual(self, totim, scaled=False):
                     inflow[idx] += flow[idx]
 
         if scaled:
-            residual_scaled = np.zeros((nlay, nrow, ncol), dtype=np.float)
+            residual_scaled = np.zeros((nlay, nrow, ncol), dtype=float)
             idx = inflow > 0.0
             residual_scaled[idx] = residual[idx] / inflow[idx]
             return residual_scaled
diff --git a/flopy/utils/check.py b/flopy/utils/check.py
index 13cd0bbb91..51ecd382ae 100644
--- a/flopy/utils/check.py
+++ b/flopy/utils/check.py
@@ -646,23 +646,23 @@ def _get_dtype(self):
             # include node column for structured grids (useful for indexing)
             return np.dtype(
                 [
-                    ("type", np.object),
-                    ("package", np.object),
-                    ("k", np.int),
-                    ("i", np.int),
-                    ("j", np.int),
-                    ("value", np.float),
-                    ("desc", np.object),
+                    ("type", object),
+                    ("package", object),
+                    ("k", int),
+                    ("i", int),
+                    ("j", int),
+                    ("value", float),
+                    ("desc", object),
                 ]
             )
         else:
             return np.dtype(
                 [
-                    ("type", np.object),
-                    ("package", np.object),
-                    ("node", np.int),
-                    ("value", np.float),
-                    ("desc", np.object),
+                    ("type", object),
+                    ("package", object),
+                    ("node", int),
+                    ("value", float),
+                    ("desc", object),
                 ]
             )
 
@@ -814,34 +814,34 @@ def _get_dtype(self):
         if mg.grid_type == "structured":
             return np.dtype(
                 [
-                    ("type", np.object),
-                    ("package", np.object),
-                    ("k", np.int),
-                    ("i", np.int),
-                    ("j", np.int),
-                    ("value", np.float),
-                    ("desc", np.object),
+                    ("type", object),
+                    ("package", object),
+                    ("k", int),
+                    ("i", int),
+                    ("j", int),
+                    ("value", float),
+                    ("desc", object),
                 ]
             )
         elif mg.grid_type == "vertex":
             return np.dtype(
                 [
-                    ("type", np.object),
-                    ("package", np.object),
-                    ("lay", np.int),
-                    ("cell", np.int),
-                    ("value", np.float),
-                    ("desc", np.object),
+                    ("type", object),
+                    ("package", object),
+                    ("lay", int),
+                    ("cell", int),
+                    ("value", float),
+                    ("desc", object),
                 ]
             )
         else:
             return np.dtype(
                 [
-                    ("type", np.object),
-                    ("package", np.object),
-                    ("node", np.int),
-                    ("value", np.float),
-                    ("desc", np.object),
+                    ("type", object),
+                    ("package", object),
+                    ("node", int),
+                    ("value", float),
+                    ("desc", object),
                 ]
             )
 
diff --git a/flopy/utils/cvfdutil.py b/flopy/utils/cvfdutil.py
index 8fde403f65..a8fa19910c 100644
--- a/flopy/utils/cvfdutil.py
+++ b/flopy/utils/cvfdutil.py
@@ -176,7 +176,7 @@ def to_cvfd(
     # In the process, filter out any duplicate vertices
     vertexdict = OrderedDict()
     vertexlist = []
-    xcyc = np.empty((ncells, 2), dtype=np.float)
+    xcyc = np.empty((ncells, 2), dtype=float)
     iv = 0
     nvertstart = 0
     if verbose:
@@ -304,7 +304,7 @@ def shapefile_to_xcyc(shp):
     sf = shapefile.Reader(shp)
     shapes = sf.shapes()
     ncells = len(shapes)
-    xcyc = np.empty((ncells, 2), dtype=np.float)
+    xcyc = np.empty((ncells, 2), dtype=float)
     for icell, shape in enumerate(shapes):
         points = shape.points
         xc, yc = centroid_of_polygon(points)
@@ -368,7 +368,7 @@ def get_disv_gridprops(verts, iverts, xcyc=None):
     nvert = verts.shape[0]
     ncpl = len(iverts)
     if xcyc is None:
-        xcyc = np.empty((ncpl, 2), dtype=np.float)
+        xcyc = np.empty((ncpl, 2), dtype=float)
         for icell in range(ncpl):
             vlist = [
                 (verts[ivert, 0], verts[ivert, 1]) for ivert in iverts[icell]
diff --git a/flopy/utils/flopy_io.py b/flopy/utils/flopy_io.py
index fc7e324eef..a7d4ac0db6 100755
--- a/flopy/utils/flopy_io.py
+++ b/flopy/utils/flopy_io.py
@@ -187,14 +187,14 @@ def write_fixed_var(v, length=10, ipos=None, free=False, comment=None):
         if free:
             write_fmt = "{} "
         else:
-            if isinstance(v[n], (float, np.float, np.float32, np.float64)):
+            if isinstance(v[n], (float, np.float32, np.float64)):
                 width = ipos[n] - 6
                 vmin, vmax = 10 ** -width, 10 ** width
                 if abs(v[n]) < vmin or abs(v[n]) > vmax:
                     ctype = "g"
                 else:
                     ctype = ".{}f".format(width)
-            elif isinstance(v[n], (int, np.int, np.int32, np.int64)):
+            elif isinstance(v[n], (int, np.int32, np.int64)):
                 ctype = "d"
             else:
                 ctype = ""
diff --git a/flopy/utils/gridgen.py b/flopy/utils/gridgen.py
index 71bc67fc81..019b4da3a9 100644
--- a/flopy/utils/gridgen.py
+++ b/flopy/utils/gridgen.py
@@ -112,10 +112,10 @@ def ndarray_to_asciigrid(fname, a, extent, nodata=1.0e30):
     header += "xllcorner {}\n".format(xmin)
     header += "yllcorner {}\n".format(ymin)
     header += "cellsize {}\n".format(dx)
-    header += "NODATA_value {}\n".format(np.float(nodata))
+    header += "NODATA_value {}\n".format(float(nodata))
     # replace nan with nodata
     idx = np.isnan(a)
-    a[idx] = np.float(nodata)
+    a[idx] = float(nodata)
     # write
     with open(fname, "wb") as f:
         f.write(header.encode("ascii"))
@@ -216,7 +216,7 @@ def __init__(
 
         self.nodes = 0
         self.nja = 0
-        self.nodelay = np.zeros((self.nlay), dtype=np.int)
+        self.nodelay = np.zeros((self.nlay), dtype=int)
         self._vertdict = {}
         self.model_ws = model_ws
         exe_name = which(exe_name)
@@ -727,14 +727,14 @@ def get_nod_recarray(self):
         f = open(fname, "r")
         dt = np.dtype(
             [
-                ("node", np.int),
-                ("layer", np.int),
-                ("x", np.float),
-                ("y", np.float),
-                ("z", np.float),
-                ("dx", np.float),
-                ("dy", np.float),
-                ("dz", np.float),
+                ("node", int),
+                ("layer", int),
+                ("x", float),
+                ("y", float),
+                ("z", float),
+                ("dx", float),
+                ("dy", float),
+                ("dz", float),
             ]
         )
         node_ra = np.genfromtxt(fname, dtype=dt, skip_header=1)
@@ -799,7 +799,7 @@ def get_disu(
         self.nodes = nodes
 
         # nodelay
-        nodelay = np.empty((nlay), dtype=np.int)
+        nodelay = np.empty((nlay), dtype=int)
         fname = os.path.join(self.model_ws, "qtg.nodesperlay.dat")
         f = open(fname, "r")
         nodelay = read1d(f, nodelay)
@@ -874,7 +874,7 @@ def get_disu(
             istart = istop
 
         # iac
-        iac = np.empty((nodes), dtype=np.int)
+        iac = np.empty((nodes), dtype=int)
         fname = os.path.join(self.model_ws, "qtg.iac.dat")
         f = open(fname, "r")
         iac = read1d(f, iac)
@@ -886,14 +886,14 @@ def get_disu(
 
         # ja -- this is being read is as one-based, which is also what is
         # expected by the ModflowDisu constructor
-        ja = np.empty((njag), dtype=np.int)
+        ja = np.empty((njag), dtype=int)
         fname = os.path.join(self.model_ws, "qtg.ja.dat")
         f = open(fname, "r")
         ja = read1d(f, ja)
         f.close()
 
         # ivc
-        fldr = np.empty((njag), dtype=np.int)
+        fldr = np.empty((njag), dtype=int)
         fname = os.path.join(self.model_ws, "qtg.fldr.dat")
         f = open(fname, "r")
         fldr = read1d(f, fldr)
@@ -988,7 +988,7 @@ def get_nodelay(self):
 
         """
         nlay = self.get_nlay()
-        nodelay = np.empty((nlay), dtype=np.int)
+        nodelay = np.empty((nlay), dtype=int)
         fname = os.path.join(self.model_ws, "qtg.nodesperlay.dat")
         f = open(fname, "r")
         nodelay = read1d(f, nodelay)
@@ -1080,7 +1080,7 @@ def get_iac(self):
 
         """
         nodes = self.get_nodes()
-        iac = np.empty((nodes), dtype=np.int)
+        iac = np.empty((nodes), dtype=int)
         fname = os.path.join(self.model_ws, "qtg.iac.dat")
         f = open(fname, "r")
         iac = read1d(f, iac)
@@ -1106,7 +1106,7 @@ def get_ja(self, nja=None):
         if nja is None:
             iac = self.get_iac()
             nja = iac.sum()
-        ja = np.empty((nja), dtype=np.int)
+        ja = np.empty((nja), dtype=int)
         fname = os.path.join(self.model_ws, "qtg.ja.dat")
         f = open(fname, "r")
         ja = read1d(f, ja)
@@ -1128,7 +1128,7 @@ def get_fldr(self):
         """
         iac = self.get_iac()
         njag = iac.sum()
-        fldr = np.empty((njag), dtype=np.int)
+        fldr = np.empty((njag), dtype=int)
         fname = os.path.join(self.model_ws, "qtg.fldr.dat")
         f = open(fname, "r")
         fldr = read1d(f, fldr)
@@ -1154,7 +1154,7 @@ def get_ivc(self, fldr=None):
         """
         if fldr is None:
             fldr = self.get_fldr()
-        ivc = np.zeros(fldr.shape, dtype=np.int)
+        ivc = np.zeros(fldr.shape, dtype=int)
         idx = abs(fldr) == 3
         ivc[idx] = 1
         return ivc
@@ -1179,7 +1179,7 @@ def get_ihc(self, nodelay=None, ia=None, fldr=None):
         """
         if fldr is None:
             fldr = self.get_fldr()
-        ihc = np.empty(fldr.shape, dtype=np.int)
+        ihc = np.empty(fldr.shape, dtype=int)
         ihc = np.where(abs(fldr) == 0, 0, ihc)
         ihc = np.where(abs(fldr) == 1, 1, ihc)
         ihc = np.where(abs(fldr) == 2, 1, ihc)
@@ -1193,7 +1193,7 @@ def get_ihc(self, nodelay=None, ia=None, fldr=None):
             ia = get_ia_from_iac(iac)
         nodes = ia.shape[0] - 1
         nlayers = nodelay.shape[0]
-        layers = -1 * np.ones(nodes, dtype=np.int)
+        layers = -1 * np.ones(nodes, dtype=int)
         node_layer_range = [0] + list(np.add.accumulate(nodelay))
         for ilay in range(nlayers):
             istart = node_layer_range[ilay]
@@ -1324,7 +1324,7 @@ def get_angldegx(self, fldr=None):
         """
         if fldr is None:
             fldr = self.get_fldr()
-        angldegx = np.zeros(fldr.shape, dtype=np.float)
+        angldegx = np.zeros(fldr.shape, dtype=float)
         angldegx = np.where(fldr == 0, 1.0e30, angldegx)
         angldegx = np.where(abs(fldr) == 3, 1.0e30, angldegx)
         angldegx = np.where(fldr == 2, 90, angldegx)
@@ -1373,7 +1373,7 @@ def get_cellxy(self, ncells):
             x and y cell centers.  Shape is (ncells, 2)
 
         """
-        cellxy = np.empty((ncells, 2), dtype=np.float)
+        cellxy = np.empty((ncells, 2), dtype=float)
         for n in range(ncells):
             x, y = self.get_center(n)
             cellxy[n, 0] = x
diff --git a/flopy/utils/lgrutil.py b/flopy/utils/lgrutil.py
index ae122ff2fe..d2eb6cec73 100644
--- a/flopy/utils/lgrutil.py
+++ b/flopy/utils/lgrutil.py
@@ -77,11 +77,11 @@ def __init__(
 
         # # child cells per parent and child cells per parent layer
         self.ncpp = ncpp
-        self.ncppl = Util2d(m, (nlayp,), np.int, ncppl, "ncppl").array
+        self.ncppl = Util2d(m, (nlayp,), int, ncppl, "ncppl").array
 
         # calculate ibcl which is the bottom child layer (one based) in each
         # parent layer
-        self.ibcl = np.zeros(self.nlayp, dtype=np.int)
+        self.ibcl = np.zeros(self.nlayp, dtype=int)
         self.ibcl[0] = self.ncppl[0]
         for k in range(1, self.nlayp):
             if self.ncppl[k] > 0:
@@ -231,7 +231,7 @@ def get_idomain(self):
             idomain array for the child model
 
         """
-        idomain = np.ones((self.nlay, self.nrow, self.ncol), dtype=np.int)
+        idomain = np.ones((self.nlay, self.nrow, self.ncol), dtype=int)
         for kc in range(self.nlay):
             for ic in range(self.nrow):
                 for jc in range(self.ncol):
diff --git a/flopy/utils/observationfile.py b/flopy/utils/observationfile.py
index 3eebc4abb5..2802a12f5c 100644
--- a/flopy/utils/observationfile.py
+++ b/flopy/utils/observationfile.py
@@ -296,7 +296,7 @@ def __init__(self, filename, verbose=False, isBinary=True):
             self.nobs = self.read_integer()
 
             # # continue reading the file
-            # self.v = np.empty(self.nobs, dtype=np.float)
+            # self.v = np.empty(self.nobs, dtype=float)
             # self.v.fill(1.0E+32)
 
             # read obsnames
@@ -401,7 +401,7 @@ def __init__(self, filename, verbose=False, hydlbl_len=20):
 
         # continue reading the file
         self.itmuni = self.read_integer()
-        self.v = np.empty(self.nobs, dtype=np.float)
+        self.v = np.empty(self.nobs, dtype=float)
         self.v.fill(1.0e32)
         ctime = self.read_text(nchar=4)
         self.hydlbl_len = int(hydlbl_len)
diff --git a/flopy/utils/optionblock.py b/flopy/utils/optionblock.py
index 5926414f3d..07e9e3fc1e 100644
--- a/flopy/utils/optionblock.py
+++ b/flopy/utils/optionblock.py
@@ -168,7 +168,7 @@ def __setattr__(self, key, value):
         elif isinstance(value, np.recarray):
             for name in value.dtype.names:
                 if self._attr_types[name] == np.bool_:
-                    if not isinstance(value, (bool, np.bool_, np.bool)):
+                    if not isinstance(value, (bool, np.bool_)):
                         raise TypeError(err_msg.format(self._attr_types[name]))
                 else:
                     try:
@@ -180,7 +180,7 @@ def __setattr__(self, key, value):
 
         elif key in self._attr_types:
             if self._attr_types[key] == np.bool_:
-                if not isinstance(value, (bool, np.bool_, np.bool)):
+                if not isinstance(value, (bool, np.bool_)):
                     raise TypeError(err_msg.format(self._attr_types[key]))
             else:
                 try:
@@ -249,14 +249,14 @@ def _set_attributes(self):
         """
         # set up all attributes for the class!
         for key, ctx in self._context.items():
-            if ctx[OptionBlock.dtype] in (np.bool_, bool, np.bool):
+            if ctx[OptionBlock.dtype] in (np.bool_, bool):
                 self.__setattr__(key, False)
             else:
                 self.__setattr__(key, None)
 
             if ctx[OptionBlock.nested]:
                 for k, d in ctx[OptionBlock.vars].items():
-                    if d[OptionBlock.dtype] in (np.bool_, bool, np.bool):
+                    if d[OptionBlock.dtype] in (np.bool_, bool):
                         self.__setattr__(k, False)
                     else:
                         self.__setattr__(k, None)
diff --git a/flopy/utils/rasters.py b/flopy/utils/rasters.py
index 8bf8bfb004..b5616610f9 100644
--- a/flopy/utils/rasters.py
+++ b/flopy/utils/rasters.py
@@ -51,11 +51,11 @@ class Raster(object):
 
     """
 
-    FLOAT32 = (float, np.float, np.float32, np.float_)
+    FLOAT32 = (float, np.float32, np.float_)
     FLOAT64 = (np.float64,)
     INT8 = (np.int8,)
     INT16 = (np.int16,)
-    INT32 = (int, np.int, np.int32, np.int_)
+    INT32 = (int, np.int32, np.int_)
     INT64 = (np.int64,)
 
     def __init__(
diff --git a/flopy/utils/reference.py b/flopy/utils/reference.py
index feb37d9b05..4c60be560f 100755
--- a/flopy/utils/reference.py
+++ b/flopy/utils/reference.py
@@ -1469,7 +1469,7 @@ def get_2d_vertex_connectivity(self):
         nrowvert = self.nrow + 1
         ncolvert = self.ncol + 1
         npoints = nrowvert * ncolvert
-        verts = np.empty((npoints, 2), dtype=np.float)
+        verts = np.empty((npoints, 2), dtype=float)
         verts[:, 0] = x
         verts[:, 1] = y
         iverts = []
@@ -1496,7 +1496,7 @@ def get_3d_shared_vertex_connectivity(self, nlay, botm, ibound=None):
         npoints = nrvncv * nlayvert
 
         # create and fill a 3d points array for the grid
-        verts = np.empty((npoints, 3), dtype=np.float)
+        verts = np.empty((npoints, 3), dtype=float)
         verts[:, 0] = np.tile(x, nlayvert)
         verts[:, 1] = np.tile(y, nlayvert)
         istart = 0
@@ -1539,7 +1539,7 @@ def get_3d_shared_vertex_connectivity(self, nlay, botm, ibound=None):
         if ibound is not None:
 
             # go through the vertex list and mark vertices that are used
-            ivertrenum = np.zeros(npoints, dtype=np.int)
+            ivertrenum = np.zeros(npoints, dtype=int)
             for vlist in iverts:
                 for iv in vlist:
                     # mark vertices that are actually used
@@ -1569,11 +1569,11 @@ def get_3d_shared_vertex_connectivity(self, nlay, botm, ibound=None):
     def get_3d_vertex_connectivity(self, nlay, top, bot, ibound=None):
         if ibound is None:
             ncells = nlay * self.nrow * self.ncol
-            ibound = np.ones((nlay, self.nrow, self.ncol), dtype=np.int)
+            ibound = np.ones((nlay, self.nrow, self.ncol), dtype=int)
         else:
             ncells = (ibound != 0).sum()
         npoints = ncells * 8
-        verts = np.empty((npoints, 3), dtype=np.float)
+        verts = np.empty((npoints, 3), dtype=float)
         iverts = []
         ipoint = 0
         for k in range(nlay):
@@ -1827,9 +1827,9 @@ def from_argus_export(cls, fname, nlay=1):
         ncells, nverts = ll[0:2]
         ncells = int(ncells)
         nverts = int(nverts)
-        verts = np.empty((nverts, 2), dtype=np.float)
-        xc = np.empty((ncells), dtype=np.float)
-        yc = np.empty((ncells), dtype=np.float)
+        verts = np.empty((nverts, 2), dtype=float)
+        xc = np.empty((ncells), dtype=float)
+        yc = np.empty((ncells), dtype=float)
 
         # read the vertices
         f.readline()
diff --git a/flopy/utils/swroutputfile.py b/flopy/utils/swroutputfile.py
index 7c36c98edf..0392f4aba7 100644
--- a/flopy/utils/swroutputfile.py
+++ b/flopy/utils/swroutputfile.py
@@ -335,7 +335,7 @@ def _read_connectivity(self):
         self.conn_dtype = np.dtype(
             [("reach", "i4"), ("from", "i4"), ("to", "i4")]
         )
-        conn = np.zeros((self.nrecord, 3), np.int)
+        conn = np.zeros((self.nrecord, 3), int)
         icount = 0
         for nrg in range(self.flowitems):
             flowitems = self.read_integer()
@@ -405,7 +405,7 @@ def _build_dtypes(self):
     def _read_header(self):
         nitems = 0
         if self.type == "exchange" or self.type == "structure":
-            itemlist = np.zeros(self.nrecord, np.int)
+            itemlist = np.zeros(self.nrecord, int)
             try:
                 for i in range(self.nrecord):
                     itemlist[i] = self.read_integer()
diff --git a/flopy/utils/triangle.py b/flopy/utils/triangle.py
index 697ee77bd2..9be313d95e 100644
--- a/flopy/utils/triangle.py
+++ b/flopy/utils/triangle.py
@@ -190,7 +190,7 @@ def build(self, verbose=False):
 
         # create verts and iverts
         self.verts = self.node[["x", "y"]]
-        self.verts = np.array(self.verts.tolist(), np.float)
+        self.verts = np.array(self.verts.tolist(), float)
         self.iverts = []
         for row in self.ele:
             self.iverts.append([row[1], row[2], row[3]])
@@ -279,7 +279,7 @@ def get_boundary_marker_array(self):
             polygon that is added with the add_polygon method.
 
         """
-        iedge = np.zeros((self.ncpl), dtype=np.int)
+        iedge = np.zeros((self.ncpl), dtype=int)
         boundary_markers = np.unique(self.edge["boundary_marker"])
         for ibm in boundary_markers:
             icells = self.get_edge_cells(ibm)
@@ -484,7 +484,7 @@ def get_xcyc(self):
 
         """
         ncpl = len(self.iverts)
-        xcyc = np.empty((ncpl, 2), dtype=np.float)
+        xcyc = np.empty((ncpl, 2), dtype=float)
         for i, icell2d in enumerate(self.iverts):
             points = []
             for iv in icell2d:
diff --git a/flopy/utils/util_array.py b/flopy/utils/util_array.py
index b9bf1c025e..f91394e10a 100644
--- a/flopy/utils/util_array.py
+++ b/flopy/utils/util_array.py
@@ -470,7 +470,7 @@ class Util3d(DataInterface):
         this package will be added.
     shape : length 3 tuple
         shape of the 3-D array, typically (nlay,nrow,ncol)
-    dtype : [np.int32, np.float32, np.bool]
+    dtype : [np.int32, np.float32, bool]
         the type of the data
     value : variable
         the data to be assigned to the 3-D array.
@@ -1016,7 +1016,7 @@ class Transient3d(DataInterface):
         this package will be added.
     shape : length 3 tuple
         shape of the 3-D transient arrays, typically (nlay,nrow,ncol)
-    dtype : [np.int32, np.float32, np.bool]
+    dtype : [np.int32, np.float32, bool]
         the type of the data
     value : variable
         the data to be assigned to the 3-D arrays. Typically a dict
@@ -1302,7 +1302,7 @@ class Transient2d(DataInterface):
         this package will be added.
     shape : length 2 tuple
         shape of the 2-D transient arrays, typically (nrow,ncol)
-    dtype : [np.int32, np.float32, np.bool]
+    dtype : [np.int32, np.float32, bool]
         the type of the data
     value : variable
         the data to be assigned to the 2-D arrays. Typically a dict
@@ -1804,7 +1804,7 @@ class Util2d(DataInterface):
         this package will be added.
     shape : tuple
         Shape of the 1- or 2-D array
-    dtype : [np.int32, np.float32, np.bool]
+    dtype : [np.int32, np.float32, bool]
         the type of the data
     value : variable
         the data to be assigned to the 1- or 2-D array.
@@ -1951,7 +1951,7 @@ def __init__(
                     )
                 )
             dtype = np.int32
-        if dtype not in [np.int32, np.float32, np.bool]:
+        if dtype not in [np.int32, np.float32, bool]:
             raise TypeError("Util2d:unsupported dtype: " + str(dtype))
 
         if name is not None:
@@ -2902,20 +2902,20 @@ def parse_value(self, value):
             value = np.array(value)
 
         if isinstance(value, bool):
-            if self._dtype == np.bool:
+            if self._dtype == bool:
                 try:
-                    self.__value = np.bool(value)
+                    self.__value = bool(value)
 
                 except:
                     raise Exception(
                         "Util2d:could not cast "
-                        + 'boolean value to type "np.bool": '
+                        + 'boolean value to type "bool": '
                         + str(value)
                     )
             else:
                 raise Exception(
                     "Util2d:value type is bool, "
-                    + " but dtype not set as np.bool"
+                    + " but dtype not set as bool"
                 )
         elif isinstance(value, str):
             if os.path.exists(value):
@@ -3161,21 +3161,21 @@ def parse_control_record(
         freefmt, cnstnt, fmtin, iprn, nunit = None, None, None, -1, None
         fname = None
         isfloat = False
-        if dtype == np.float or dtype == np.float32:
+        if dtype == float or dtype == np.float32:
             isfloat = True
             # if free format keywords
         if str(raw[0].lower()) in str(free_fmt):
             freefmt = raw[0].lower()
             if raw[0].lower() == "constant":
                 if isfloat:
-                    cnstnt = np.float(raw[1].lower().replace("d", "e"))
+                    cnstnt = float(raw[1].lower().replace("d", "e"))
                 else:
-                    cnstnt = np.int(raw[1].lower())
+                    cnstnt = int(raw[1].lower())
             if raw[0].lower() == "internal":
                 if isfloat:
-                    cnstnt = np.float(raw[1].lower().replace("d", "e"))
+                    cnstnt = float(raw[1].lower().replace("d", "e"))
                 else:
-                    cnstnt = np.int(raw[1].lower())
+                    cnstnt = int(raw[1].lower())
                 fmtin = raw[2].strip()
                 iprn = 0
                 if len(raw) >= 4:
@@ -3193,9 +3193,9 @@ def parse_control_record(
 
                 nunit = int(raw[1])
                 if isfloat:
-                    cnstnt = np.float(raw[2].lower().replace("d", "e"))
+                    cnstnt = float(raw[2].lower().replace("d", "e"))
                 else:
-                    cnstnt = np.int(raw[2].lower())
+                    cnstnt = int(raw[2].lower())
                 fmtin = raw[3].strip()
                 iprn = 0
                 if len(raw) >= 5:
@@ -3203,26 +3203,26 @@ def parse_control_record(
             elif raw[0].lower() == "open/close":
                 fname = raw[1].strip()
                 if isfloat:
-                    cnstnt = np.float(raw[2].lower().replace("d", "e"))
+                    cnstnt = float(raw[2].lower().replace("d", "e"))
                 else:
-                    cnstnt = np.int(raw[2].lower())
+                    cnstnt = int(raw[2].lower())
                 fmtin = raw[3].strip()
                 iprn = 0
                 if len(raw) >= 5:
                     iprn = int(raw[4])
                 npl, fmt, width, decimal = None, None, None, None
         else:
-            locat = np.int(line[0:10].strip())
+            locat = int(line[0:10].strip())
             if isfloat:
                 if len(line) >= 20:
-                    cnstnt = np.float(
+                    cnstnt = float(
                         line[10:20].strip().lower().replace("d", "e")
                     )
                 else:
                     cnstnt = 0.0
             else:
                 if len(line) >= 20:
-                    cnstnt = np.int(line[10:20].strip())
+                    cnstnt = int(line[10:20].strip())
                 else:
                     cnstnt = 0
                 # if cnstnt == 0:
@@ -3233,7 +3233,7 @@ def parse_control_record(
                 else:
                     fmtin = ""
                 try:
-                    iprn = np.int(line[40:50].strip())
+                    iprn = int(line[40:50].strip())
                 except:
                     iprn = 0
             # locat = int(raw[0])
@@ -3244,7 +3244,7 @@ def parse_control_record(
                 freefmt = "constant"
             elif locat < 0:
                 freefmt = "external"
-                nunit = np.int(locat) * -1
+                nunit = int(locat) * -1
                 fmtin = "(binary)"
             elif locat > 0:
                 # if the unit number matches the current file, it's internal
@@ -3252,7 +3252,7 @@ def parse_control_record(
                     freefmt = "internal"
                 else:
                     freefmt = "external"
-                nunit = np.int(locat)
+                nunit = int(locat)
 
             # Reset for special MT3D control flags
             if array_format == "mt3d":
diff --git a/flopy/utils/util_list.py b/flopy/utils/util_list.py
index 71ef5db33b..9e166e5e94 100644
--- a/flopy/utils/util_list.py
+++ b/flopy/utils/util_list.py
@@ -1143,12 +1143,12 @@ def to_array(self, kper=0, mask=False):
         for name, arr in arrays.items():
             if unstructured:
                 cnt = np.zeros(
-                    (self._model.nlay * self._model.ncpl,), dtype=np.float
+                    (self._model.nlay * self._model.ncpl,), dtype=float
                 )
             else:
                 cnt = np.zeros(
                     (self._model.nlay, self._model.nrow, self._model.ncol),
-                    dtype=np.float,
+                    dtype=float,
                 )
             # print(name,kper)
             for rec in sarr:
diff --git a/flopy/utils/voronoi.py b/flopy/utils/voronoi.py
index 9d87bcb9b7..eee2674835 100644
--- a/flopy/utils/voronoi.py
+++ b/flopy/utils/voronoi.py
@@ -18,7 +18,7 @@ def get_sorted_vertices(icell_vertices, vertices):
 
 
 def get_valid_faces(vor):
-    nvalid_faces = np.zeros(vor.npoints, dtype=np.int)
+    nvalid_faces = np.zeros(vor.npoints, dtype=int)
     for pointidx, simplex in zip(vor.ridge_points, vor.ridge_vertices):
         simplex = np.asarray(simplex)
         if np.all(simplex >= 0):