diff --git a/pyproject.toml b/pyproject.toml
index 0e9a5b4ef..ca9ff401b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -164,10 +164,11 @@ markers = [
 ignore_directives = ["argparse", "automodule"]
 
 [tool.ruff]
-lint.ignore = [
-    "C901",
-]
-lint.select = [
+src = ["src"]
+line-length = 88
+
+[tool.ruff.lint]
+select = [
     "C",
     "E",
     "F",
@@ -178,8 +179,12 @@ lint.select = [
     "RSE",
     "SIM",
     "W",
+    "PD",
+]
+ignore = [
+    "C901",
+    "PD013",
 ]
-line-length = 88
 
 [tool.ruff.lint.isort]
 combine-as-imports = true
diff --git a/src/subscript/check_swatinit/check_swatinit.py b/src/subscript/check_swatinit/check_swatinit.py
index 8bc055df2..1b3fe872a 100644
--- a/src/subscript/check_swatinit/check_swatinit.py
+++ b/src/subscript/check_swatinit/check_swatinit.py
@@ -80,7 +80,9 @@ def main() -> None:
         print(f"Dumping volume plot to {args.volplotfile}")
         pyplot.savefig(args.volplotfile)
 
-    if (args.plotfile or args.plot) and args.eqlnum not in qc_frame["EQLNUM"].values:
+    if (args.plotfile or args.plot) and args.eqlnum not in qc_frame[
+        "EQLNUM"
+    ].to_numpy():
         sys.exit(f"Error: EQLNUM {args.eqlnum} does not exist in grid. No plotting.")
     if args.plot or args.plotfile:
         plotter.plot_qc_panels(qc_frame[qc_frame["EQLNUM"] == args.eqlnum])
@@ -129,7 +131,7 @@ def reorder_dframe_for_nonnans(dframe: pd.DataFrame) -> pd.DataFrame:
     will aid data analysis application to deduce correct datatypes for
     columns"""
     null_count = "__NULL_COUNT__"
-    dframe[null_count] = dframe.isnull().sum(axis=1)
+    dframe[null_count] = dframe.isna().sum(axis=1)
     return (
         dframe.sort_values(null_count).drop(null_count, axis=1).reset_index(drop=True)
     )
@@ -253,11 +255,11 @@ def make_qc_gridframe(eclfiles: res2df.ResdataFiles) -> pd.DataFrame:
         # GLOBAL_INDEX is 0-indexed.
         grid_df["SWATINIT_DECK"] = pd.Series(swatinit_deckdata)[
             grid_df["GLOBAL_INDEX"].astype(int).tolist()
-        ].values
+        ].to_numpy()
 
     if "SWATINIT" not in grid_df:
         # OPM-flow does not include SWATINIT in the INIT file.
-        grid_df.rename({"SWATINIT_DECK": "SWATINIT"}, axis="columns", inplace=True)
+        grid_df = grid_df.rename({"SWATINIT_DECK": "SWATINIT"}, axis="columns")
     elif "SWATINIT_DECK" in grid_df:
         # (if SWATINIT is inputted using binary data in Eclipse deck, the code above
         # is not able to extract it)
@@ -352,7 +354,7 @@ def qc_flag(qc_frame: pd.DataFrame) -> pd.DataFrame:
         qc_col[
             (~np.isclose(qc_frame["OIP_INIT"], 0))
             & (~np.isclose(qc_frame["SWAT"], qc_frame["SWATINIT"], atol=1e-6))
-            & (~pd.isnull(qc_frame["PC_SCALING"]))
+            & (~pd.isna(qc_frame["PC_SCALING"]))
         ] = __FINE_EQUIL__
 
     # SWATINIT=1 above contact:
@@ -409,9 +411,7 @@ def qc_flag(qc_frame: pd.DataFrame) -> pd.DataFrame:
 
     # Tag the remainder with "unknown", when/if this happens, it is a bug or a
     # feature request:
-    qc_col.fillna(__UNKNOWN__, inplace=True)
-
-    return qc_col
+    return qc_col.fillna(__UNKNOWN__)
 
 
 def qc_volumes(qc_frame: pd.DataFrame) -> Dict[str, float]:
@@ -496,8 +496,10 @@ def _evaluate_pc(
             np.interp(
                 swat,
                 swl
-                + (satfunc[sat_name].values - sw_min) / (sw_max - sw_min) * (swu - swl),
-                satfunc[pc_name].values * pc_scaling,
+                + (satfunc[sat_name].to_numpy() - sw_min)
+                / (sw_max - sw_min)
+                * (swu - swl),
+                satfunc[pc_name].to_numpy() * pc_scaling,
             )
         )
     return p_cap
@@ -534,15 +536,15 @@ def compute_pc(qc_frame: pd.DataFrame, satfunc_df: pd.DataFrame) -> pd.Series:
 
     for satnum, satnum_frame in qc_frame.groupby("SATNUM"):
         if "SWLPC" in satnum_frame:
-            swls = satnum_frame["SWLPC"].values
+            swls = satnum_frame["SWLPC"].to_numpy()
         elif "SWL" in satnum_frame:
-            swls = satnum_frame["SWL"].values
+            swls = satnum_frame["SWL"].to_numpy()
         else:
             swls = None
-        swus = satnum_frame["SWU"].values if "SWU" in satnum_frame else None
+        swus = satnum_frame["SWU"].to_numpy() if "SWU" in satnum_frame else None
         p_cap[satnum_frame.index] = _evaluate_pc(
-            satnum_frame["SWAT"].values,
-            satnum_frame["PC_SCALING"].values,
+            satnum_frame["SWAT"].to_numpy(),
+            satnum_frame["PC_SCALING"].to_numpy(),
             swls,
             swus,
             satfunc_df[satfunc_df["SATNUM"] == satnum],
@@ -596,7 +598,7 @@ def merge_equil(grid_df: pd.DataFrame, equil_df: pd.DataFrame) -> pd.DataFrame:
     assert "PRESSURE" in equil_df
 
     # Be compatible with future change in res2df:
-    equil_df.rename({"ACCURACY": "OIP_INIT"}, axis="columns", inplace=True)
+    equil_df = equil_df.rename({"ACCURACY": "OIP_INIT"}, axis="columns")
 
     contacts = list({"OWC", "GOC", "GWC"}.intersection(set(equil_df.columns)))
     # Rename and slice the equil dataframe:
@@ -607,7 +609,7 @@ def merge_equil(grid_df: pd.DataFrame, equil_df: pd.DataFrame) -> pd.DataFrame:
         equil_df = equil_df[equil_df["KEYWORD"] == "EQUIL"]
     equil_df = equil_df[["Z_DATUM", "PRESSURE_DATUM", "EQLNUM", "OIP_INIT"] + contacts]
     equil_df["EQLNUM"] = equil_df["EQLNUM"].astype(int)
-    assert not pd.isnull(equil_df).any().any(), (
+    assert not pd.isna(equil_df).any().any(), (
         f"BUG: NaNs in equil dataframe:\n{equil_df}"
     )
     return grid_df.merge(equil_df, on="EQLNUM", how="left")
diff --git a/src/subscript/check_swatinit/plotter.py b/src/subscript/check_swatinit/plotter.py
index 0c27d60d8..69d0702e5 100644
--- a/src/subscript/check_swatinit/plotter.py
+++ b/src/subscript/check_swatinit/plotter.py
@@ -63,8 +63,8 @@ def plot_qc_panels(
     pyplot.subplot(2, 2, 4)
     pc_depth(qc_frame)
 
-    oip_init = qc_frame["OIP_INIT"].values[0]
-    eqlnum = qc_frame["EQLNUM"].values[0]
+    oip_init = qc_frame["OIP_INIT"].to_numpy()[0]
+    eqlnum = qc_frame["EQLNUM"].to_numpy()[0]
     pyplot.suptitle(f"EQLNUM: {eqlnum}, OIP_INIT: {oip_init}")
     if show:
         pyplot.show()
@@ -150,15 +150,15 @@ def pc_depth(
 def add_contacts_to_plot(qc_frame: pd.DataFrame, axis: pyplot.Axes) -> None:
     """Annotate axes with named horizontal lines for contacts."""
     if "OWC" in qc_frame:
-        owc = qc_frame["OWC"].values[0]  # OWC is assumed constant in the dataframe
+        owc = qc_frame["OWC"].to_numpy()[0]  # OWC is assumed constant in the dataframe
         axis.axhline(owc, color="black", linestyle="--", linewidth=1)
         axis.annotate(f"OWC={owc:g}", (0, owc))
     if "GOC" in qc_frame:
-        goc = qc_frame["GOC"].values[0]
+        goc = qc_frame["GOC"].to_numpy()[0]
         axis.axhline(goc, color="black", linestyle="--", linewidth=1)
         axis.annotate(f"GOC={goc:g}", (0, goc))
     if "GWC" in qc_frame:
-        gwc = qc_frame["GWC"].values[0]
+        gwc = qc_frame["GWC"].to_numpy()[0]
         axis.axhline(gwc, color="black", linestyle="--", linewidth=1)
         axis.annotate(f"GWC={gwc:g}", (0, gwc))
 
@@ -199,10 +199,10 @@ def wvol_waterfall(qc_vols: Dict[str, float]) -> None:
     blank.loc["SWAT_WVOL"] = 0
 
     fig = trans.plot(kind="bar", alpha=0.7, stacked=True, legend=None, bottom=blank)
-    fig.plot(step.index, step.values, "k")
+    fig.plot(step.index, step.to_numpy(), "k")
     pyplot.gcf().subplots_adjust(bottom=0.25)
 
-    blanktrans = blank.values + trans["volume"].values
+    blanktrans = blank.to_numpy() + trans["volume"].to_numpy()
     span = blank.max() - blanktrans[1:-1].min()
 
     if np.isclose(span, 0.0):
diff --git a/src/subscript/csv2ofmvol/csv2ofmvol.py b/src/subscript/csv2ofmvol/csv2ofmvol.py
index 054497ded..cc635650f 100644
--- a/src/subscript/csv2ofmvol/csv2ofmvol.py
+++ b/src/subscript/csv2ofmvol/csv2ofmvol.py
@@ -116,7 +116,7 @@ def read_pdm_csv_files(
     # Reindex:
     data = data.set_index(["WELL", "DATE"])
 
-    if not [data.columns.values]:
+    if not [data.columns.to_numpy()]:
         raise ValueError("No data columns found")
 
     # Drop duplicate multiindices (WELL, DATE)
@@ -220,7 +220,7 @@ def df2vol(data: pd.DataFrame) -> str:
 
     # Fill empty cells with zeros, empty cells can stem from concatenation
     # of dataframes with gas and water injectors.
-    voldata.fillna(value=0.0, inplace=True)
+    voldata = voldata.fillna(value=0.0)
 
     volstr = ""
     volstr += "*METRIC\n"
diff --git a/src/subscript/csv_merge/csv_merge.py b/src/subscript/csv_merge/csv_merge.py
index aa15b08e3..86eb827ae 100755
--- a/src/subscript/csv_merge/csv_merge.py
+++ b/src/subscript/csv_merge/csv_merge.py
@@ -307,7 +307,7 @@ def csv_merge_main(
             if len(merged_df[col].unique()) == 1:
                 columnstodelete.append(col)
         logger.info("Dropping constant columns %s", str(columnstodelete))
-        merged_df.drop(columnstodelete, inplace=True, axis=1)
+        merged_df = merged_df.drop(columnstodelete, axis=1)
 
     if merged_df.empty:
         logger.error("No data to output")
diff --git a/src/subscript/csv_stack/csv_stack.py b/src/subscript/csv_stack/csv_stack.py
index 1e42519d5..242016513 100755
--- a/src/subscript/csv_stack/csv_stack.py
+++ b/src/subscript/csv_stack/csv_stack.py
@@ -311,7 +311,7 @@ def csv_stack(
 
         # Stack the multiindex columns, this will add a lot of rows to
         # our ensemble, and condense the number of columns
-        dframe = dframe.stack()
+        dframe = dframe.stack(future_stack=True)
 
         # The values from non-multiindex-columns must be propagated to
         # the rows that emerged from the stacking. If you use the
diff --git a/src/subscript/fmuobs/parsers.py b/src/subscript/fmuobs/parsers.py
index 29527ead5..737302380 100644
--- a/src/subscript/fmuobs/parsers.py
+++ b/src/subscript/fmuobs/parsers.py
@@ -458,8 +458,7 @@ def resinsight_df2df(ri_dframe: pd.DataFrame) -> pd.DataFrame:
     if ri_dframe.empty:
         return pd.DataFrame()
 
-    dframe = ri_dframe.copy()
-    dframe.rename({"VECTOR": "KEY"}, axis="columns", inplace=True)
+    dframe = ri_dframe.copy().rename({"VECTOR": "KEY"}, axis="columns")
     dframe["LABEL"] = (
         dframe["KEY"].astype(str)
         + "-"
diff --git a/src/subscript/fmuobs/writers.py b/src/subscript/fmuobs/writers.py
index 321eee9f1..b136b97af 100644
--- a/src/subscript/fmuobs/writers.py
+++ b/src/subscript/fmuobs/writers.py
@@ -33,20 +33,20 @@ def dfsummary2ertobs(obs_df: pd.DataFrame) -> str:
     for _, row in smry_df.iterrows():
         ertobs_str += "SUMMARY_OBSERVATION " + str(row["LABEL"]) + "\n"
         ertobs_str += "{\n"
-        if "COMMENT" in row and not pd.isnull(row["COMMENT"]):
+        if "COMMENT" in row and not pd.isna(row["COMMENT"]):
             ertobs_str += (
                 "    -- "
                 + str(row["COMMENT"]).replace("\n", "\n    -- ").strip()
                 + "\n"
             )
-        if "DATE" in row and not pd.isnull(row["DATE"]):
+        if "DATE" in row and not pd.isna(row["DATE"]):
             ertobs_str += (
                 "    DATE = "
                 + str(pd.to_datetime(row["DATE"]).strftime(ERT_ISO_DATE_FORMAT))
                 + ";\n"
             )
         for dataname in ["KEY", "DAYS", "RESTART", "VALUE", "ERROR", "SOURCE"]:
-            if dataname in row and not pd.isnull(row[dataname]):
+            if dataname in row and not pd.isna(row[dataname]):
                 ertobs_str += "    " + dataname + " = " + str(row[dataname]) + ";\n"
         ertobs_str += "};\n"
     return ertobs_str
@@ -70,12 +70,14 @@ def dfblock2ertobs(obs_df: pd.DataFrame) -> str:
         )
     for obslabel, block_df in block_obs_df.groupby("LABEL"):
         ertobs_str += "BLOCK_OBSERVATION " + obslabel + "\n{\n"
-        if "COMMENT" in block_df and not pd.isnull(block_df["COMMENT"]).any():
+        if "COMMENT" in block_df and not pd.isna(block_df["COMMENT"]).any():
             if len(block_df["COMMENT"].dropna().unique()) != 1:
                 logger.warning("Inconsistency in COMMENT in block dataframe")
             ertobs_str += (
                 "    -- "
-                + str(block_df["COMMENT"].values[0]).replace("\n", "\n    -- ").strip()
+                + str(block_df["COMMENT"].to_numpy()[0])
+                .replace("\n", "\n    -- ")
+                .strip()
                 + "\n"
             )
         for dataname in ["FIELD", "DATE"]:
@@ -85,7 +87,7 @@ def dfblock2ertobs(obs_df: pd.DataFrame) -> str:
                         "    "
                         + dataname
                         + " = "
-                        + str(block_df[dataname].values[0])
+                        + str(block_df[dataname].to_numpy()[0])
                         + ";\n"
                     )
                 else:
@@ -94,7 +96,7 @@ def dfblock2ertobs(obs_df: pd.DataFrame) -> str:
                         f"block dataframe for one label has multiple {dataname}"
                     )
         for _, row in block_df.iterrows():
-            if "SUBCOMMENT" in row and not pd.isnull("SUBCOMMENT"):
+            if "SUBCOMMENT" in row and not pd.isna("SUBCOMMENT"):
                 ertobs_str += (
                     "    -- "
                     + str(row["SUBCOMMENT"]).strip().replace("\n", "\n    -- ").strip()
@@ -102,7 +104,7 @@ def dfblock2ertobs(obs_df: pd.DataFrame) -> str:
                 )
             ertobs_str += "    OBS " + row["OBS"] + " {"
             for dataname in ["I", "J", "K", "VALUE", "ERROR", "SOURCE"]:
-                if dataname in row and not pd.isnull(row[dataname]):
+                if dataname in row and not pd.isna(row[dataname]):
                     ertobs_str += " " + dataname + " = " + str(row[dataname]) + ";"
             ertobs_str += "};\n"
         ertobs_str += "};\n"
@@ -138,17 +140,17 @@ def dfhistory2ertobs(obs_df: pd.DataFrame) -> str:
             .to_dict(orient="records")[0]
         )
         for dataname in ["ERROR", "ERROR_MODE", "ERROR_MIN"]:
-            if dataname in default_row and not pd.isnull(default_row[dataname]):
+            if dataname in default_row and not pd.isna(default_row[dataname]):
                 ertobs_str += (
                     "  " + dataname + " = " + str(default_row[dataname]) + ";\n"
                 )
         for _, row in history_df.iterrows():
-            if "SEGMENT" in row and not pd.isnull(row["SEGMENT"]):
+            if "SEGMENT" in row and not pd.isna(row["SEGMENT"]):
                 if row["SEGMENT"] == "DEFAULT":
                     continue
                 ertobs_str += "  SEGMENT " + row["SEGMENT"] + " {"
                 for dataname in ["START", "STOP", "ERROR", "ERROR_MODE"]:
-                    if dataname in row and not pd.isnull(row[dataname]):
+                    if dataname in row and not pd.isna(row[dataname]):
                         ertobs_str += " " + dataname + " = " + str(row[dataname]) + ";"
                 ertobs_str += "};\n"
         ertobs_str += "};\n"
@@ -185,7 +187,7 @@ def dfgeneral2ertobs(obs_df: pd.DataFrame) -> str:
             "INDEX_LIST",
             "ERROR_COVAR",
         ]:
-            if dataname in row and not pd.isnull(row[dataname]):
+            if dataname in row and not pd.isna(row[dataname]):
                 ertobs_str += "    " + dataname + " = " + str(row[dataname]) + ";\n"
         ertobs_str += "};\n"
 
@@ -240,11 +242,11 @@ def summary_df2obsdict(smry_df: pd.DataFrame) -> List[dict]:
     assert isinstance(smry_df, pd.DataFrame)
     if "CLASS" in smry_df:
         assert len(smry_df["CLASS"].unique()) == 1
-        smry_df.drop("CLASS", axis=1, inplace=True)
+        smry_df = smry_df.drop("CLASS", axis=1)
 
     smry_obs_list = []
     if isinstance(smry_df, pd.DataFrame):
-        smry_df.dropna(axis=1, how="all", inplace=True)
+        smry_df = smry_df.dropna(axis=1, how="all")
 
     if "DATE" not in smry_df:
         raise ValueError("Can't have summary observation without a date")
@@ -258,10 +260,10 @@ def summary_df2obsdict(smry_df: pd.DataFrame) -> List[dict]:
     for smrykey, smrykey_df in smry_df.groupby("KEY"):
         smry_obs_element = {}
         smry_obs_element["key"] = smrykey
-        if "COMMENT" in smrykey_df and not pd.isnull(smrykey_df["COMMENT"]).all():
+        if "COMMENT" in smrykey_df and not pd.isna(smrykey_df["COMMENT"]).all():
             smry_obs_element["comment"] = smrykey_df["COMMENT"].unique()[0]
         if isinstance(smrykey_df, pd.DataFrame):
-            smrykey_df.drop("KEY", axis=1, inplace=True)
+            smrykey_df = smrykey_df.drop("KEY", axis=1)
         if "SUBCOMMENT" in smrykey_df:
             smrykey_df["COMMENT"] = smrykey_df["SUBCOMMENT"]
             del smrykey_df["SUBCOMMENT"]
@@ -317,14 +319,12 @@ def block_df2obsdict(block_df: pd.DataFrame) -> List[dict]:
     block_obs_list = []
     if "CLASS" in block_df:
         assert len(block_df["CLASS"].unique()) == 1
-        block_df.drop("CLASS", axis=1, inplace=True)
+        block_df = block_df.drop("CLASS", axis=1)
 
     if "DATE" not in block_df:
         raise ValueError("Can't have rft/block observation without a date")
 
-    block_df = convert_dframe_date_to_str(block_df)
-
-    block_df.dropna(axis=1, how="all", inplace=True)
+    block_df = convert_dframe_date_to_str(block_df).dropna(axis=1, how="all")
 
     for blocklabel, blocklabel_df in block_df.groupby(["LABEL", "DATE"]):
         blocklabel_dict = {}
@@ -372,13 +372,13 @@ def df2obsdict(obs_df: pd.DataFrame) -> dict:
         return {}
 
     # Process SUMMARY_OBSERVATION:
-    if "SUMMARY_OBSERVATION" in obs_df["CLASS"].values:
+    if "SUMMARY_OBSERVATION" in obs_df["CLASS"].to_numpy():
         obsdict[CLASS_SHORTNAME["SUMMARY_OBSERVATION"]] = summary_df2obsdict(
             obs_df.set_index("CLASS").loc[["SUMMARY_OBSERVATION"]]
         )
 
     # Process BLOCK_OBSERVATION:
-    if "BLOCK_OBSERVATION" in obs_df["CLASS"].values:
+    if "BLOCK_OBSERVATION" in obs_df["CLASS"].to_numpy():
         obsdict[CLASS_SHORTNAME["BLOCK_OBSERVATION"]] = block_df2obsdict(
             obs_df.set_index("CLASS").loc[["BLOCK_OBSERVATION"]]
         )
@@ -408,9 +408,9 @@ def df2resinsight_df(obs_df: pd.DataFrame) -> pd.DataFrame:
     ri_dframe = obs_df.copy()
 
     # Only SUMMARY_OBSERVATION is supported:
-    ri_dframe = ri_dframe[ri_dframe["CLASS"] == "SUMMARY_OBSERVATION"]
-
-    ri_dframe.rename({"KEY": "VECTOR"}, axis="columns", inplace=True)
+    ri_dframe = ri_dframe[ri_dframe["CLASS"] == "SUMMARY_OBSERVATION"].rename(
+        {"KEY": "VECTOR"}, axis="columns"
+    )
 
     # Ensure all vectors are present:
     for ri_vec in ri_column_names:
diff --git a/src/subscript/interp_relperm/interp_relperm.py b/src/subscript/interp_relperm/interp_relperm.py
index 66bb0df6b..46649376c 100755
--- a/src/subscript/interp_relperm/interp_relperm.py
+++ b/src/subscript/interp_relperm/interp_relperm.py
@@ -169,7 +169,7 @@ def make_wateroilgas(dframe: pd.DataFrame, delta_s: float) -> pyscal.WaterOilGas
 
     # If we have parsed family 2 input, KRO and KROW are not
     # on the same row. Merge the rows into family 1 style:
-    if "KEYWORD" in dframe and "SOF3" in dframe["KEYWORD"].values:
+    if "KEYWORD" in dframe and "SOF3" in dframe["KEYWORD"].to_numpy():
         sof3_rows = dframe["KEYWORD"] == "SOF3"
         dframe.loc[sof3_rows, "SW"] = 1 - dframe[sof3_rows]["SO"]
         swl = dframe["SW"].min()
@@ -406,9 +406,9 @@ def process_config(cfg: Dict[str, Any], root_path: Optional[Path] = None) -> Non
         set(base_df.columns) == set(low_df.columns)
         and set(base_df.columns) == set(high_df.columns)
     ):
-        logger.error("Base input had columns: %s", str(base_df.columns.values))
-        logger.error("Low input had columns: %s", str(low_df.columns.values))
-        logger.error("High input had columns: %s", str(high_df.columns.values))
+        logger.error("Base input had columns: %s", str(base_df.columns.to_numpy()))
+        logger.error("Low input had columns: %s", str(low_df.columns.to_numpy()))
+        logger.error("High input had columns: %s", str(high_df.columns.to_numpy()))
         logger.error("Inconsistent input data, check keywords in input files")
         sys.exit(1)
 
diff --git a/src/subscript/merge_rft_ertobs/merge_rft_ertobs.py b/src/subscript/merge_rft_ertobs/merge_rft_ertobs.py
index bf8164555..660a8f91c 100644
--- a/src/subscript/merge_rft_ertobs/merge_rft_ertobs.py
+++ b/src/subscript/merge_rft_ertobs/merge_rft_ertobs.py
@@ -237,9 +237,9 @@ def merge_rft_ertobs(gendatacsv: str, obsdir: str) -> pd.DataFrame:
         logger.info("Found %s active observation points", str(len(obs_df)))
 
     if "report_step" in sim_df.columns:
-        return pd.merge(sim_df, obs_df, how="left", on=["well", "order", "report_step"])
+        return sim_df.merge(obs_df, how="left", on=["well", "order", "report_step"])
     # Ensure backward compatibility where gendata_rft doesn't have report_step
-    return pd.merge(sim_df, obs_df, how="left", on=["well", "order"])
+    return sim_df.merge(obs_df, how="left", on=["well", "order"])
 
 
 def main() -> None:
diff --git a/src/subscript/params2csv/params2csv.py b/src/subscript/params2csv/params2csv.py
index 6661045ba..eee8486ff 100755
--- a/src/subscript/params2csv/params2csv.py
+++ b/src/subscript/params2csv/params2csv.py
@@ -176,10 +176,9 @@ def params2csv_main(args: argparse.Namespace) -> None:
             header=None,
             usecols=[0, 1],
             sep=r"\s+",
-        )
-
-        paramtable.drop_duplicates(
-            "key", keep="last", inplace=True
+        ).drop_duplicates(
+            "key",
+            keep="last",
         )  # if key is repeated, keep the last one.
         transposed = paramtable.set_index("key").transpose()
         if args.filenamecolumnname in transposed.columns:
diff --git a/src/subscript/presentvalue/presentvalue.py b/src/subscript/presentvalue/presentvalue.py
index b46f0cf0d..fe708f415 100755
--- a/src/subscript/presentvalue/presentvalue.py
+++ b/src/subscript/presentvalue/presentvalue.py
@@ -373,7 +373,7 @@ def calc_presentvalue_df(
     prodecon[["oilprice", "gasprice", "usdtonok", "discountrate"]] = prodecon[
         ["oilprice", "gasprice", "usdtonok", "discountrate"]
     ].bfill()
-    prodecon.fillna(value=0, inplace=True)  # Zero-pad other data (costs)
+    prodecon = prodecon.fillna(value=0)  # Zero-pad other data (costs)
 
     prodecon["deltayears"] = prodecon.index - discountto
 
@@ -422,15 +422,18 @@ def get_yearly_summary(
     ):
         raise ValueError("Only cumulative Eclipse vectors can be used")
     eclfiles = res2df.ResdataFiles(eclfile)
-    sum_df = res2df.summary.df(
-        eclfiles, column_keys=[oilvector, gasvector, gasinjvector], time_index="yearly"
-    )
-    sum_df.rename(
-        {oilvector: "OPT", gasvector: "GPT", gasinjvector: "GIT"},
-        axis="columns",
-        inplace=True,
+    sum_df = (
+        res2df.summary.df(
+            eclfiles,
+            column_keys=[oilvector, gasvector, gasinjvector],
+            time_index="yearly",
+        )
+        .rename(
+            {oilvector: "OPT", gasvector: "GPT", gasinjvector: "GIT"},
+            axis="columns",
+        )
+        .reset_index()
     )
-    sum_df = sum_df.reset_index()
 
     if "GIT" not in sum_df:
         sum_df["GIT"] = 0
diff --git a/src/subscript/prtvol2csv/prtvol2csv.py b/src/subscript/prtvol2csv/prtvol2csv.py
index 24d19b7ad..148eaa0e6 100644
--- a/src/subscript/prtvol2csv/prtvol2csv.py
+++ b/src/subscript/prtvol2csv/prtvol2csv.py
@@ -193,10 +193,9 @@ def currently_in_place_from_prt(
     inplace_df = inplace_df[inplace_df["DATATYPE"] == "CURRENTLY IN PLACE"]
 
     # Cleanup dataframe:
-    inplace_df.drop(
-        ["DATATYPE", "TO_REGION", "FIPNAME", "DATE"], axis="columns", inplace=True
-    )
-    inplace_df.set_index("REGION", inplace=True)
+    inplace_df = inplace_df.drop(
+        ["DATATYPE", "TO_REGION", "FIPNAME", "DATE"], axis="columns"
+    ).set_index("REGION")
     inplace_df.index.name = fipname  # Use "FIPNUM" if not handled by Webviz
 
     logger.info("Extracted CURRENTLY IN PLACE from %s at date %s", prt_file, date_str)
diff --git a/src/subscript/restartthinner/restartthinner.py b/src/subscript/restartthinner/restartthinner.py
index e0f27428a..4cdfd4fe9 100644
--- a/src/subscript/restartthinner/restartthinner.py
+++ b/src/subscript/restartthinner/restartthinner.py
@@ -145,7 +145,7 @@ def restartthinner(
                 pandas.Timestamp(restart_dates[-1]).value,
                 int(numberofslices),
             )
-        ).values
+        ).to_numpy()
     else:
         slicedates = [restart_dates[-1]]  # Only return last date if only one is wanted
 
diff --git a/src/subscript/vfp2csv/vfp2csv.py b/src/subscript/vfp2csv/vfp2csv.py
index 2d8062c67..0b5f5b73a 100644
--- a/src/subscript/vfp2csv/vfp2csv.py
+++ b/src/subscript/vfp2csv/vfp2csv.py
@@ -171,16 +171,16 @@ def vfpfile2df(filename: str) -> pd.DataFrame:
         bhp_values.columns = pd.MultiIndex.from_tuples(indextuples)
 
         # Now stack
-        bhp_values_stacked = bhp_values.stack()
+        bhp_values_stacked = bhp_values.stack(future_stack=True)
 
         # In order to propagate the gfr, thp, wct values after
         # stacking to the correct rows, we should either understand
         # how to do that properly using pandas, but for now, we try a
         # backwards fill, hopefully that is robust enough
-        bhp_values_stacked.bfill(inplace=True)
-        # Also reset the index:
-        bhp_values_stacked.reset_index(inplace=True)
-        bhp_values_stacked.drop("level_0", axis="columns", inplace=True)
+        # Also reset the index
+        bhp_values_stacked = (
+            bhp_values_stacked.bfill().reset_index().drop("level_0", axis="columns")
+        )
         # This column is not meaningful (it is the old index)
 
         # Delete rows that does not belong to any flow rate (this is
diff --git a/tests/test_check_swatinit.py b/tests/test_check_swatinit.py
index c56ca52c6..9e7f4ed4f 100644
--- a/tests/test_check_swatinit.py
+++ b/tests/test_check_swatinit.py
@@ -501,12 +501,12 @@ def test_compute_pc(propslist, satfunc_df, expected_pc):
     pc_series = compute_pc(qc_frame, satfunc_df)
     if qc_frame.empty:
         if not pc_series.empty:
-            assert all(pd.isnull(pc_series))
+            assert all(pd.isna(pc_series))
     else:
-        if pd.isnull(expected_pc):
-            assert pd.isnull(pc_series.values[0])
+        if pd.isna(expected_pc):
+            assert pd.isna(pc_series.to_numpy()[0])
         else:
-            assert pc_series.values[0] == expected_pc
+            assert pc_series.to_numpy()[0] == expected_pc
 
 
 def test_eqlnum2(tmp_path, mocker):
diff --git a/tests/test_check_swatinit_simulators.py b/tests/test_check_swatinit_simulators.py
index 893ae87c2..4c42653a4 100644
--- a/tests/test_check_swatinit_simulators.py
+++ b/tests/test_check_swatinit_simulators.py
@@ -90,8 +90,8 @@ def test_swat_higher_than_swatinit_via_swl_above_contact(simulator, tmp_path):
 
     # When SWL is truncated, we cannot trust PC_SCALING to be used to
     # compute PC, so it is removed from the dataframe.
-    assert pd.isnull(qc_frame["PC_SCALING"][0])
-    assert pd.isnull(qc_frame["PC"][0])
+    assert pd.isna(qc_frame["PC_SCALING"][0])
+    assert pd.isna(qc_frame["PC"][0])
 
 
 def test_swat_limited_by_ppcwmax_above_contact(simulator, tmp_path):
@@ -446,11 +446,11 @@ def test_swatinit_less_than_1_below_contact(simulator, tmp_path):
         # E100 will not report a PPCW in this case, resdata gives -1e20,
         # which becomes a NaN through res2df and then NaN columns are dropped.
         if "PPCW" in qc_frame:
-            assert pd.isnull(qc_frame["PPCW"][0])
+            assert pd.isna(qc_frame["PPCW"][0])
         if "PC_SCALING" in qc_frame:
-            assert pd.isnull(qc_frame["PC_SCALING"][0])
+            assert pd.isna(qc_frame["PC_SCALING"][0])
         if "PC" in qc_frame:
-            assert pd.isnull(qc_frame["PC"][0])
+            assert pd.isna(qc_frame["PC"][0])
 
 
 @pytest.mark.skipif(IN_SUBSCRIPT_GITHUB_ACTIONS, reason="Test require flow dev version")
@@ -502,11 +502,11 @@ def test_swatinit_less_than_1_below_contact_neg_pc(simulator, tmp_path):
         assert np.isclose(qc_frame["SWAT"][0], expected_swat)
         # PPCW is set to NaN, so we don't have that column
         if "PPCW" in qc_frame:
-            assert pd.isnull(qc_frame["PPCW"][0])
+            assert pd.isna(qc_frame["PPCW"][0])
         if "PC_SCALING" in qc_frame:
-            assert pd.isnull(qc_frame["PC_SCALING"][0])
+            assert pd.isna(qc_frame["PC_SCALING"][0])
         if "PC" in qc_frame:
-            assert pd.isnull(qc_frame["PC"][0])
+            assert pd.isna(qc_frame["PC"][0])
         assert np.isclose(qc_frame["PCW"][0], 3.0)  # Untouched input
 
         assert np.isclose(
@@ -631,7 +631,7 @@ def test_swatinit_1_below_contact(simulator, tmp_path):
         assert np.isclose(qc_frame["PC"][0], 0)
     else:
         if "PPCW" in qc_frame:
-            assert pd.isnull(qc_frame["PPCW"][0])
+            assert pd.isna(qc_frame["PPCW"][0])
 
     qc_vols = qc_volumes(qc_frame)
     assert np.isclose(qc_vols[__WATER__], 0.0)
diff --git a/tests/test_csv_stack.py b/tests/test_csv_stack.py
index 609de8916..ed102fce6 100644
--- a/tests/test_csv_stack.py
+++ b/tests/test_csv_stack.py
@@ -157,8 +157,8 @@ def test_commandlinetool(tmp_path, mocker):
     assert "WELL" in stacked
     assert "WOPT:A1" not in stacked
     assert "WOPT" in stacked
-    assert "A1" in stacked["WELL"].values
-    assert "A2" in stacked["WELL"].values
+    assert "A1" in stacked["WELL"].to_numpy()
+    assert "A2" in stacked["WELL"].to_numpy()
     assert "CONST" not in stacked
 
     mocker.patch(
@@ -197,8 +197,8 @@ def test_commandlinetool(tmp_path, mocker):
     assert "REGION" in stacked
     assert "CONST" not in stacked
     assert "RPR" in stacked
-    assert 1 in stacked["REGION"].astype(int).values
-    assert 2 in stacked["REGION"].astype(int).values
+    assert 1 in stacked["REGION"].astype(int).to_numpy()
+    assert 2 in stacked["REGION"].astype(int).to_numpy()
 
 
 @pytest.mark.parametrize("verbose", [False, True])
diff --git a/tests/test_fmuobs.py b/tests/test_fmuobs.py
index 6f4ffaa77..0b760eb8c 100644
--- a/tests/test_fmuobs.py
+++ b/tests/test_fmuobs.py
@@ -112,9 +112,8 @@ def test_roundtrip_ertobs(filename, readonly_testdata_dir):
 
     # Convert to ERT obs format and back again:
     ertobs_str = df2ertobs(dframe)
-    ert_roundtrip_dframe = ertobs2df(ertobs_str)
-    ert_roundtrip_dframe.set_index("CLASS", inplace=True)
-    dframe.set_index("CLASS", inplace=True)
+    ert_roundtrip_dframe = ertobs2df(ertobs_str).set_index("CLASS")
+    dframe = dframe.set_index("CLASS")
 
     # This big loop is only here to aid in debugging when
     # the dataframes do not match, asserting equivalence of
@@ -126,25 +125,24 @@ def test_roundtrip_ertobs(filename, readonly_testdata_dir):
             .sort_index(axis=1)
         )
         subframe = dframe.loc[[_class]].dropna(axis=1, how="all").sort_index(axis=1)
-        roundtrip_subframe.set_index(
+        roundtrip_subframe = roundtrip_subframe.set_index(
             list(
                 {"CLASS", "LABEL", "OBS", "SEGMENT"}.intersection(
                     set(roundtrip_subframe.columns)
                 )
             ),
-            inplace=True,
-        )
-        roundtrip_subframe.sort_index(inplace=True)
-        subframe.set_index(
-            list(
-                {"CLASS", "LABEL", "OBS", "SEGMENT"}.intersection(set(subframe.columns))
-            ),
-            inplace=True,
-        )
-        subframe.sort_index(inplace=True)
+        ).sort_index()
         # Comments are not preservable through ertobs roundtrips:
-        subframe.drop(
-            ["COMMENT", "SUBCOMMENT"], axis="columns", errors="ignore", inplace=True
+        subframe = (
+            subframe.set_index(
+                list(
+                    {"CLASS", "LABEL", "OBS", "SEGMENT"}.intersection(
+                        set(subframe.columns)
+                    )
+                ),
+            )
+            .sort_index()
+            .drop(["COMMENT", "SUBCOMMENT"], axis="columns", errors="ignore")
         )
         if _class == "BLOCK_OBSERVATION" and "WELL" in subframe:
             # WELL as used in yaml is not preservable in roundtrips
@@ -189,9 +187,8 @@ def test_roundtrip_yaml(filename, readonly_testdata_dir):
     ].dropna(axis="columns", how="all")
     # Convert to YAML (really dict) format and back again:
     obsdict = df2obsdict(dframe)
-    yaml_roundtrip_dframe = obsdict2df(obsdict)
-    yaml_roundtrip_dframe.set_index("CLASS", inplace=True)
-    dframe.set_index("CLASS", inplace=True)
+    yaml_roundtrip_dframe = obsdict2df(obsdict).set_index("CLASS")
+    dframe = dframe.set_index("CLASS")
     if "WELL" in yaml_roundtrip_dframe:
         # WELL as used in yaml is not preservable in roundtrips
         del yaml_roundtrip_dframe["WELL"]
diff --git a/tests/test_merge_rft_ertobs.py b/tests/test_merge_rft_ertobs.py
index 0fb0123b0..7deadd1bc 100644
--- a/tests/test_merge_rft_ertobs.py
+++ b/tests/test_merge_rft_ertobs.py
@@ -125,7 +125,7 @@ def test_merge_drogon_inactive(drogondata):
     dframe = merge_rft_ertobs("gendata_rft.csv", "rft")
     assert not dframe.empty
     assert {"pressure", "observed", "error", "well", "time"}.issubset(dframe.columns)
-    assert sum(dframe["pressure"].isnull()) == 1
+    assert sum(dframe["pressure"].isna()) == 1
     assert not np.isclose(
         (dframe["observed"] - dframe["pressure"]).abs().mean(), 6.2141156
     )
@@ -138,7 +138,7 @@ def test_merge_drogon_missing_observation(drogondata):
     dframe = merge_rft_ertobs("gendata_rft.csv", "rft")
     assert not dframe.empty
     assert {"pressure", "observed", "error", "well", "time"}.issubset(dframe.columns)
-    assert sum(dframe["observed"].isnull()) == 1
+    assert sum(dframe["observed"].isna()) == 1
     assert not np.isclose(
         (dframe["observed"] - dframe["pressure"]).abs().mean(), 6.2141156
     )
diff --git a/tests/test_ofmvol2csv.py b/tests/test_ofmvol2csv.py
index b4f291e05..2fb7a899d 100644
--- a/tests/test_ofmvol2csv.py
+++ b/tests/test_ofmvol2csv.py
@@ -267,7 +267,7 @@ def test_parse_well(inputlines, expected):
     """Test parsing well data"""
     if "DATE" in expected:
         expected["DATE"] = pd.to_datetime(expected["DATE"])
-        expected.set_index(["WELL", "DATE"], inplace=True)
+        expected = expected.set_index(["WELL", "DATE"])
     # Assume there is DATE line in the test input
     inputlines = ofmvol2csv.cleanse_ofm_lines(inputlines)
     colnames = ofmvol2csv.extract_columnnames(inputlines)
@@ -362,7 +362,7 @@ def test_process_volstr(inputlines, expected):
     and well data"""
     if "DATE" in expected:
         expected["DATE"] = pd.to_datetime(expected["DATE"])
-        expected.set_index(["WELL", "DATE"], inplace=True)
+        expected = expected.set_index(["WELL", "DATE"])
     dframe = ofmvol2csv.process_volstr("\n".join(inputlines))
     pd.testing.assert_frame_equal(dframe, expected)
 
diff --git a/tests/test_params2csv.py b/tests/test_params2csv.py
index 45bc351f5..90ec05026 100644
--- a/tests/test_params2csv.py
+++ b/tests/test_params2csv.py
@@ -36,7 +36,7 @@ def test_main(tmp_path, mocker):
     assert "CONSTANT" not in result
     assert "BOGUS" not in result
     assert "filename" in result
-    assert set(result["filename"].values) == {"parameters1.txt", "parameters2.txt"}
+    assert set(result["filename"].to_numpy()) == {"parameters1.txt", "parameters2.txt"}
 
     # Test the cleaning mode:
     mocker.patch(
@@ -69,7 +69,7 @@ def test_main(tmp_path, mocker):
     assert "CONSTANT" not in result
     assert "BOGUS" not in result
     assert "filename" in result
-    assert set(result["filename"].values) == {"parameters1.txt", "parameters2.txt"}
+    assert set(result["filename"].to_numpy()) == {"parameters1.txt", "parameters2.txt"}
 
 
 def test_spaces_in_values(tmp_path, mocker):
@@ -84,7 +84,7 @@ def test_spaces_in_values(tmp_path, mocker):
     params2csv.main()
     result = pd.read_csv("params.csv")
     assert "somekey" in result
-    assert result["somekey"].values[0] == "value with spaces"
+    assert result["somekey"].to_numpy()[0] == "value with spaces"
 
 
 def test_spaces_in_values_single_quotes(tmp_path, mocker):
@@ -96,7 +96,7 @@ def test_spaces_in_values_single_quotes(tmp_path, mocker):
     params2csv.main()
     result = pd.read_csv("params.csv")
     assert "somekey" in result
-    assert result["somekey"].values[0] == "value with spaces"
+    assert result["somekey"].to_numpy()[0] == "value with spaces"
 
 
 @pytest.mark.integration
diff --git a/tests/test_presentvalue.py b/tests/test_presentvalue.py
index 541cd1366..b09b69d77 100644
--- a/tests/test_presentvalue.py
+++ b/tests/test_presentvalue.py
@@ -99,7 +99,7 @@ def test_prepare_econ_table_csv(tmp_path):
     assert econ_df["discountrate"].unique() == [8]  # defaulted
     assert econ_df["usdtonok"].unique() == [7]
     assert econ_df["costs"].unique() == [100]
-    assert econ_df.index.values == [2030]
+    assert econ_df.index.to_numpy() == [2030]
 
 
 ECONCOLS = ["year", "oilprice", "gasprice", "usdtonok", "costs", "discountrate"]
diff --git a/tests/test_prtvol2csv.py b/tests/test_prtvol2csv.py
index 5b424a23d..024fbb37b 100644
--- a/tests/test_prtvol2csv.py
+++ b/tests/test_prtvol2csv.py
@@ -644,20 +644,20 @@ def test_prtvol2df(tmp_path):
 
     assert prtvol2csv.prtvol2df(
         simv, resv, FipMapper(mapdata={"region2fipnum": {"West": 1}})
-    )["REGION"].values == ["West"]
+    )["REGION"].to_numpy() == ["West"]
 
     # Reverse the supplied map, should give the same:
     assert prtvol2csv.prtvol2df(
         simv, resv, FipMapper(mapdata={"fipnum2region": {1: "West"}})
-    )["REGION"].values == ["West"]
+    )["REGION"].to_numpy() == ["West"]
 
     # And then for zones:
     assert prtvol2csv.prtvol2df(
         simv, resv, FipMapper(mapdata={"fipnum2zone": {1: "Upper"}})
-    )["ZONE"].values == ["Upper"]
+    )["ZONE"].to_numpy() == ["Upper"]
     assert prtvol2csv.prtvol2df(
         simv, resv, FipMapper(mapdata={"zone2fipnum": {"Upper": 1}})
-    )["ZONE"].values == ["Upper"]
+    )["ZONE"].to_numpy() == ["Upper"]
     # if we use {"Upper": "1"} it will fail, but no pytest.raises on
     # that yet, perhaps it will be fixed later.
 
@@ -668,7 +668,7 @@ def test_prtvol2df(tmp_path):
     )
     assert prtvol2csv.prtvol2df(simv, resv, FipMapper(yamlfile="z2f_int.yml"))[
         "ZONE"
-    ].values == ["Upper"]
+    ].to_numpy() == ["Upper"]
     prtvol2csv.prtvol2df(simv, resv, FipMapper(yamlfile="z2f_int.yml")).to_csv(
         "foo.csv"
     )
@@ -679,8 +679,8 @@ def test_prtvol2df(tmp_path):
         resv,
         FipMapper(mapdata={"fipnum2region": {1: "West"}, "zone2fipnum": {"Upper": 1}}),
     )
-    assert volumes["REGION"].values == ["West"]
-    assert volumes["ZONE"].values == ["Upper"]
+    assert volumes["REGION"].to_numpy() == ["West"]
+    assert volumes["ZONE"].to_numpy() == ["Upper"]
 
     # fipnummaps referring to non-existing fipnums:
     volumes = prtvol2csv.prtvol2df(
@@ -702,7 +702,7 @@ def test_prtvol2df(tmp_path):
     )
     assert prtvol2csv.prtvol2df(
         simv, resv, FipMapper(yamlfile="global_master_config.yml")
-    )["ZONE"].values == ["Upper"]
+    )["ZONE"].to_numpy() == ["Upper"]
 
 
 def test_webviz_regiontofipnum_format(tmp_path):
@@ -717,8 +717,8 @@ def test_webviz_regiontofipnum_format(tmp_path):
         encoding="utf8",
     )
     dframe = prtvol2csv.prtvol2df(simv, resv, FipMapper(yamlfile="webviz_fip.yml"))
-    assert dframe["ZONE"].values == ["Volon"]
-    assert dframe["REGION"].values == ["West"]
+    assert dframe["ZONE"].to_numpy() == ["Volon"]
+    assert dframe["REGION"].to_numpy() == ["West"]
 
 
 @pytest.mark.integration
@@ -819,8 +819,8 @@ def test_prtvol2csv_regions_typemix(tmp_path, mocker):
     assert not dframe.empty
     assert "REGION" in dframe
     assert "ZONE" not in dframe
-    assert "RegionA" in dframe["REGION"].values
-    assert "8" in dframe["REGION"].values
+    assert "RegionA" in dframe["REGION"].to_numpy()
+    assert "8" in dframe["REGION"].to_numpy()
     assert len(dframe) == 6