diff --git a/setup.py b/setup.py
index debf4412..29aa9a3d 100644
--- a/setup.py
+++ b/setup.py
@@ -13,7 +13,7 @@
 REQUIREMENTS = [
     "resdata>=4.0.0",
     "numpy",
-    "pandas<2",
+    "pandas",
     "pyyaml>=5.1",
 ]
 
diff --git a/src/fmu/ensemble/observations.py b/src/fmu/ensemble/observations.py
index 61734d1d..e711a3a2 100644
--- a/src/fmu/ensemble/observations.py
+++ b/src/fmu/ensemble/observations.py
@@ -194,7 +194,7 @@ def load_smry(self, realization, smryvector, time_index="yearly", smryerror=None
             realization
         )
         virtobs["observations"] = []
-        for date, value in dataseries.iteritems():
+        for date, value in dataseries.items():
             virtobs["observations"].append(
                 {"value": value, "error": smryerror, "date": date}
             )
diff --git a/src/fmu/ensemble/realization.py b/src/fmu/ensemble/realization.py
index d05cae4e..95c1f21c 100644
--- a/src/fmu/ensemble/realization.py
+++ b/src/fmu/ensemble/realization.py
@@ -512,14 +512,7 @@ def load_status(self):
             on_bad_lines="skip",
         )
 
-        # dtype str messes up a little bit, pre-Pandas 0.24.1 gives 'None' as
-        # a string where data is missing.
-        status.replace("None", "", inplace=True)
-        # While Pandas 0.24.1 will insert proper Null values in those cells,
-        # we fill them with the empty string for the rest of this code to work
         status.fillna("", inplace=True)
-        # It should be ok to have both of these statements running, but the
-        # replace() is probably superfluous when pandas 0.23 is gone.
 
         errorjobs = status[errorcolumns[0]] != ""
 
diff --git a/src/fmu/ensemble/virtualensemble.py b/src/fmu/ensemble/virtualensemble.py
index 96bc867e..7db04dd5 100644
--- a/src/fmu/ensemble/virtualensemble.py
+++ b/src/fmu/ensemble/virtualensemble.py
@@ -255,8 +255,8 @@ def add_realization(self, realization, realidx=None, overwrite=False):
             if key not in self.data.keys():
                 self.data[key] = dframe
             else:
-                self.data[key] = self.data[key].append(
-                    dframe, ignore_index=True, sort=True
+                self.data[key] = pd.concat(
+                    [self.data[key], dframe], ignore_index=True, sort=True
                 )
         self.update_realindices()
 
@@ -1030,7 +1030,7 @@ def get_smry_meta(self, column_keys=None):
         return (
             self.get_df("__smry_metadata")
             .set_index("SMRYCOLUMN")
-            .loc[matches, :]
+            .loc[list(matches), :]
             .replace({np.nan: None})
             .to_dict(orient="index")
         )
diff --git a/src/fmu/ensemble/virtualrealization.py b/src/fmu/ensemble/virtualrealization.py
index ccdfc5a0..95f9a17f 100644
--- a/src/fmu/ensemble/virtualrealization.py
+++ b/src/fmu/ensemble/virtualrealization.py
@@ -480,7 +480,7 @@ def get_smry_meta(self, column_keys=None):
         return (
             self.get_df("__smry_metadata")
             .set_index("SMRYCOLUMN")
-            .loc[matches, :]
+            .loc[list(matches), :]
             .to_dict(orient="index")
         )
 
diff --git a/tests/test_ensemblecombination.py b/tests/test_ensemblecombination.py
index c8ee9353..8506851a 100644
--- a/tests/test_ensemblecombination.py
+++ b/tests/test_ensemblecombination.py
@@ -157,7 +157,7 @@ def test_ensemblecomb_observations():
     obs.load_smry(mean, "FOPT", time_index="yearly")
     mis = obs.mismatch(delta_ens)
     # Realization 4 is best representing the mean delta FOPT:
-    assert mis.groupby("REAL").sum()["L2"].sort_values().index.values[0] == 4
+    assert mis.groupby("REAL")["L2"].sum().sort_values().index.values[0] == 4
     # (this is the same as the representative realization for mean as
     # found in test_observations.py::test_virtual_observations)
 
diff --git a/tests/test_observations.py b/tests/test_observations.py
index 220f5d06..d7c3ec1e 100644
--- a/tests/test_observations.py
+++ b/tests/test_observations.py
@@ -679,7 +679,7 @@ def test_virtual_observations():
         mis = obs.mismatch(ens)
 
         closest_realization = (
-            mis.groupby("REAL").sum()["L2"].sort_values().index.values[0]
+            mis.groupby("REAL")["L2"].sum().sort_values().index.values[0]
         )
         representative_realizations[virtrealname] = closest_realization
 
@@ -707,7 +707,7 @@ def test_virtual_observations():
         mis = obs.mismatch(ens)
 
         closest_realization = (
-            mis.groupby("REAL").sum()["L2"].sort_values().index.values[0]
+            mis.groupby("REAL")["L2"].sum().sort_values().index.values[0]
         )
         vrepresentative_realizations[virtrealname] = closest_realization
 
diff --git a/tests/test_virtualensemble.py b/tests/test_virtualensemble.py
index eeaf7261..d299e8c2 100644
--- a/tests/test_virtualensemble.py
+++ b/tests/test_virtualensemble.py
@@ -258,8 +258,8 @@ def test_todisk(tmpdir):
         # Columns that only contains NaN will not have their
         # type preserved, this is too much to ask for, especially
         # with CSV files. So we drop columns with NaN
-        virtframe = vens.get_df(frame).dropna("columns")
-        diskframe = fromdisk.get_df(frame).dropna("columns")
+        virtframe = vens.get_df(frame).dropna(axis="columns")
+        diskframe = fromdisk.get_df(frame).dropna(axis="columns")
 
         # It would be nice to be able to use pd.Dataframe.equals,
         # but it is too strict, as columns with mixed type number/strings