diff --git a/SOAP/compute_halo_properties.py b/SOAP/compute_halo_properties.py
index cb8ff33a..7ee8f791 100644
--- a/SOAP/compute_halo_properties.py
+++ b/SOAP/compute_halo_properties.py
@@ -540,6 +540,12 @@ def compute_halo_properties():
                 "BoundSubhalo/EncloseRadius is not enabled. This means apertures with r > r_enclose will be calculated explicitly, rather than copying over values from smaller apertures"
             )
         category_filter.print_filters()
+        cosmology_errors = parameter_file.check_cosmology(cellgrid.cosmology)
+        if len(cosmology_errors):
+            print("The snapshot cosmology is incompatible with some properties:")
+            for error in cosmology_errors:
+                print(f"  {error}", flush=True)
+            comm_world.Abort(1)
 
         # Properties enabled in the parameter file must be computed, so abort
         # if the input files do not contain the datasets they require
diff --git a/SOAP/core/parameter_file.py b/SOAP/core/parameter_file.py
index da458cf3..c02740ec 100644
--- a/SOAP/core/parameter_file.py
+++ b/SOAP/core/parameter_file.py
@@ -588,6 +588,30 @@ def print_invalid_properties(self, halo_prop_list) -> None:
             for base_halo_type, prop in invalid_properties:
                 print(f"  {base_halo_type}  {prop}")
 
+    def check_cosmology(self, cosmology: Dict) -> List[str]:
+        """
+        Check that the cosmology is compatible with the properties that will
+        be calculated. This must be called after all the halo types have been
+        created, since that is when the property filters are set.
+
+        Parameters:
+         - cosmology: Dict
+           Cosmology attributes read from the snapshot.
+
+        Returns a list of error messages, which is empty if there are no problems.
+        """
+        errors = []
+
+        # The pseudo-evolution correction for the flow rates assumes flat LCDM
+        SO_filters = self.property_filters.get("SOProperties", {})
+        if any(f for name, f in SO_filters.items() if name.endswith("FlowRate")):
+            if abs(cosmology["Omega_k"]) > 1e-6:
+                errors.append("SO flow rates can only be computed if Omega_k=0")
+            if (cosmology["w_0"] != -1) or (cosmology["w_a"] != 0):
+                errors.append("SO flow rates can only be computed if w_0=-1, w_a=0")
+
+        return errors
+
     def has_enabled_properties(self, base_halo_type: str) -> bool:
         """
         Return True if the parameter file enables at least one property for the
diff --git a/SOAP/core/swift_cells.py b/SOAP/core/swift_cells.py
index bb3fbcbc..1f93d24c 100644
--- a/SOAP/core/swift_cells.py
+++ b/SOAP/core/swift_cells.py
@@ -142,6 +142,51 @@ def identify_datasets(filename, nr_files, ptypes, registry):
     return metadata
 
 
+def compute_virBN98(cosmology, a):
+    """
+    Compute the Bryan & Norman (1998) critical density multiple at the
+    given scale factor.
+
+    Parameters:
+     - cosmology: dict
+       Cosmology attributes read from the snapshot.
+     - a: float
+       Scale factor.
+
+    Returns the critical density multiple.
+    """
+    Omega_k = cosmology["Omega_k"]
+    Omega_Lambda = cosmology["Omega_lambda"]
+    Omega_m = cosmology["Omega_m"]
+    bnx = -(Omega_k / a**2 + Omega_Lambda) / (
+        Omega_k / a**2 + Omega_m / a**3 + Omega_Lambda
+    )
+    return 18.0 * np.pi**2 + 82.0 * bnx - 39.0 * bnx**2
+
+
+def compute_dlog_virBN98_dloga(cosmology, a, eps=1e-5):
+    """
+    Compute the logarithmic derivative of the Bryan & Norman (1998) critical
+    density multiple with respect to the scale factor. We difference the
+    expression used to set the multiple itself, so that the two cannot
+    become inconsistent.
+
+    Parameters:
+     - cosmology: dict
+       Cosmology attributes read from the snapshot.
+     - a: float
+       Scale factor.
+     - eps: float
+       Step size in log(a) used for the central difference.
+
+    Returns dlog(virBN98)/dlog(a).
+    """
+    return (
+        np.log(compute_virBN98(cosmology, a * np.exp(eps)))
+        - np.log(compute_virBN98(cosmology, a * np.exp(-eps)))
+    ) / (2 * eps)
+
+
 class SWIFTCellGrid:
     def get_unit(self, name):
         return unyt.Unit(name, registry=self.snap_unit_registry)
@@ -270,16 +315,14 @@ def __init__(
             )
 
             # Compute the BN98 critical density multiple
-            Omega_k = self.cosmology["Omega_k"]
-            Omega_Lambda = self.cosmology["Omega_lambda"]
-            Omega_m = self.cosmology["Omega_m"]
-            bnx = -(Omega_k / self.a**2 + Omega_Lambda) / (
-                Omega_k / self.a**2 + Omega_m / self.a**3 + Omega_Lambda
-            )
-            self.virBN98 = 18.0 * np.pi**2 + 82.0 * bnx - 39.0 * bnx**2
+            self.virBN98 = compute_virBN98(self.cosmology, self.a)
             if self.virBN98 < 50.0 or self.virBN98 > 1000.0:
                 raise RuntimeError("Invalid value for virBN98!")
 
+            # The BN98 density multiple is time dependent, so its logarithmic
+            # derivative is required to calculate the pseudo-evolution of R_BN98.
+            self.dlog_virBN98_dloga = compute_dlog_virBN98_dloga(self.cosmology, self.a)
+
             # Get the box size. Assume it's comoving with no h factors.
             comoving_length_unit = self.get_unit("snap_length") * self.a_unit
             self.boxsize = unyt.unyt_quantity(
diff --git a/SOAP/particle_selection/SO_properties.py b/SOAP/particle_selection/SO_properties.py
index cd9d5d73..aa0790ba 100644
--- a/SOAP/particle_selection/SO_properties.py
+++ b/SOAP/particle_selection/SO_properties.py
@@ -249,6 +249,7 @@ def __init__(
         observer_position: unyt.unyt_array,
         core_excision_fraction: float,
         virial_definition: bool,
+        compute_flow_rates: bool,
         search_radius: unyt.unyt_quantity,
         cosmology: dict,
         boxsize: unyt.unyt_quantity,
@@ -273,6 +274,9 @@ def __init__(
          - virial_definition: bool
            Whether to calculate the properties that are only valid for virial SO
            definitions
+         - compute_flow_rates: bool
+           Whether to calculate the flow rates. These are not valid for SO
+           definitions with a fixed physical radius.
          - search_radius: unyt.unyt_quantity
            Current search radius. Particles are guaranteed to be included up to
            this radius.
@@ -311,6 +315,7 @@ def __init__(
         self.observer_position = observer_position
         self.core_excision_fraction = core_excision_fraction
         self.virial_definition = virial_definition
+        self.compute_flow_rates = compute_flow_rates
         self.search_radius = search_radius
 
     def get_dataset(self, name: str) -> unyt.unyt_array:
@@ -2796,17 +2801,12 @@ def calculate_flow_rate(
             # Adding Hubble flow term
             if hubble:
                 v_r += radii[r_mask] * self.cosmology["H"]
-            # Account for expansion of R_SO
+            # Account for expansion of R_SO. The coefficient depends on the
+            # SO definition, and is set when this calculation is constructed.
             if pseudo_evolve:
-                G = unyt.Unit("newton_G", registry=masses.units.registry)
-                R_dot = (2 / 3) * (G * self.SO_mass * self.cosmology["H"] / 100) ** (
-                    1 / 3
+                v_r -= (
+                    R * self.cosmology["H"] * self.cosmology["pseudo_evolution_coeff"]
                 )
-                R_dot *= (
-                    2 * self.cosmology["Omega_g"] + (3 / 2) * self.cosmology["Omega_m"]
-                )
-                R_dot *= R_frac
-                v_r -= R_dot
 
             # Calculate different flow types
             # We want both the inflow and outflow rates to be positive values
@@ -2845,7 +2845,7 @@ def DarkMatterMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of dark matter through 3 spherical shells
         """
-        if (self.Ndm == 0) or (not self.virial_definition):
+        if (self.Ndm == 0) or (not self.compute_flow_rates):
             return None
 
         # Particles outside the SO radius are required to calculate the
@@ -2861,7 +2861,7 @@ def StellarMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of stars through 3 spherical shells
         """
-        if (self.Nstar == 0) or (not self.virial_definition):
+        if (self.Nstar == 0) or (not self.compute_flow_rates):
             return None
 
         # Particles outside the SO radius are required to calculate the
@@ -2877,7 +2877,7 @@ def HIMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of HI through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         # Particles outside the SO radius are required to calculate the
@@ -2903,7 +2903,7 @@ def H2MassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of H2 through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         # Particles outside the SO radius are required to calculate the
@@ -2931,7 +2931,7 @@ def MetalMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of metals through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         # Particles outside the SO radius are required to calculate the
@@ -2983,7 +2983,7 @@ def ColdGasMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of cold gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmax = 1.0e3 * unyt.K
@@ -2994,7 +2994,7 @@ def CoolGasMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of cool gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e3 * unyt.K
@@ -3008,7 +3008,7 @@ def WarmGasMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of warm gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e5 * unyt.K
@@ -3022,7 +3022,7 @@ def HotGasMassFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the mass flow rate of hot gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e7 * unyt.K
@@ -3033,7 +3033,7 @@ def ColdGasEnergyFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the energy flow rate of cold gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmax = 1.0e3 * unyt.K
@@ -3046,7 +3046,7 @@ def CoolGasEnergyFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the energy flow rate of cool gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e3 * unyt.K
@@ -3060,7 +3060,7 @@ def WarmGasEnergyFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the energy flow rate of warm gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e5 * unyt.K
@@ -3074,7 +3074,7 @@ def HotGasEnergyFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the energy flow rate of hot gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e7 * unyt.K
@@ -3087,7 +3087,7 @@ def ColdGasMomentumFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the momentum flow rate of cold gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmax = 1.0e3 * unyt.K
@@ -3100,7 +3100,7 @@ def CoolGasMomentumFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the momentum flow rate of cool gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e3 * unyt.K
@@ -3114,7 +3114,7 @@ def WarmGasMomentumFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the momentum flow rate of warm gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e5 * unyt.K
@@ -3128,7 +3128,7 @@ def HotGasMomentumFlowRate(self) -> unyt.unyt_array:
         """
         Calculate the momentum flow rate of hot gas through 3 spherical shells
         """
-        if (self.Ngas == 0) or (not self.virial_definition):
+        if (self.Ngas == 0) or (not self.compute_flow_rates):
             return None
 
         Tmin = 1.0e7 * unyt.K
@@ -3375,8 +3375,6 @@ def __init__(
         self.cosmology["H"] = cellgrid.cosmology[
             "H [internal units]"
         ] / cellgrid.get_unit("code_time")
-        self.cosmology["Omega_g"] = cellgrid.cosmology["Omega_g"]
-        self.cosmology["Omega_m"] = cellgrid.cosmology["Omega_m"]
 
         # This specifies how large a sphere is read in:
         # we use default values that are sufficiently small/large to avoid reading in too many particles
@@ -3397,6 +3395,31 @@ def __init__(
             self.virial_definition = True
         elif type == "physical":
             self.physical_radius_mpc = 0.001 * SOval
+        # Flow rates are not computed for a fixed physical radius, since it
+        # does not pseudo-evolve
+        self.compute_flow_rates = type != "physical"
+
+        # Coefficient used to correct the flow rates for the pseudo-evolution
+        # of the SO radius: Rdot = coeff * R * H, where
+        # coeff = -(1/3) dln(rho_ref)/dln(a) at fixed SO mass, and rho_ref is
+        # the reference density used to define the SO radius.
+        # Derivation is in documentation/pseudo_evolution.pdf
+        # The Omega values in the snapshot are z=0 values, so we scale them.
+        H0_over_H_sq = (
+            cellgrid.cosmology["H0 [internal units]"]
+            / cellgrid.cosmology["H [internal units]"]
+        ) ** 2
+        Omega_m = float(cellgrid.mean_density / cellgrid.critical_density)
+        Omega_r = cellgrid.cosmology["Omega_r"] * H0_over_H_sq / cellgrid.a**4
+        one_plus_q = 2 * Omega_r + 1.5 * Omega_m
+        if type == "mean":
+            self.cosmology["pseudo_evolution_coeff"] = 1.0
+        elif type == "crit":
+            self.cosmology["pseudo_evolution_coeff"] = (2 / 3) * one_plus_q
+        elif type == "BN98":
+            self.cosmology["pseudo_evolution_coeff"] = (
+                2 * one_plus_q - cellgrid.dlog_virBN98_dloga
+            ) / 3
 
         # Give this calculation a name so we can select it on the command line
         if type in ["mean", "crit"]:
@@ -3573,6 +3596,7 @@ def calculate(
                 self.observer_position,
                 self.core_excision_fraction,
                 self.virial_definition,
+                self.compute_flow_rates,
                 search_radius,
                 self.cosmology,
                 self.boxsize,
diff --git a/SOAP/property_table.py b/SOAP/property_table.py
index ca4e52b6..0a8c004b 100644
--- a/SOAP/property_table.py
+++ b/SOAP/property_table.py
@@ -3417,7 +3417,7 @@ class PropertyTable:
             shape=6,
             dtype=np.float32,
             unit="snap_mass / snap_time",
-            description="Mass flow rate of dark matter particles through spherical shells. Contains 6 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R.",
+            description="Mass flow rate of dark matter particles through spherical shells. Contains 6 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=True,
             particle_properties=[
@@ -3433,7 +3433,7 @@ class PropertyTable:
             shape=9,
             dtype=np.float32,
             unit="snap_mass / snap_time",
-            description="Mass flow rate of cold gas particles ($\\log T < 3$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R.",
+            description="Mass flow rate of cold gas particles ($\\log T < 3$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3450,7 +3450,7 @@ class PropertyTable:
             shape=9,
             dtype=np.float32,
             unit="snap_mass / snap_time",
-            description="Mass flow rate of cool gas particles ($3 < \\log T < 5$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R.",
+            description="Mass flow rate of cool gas particles ($3 < \\log T < 5$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3467,7 +3467,7 @@ class PropertyTable:
             shape=9,
             dtype=np.float32,
             unit="snap_mass / snap_time",
-            description="Mass flow rate of warm gas particles ($5 < \\log T < 7$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R.",
+            description="Mass flow rate of warm gas particles ($5 < \\log T < 7$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3484,7 +3484,7 @@ class PropertyTable:
             shape=9,
             dtype=np.float32,
             unit="snap_mass / snap_time",
-            description="Mass flow rate of hot gas particles ($7 < \\log T$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R.",
+            description="Mass flow rate of hot gas particles ($7 < \\log T$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3501,7 +3501,7 @@ class PropertyTable:
             shape=6,
             dtype=np.float32,
             unit="snap_mass / snap_time",
-            description="Mass flow rate of gas particles through spherical shells weighted by HI fraction. Contains 6 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R.",
+            description="Mass flow rate of gas particles through spherical shells weighted by HI fraction. Contains 6 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3519,7 +3519,7 @@ class PropertyTable:
             shape=6,
             dtype=np.float32,
             unit="snap_mass / snap_time",
-            description="Mass flow rate of gas particles through spherical shells weighted by H2 fraction. Does not include Helium. Contains 6 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R.",
+            description="Mass flow rate of gas particles through spherical shells weighted by H2 fraction. Does not include Helium. Contains 6 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3537,7 +3537,7 @@ class PropertyTable:
             shape=6,
             dtype=np.float32,
             unit="snap_mass / snap_time",
-            description="Mass flow rate of gas particles through spherical shells weighted by metal fraction. Contains 6 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R.",
+            description="Mass flow rate of gas particles through spherical shells weighted by metal fraction. Contains 6 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3554,7 +3554,7 @@ class PropertyTable:
             shape=6,
             dtype=np.float32,
             unit="snap_mass / snap_time",
-            description="Mass flow rate of star particles through spherical shells. Contains 6 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R.",
+            description="Mass flow rate of star particles through spherical shells. Contains 6 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3570,7 +3570,7 @@ class PropertyTable:
             shape=9,
             dtype=np.float32,
             unit="snap_mass*snap_length**2/snap_time**3",
-            description="Energy flow rate of cold gas particles ($\\log T < 3$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R.",
+            description="Energy flow rate of cold gas particles ($\\log T < 3$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3588,7 +3588,7 @@ class PropertyTable:
             shape=9,
             dtype=np.float32,
             unit="snap_mass*snap_length**2/snap_time**3",
-            description="Energy flow rate of cool gas particles ($3 < \\log T < 5$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R.",
+            description="Energy flow rate of cool gas particles ($3 < \\log T < 5$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3606,7 +3606,7 @@ class PropertyTable:
             shape=9,
             dtype=np.float32,
             unit="snap_mass*snap_length**2/snap_time**3",
-            description="Energy flow rate of warm gas particles ($5 < \\log T < 7$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R.",
+            description="Energy flow rate of warm gas particles ($5 < \\log T < 7$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3624,7 +3624,7 @@ class PropertyTable:
             shape=9,
             dtype=np.float32,
             unit="snap_mass*snap_length**2/snap_time**3",
-            description="Energy flow rate of hot gas particles ($7 < \\log T$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R.",
+            description="Energy flow rate of hot gas particles ($7 < \\log T$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3642,7 +3642,7 @@ class PropertyTable:
             shape=9,
             dtype=np.float32,
             unit="snap_mass*snap_length/snap_time**2",
-            description="Momentum flow rate of cold gas particles ($\\log T < 3$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R.",
+            description="Momentum flow rate of cold gas particles ($\\log T < 3$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3660,7 +3660,7 @@ class PropertyTable:
             shape=9,
             dtype=np.float32,
             unit="snap_mass*snap_length/snap_time**2",
-            description="Momentum flow rate of cool gas particles ($3 < \\log T < 5$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R.",
+            description="Momentum flow rate of cool gas particles ($3 < \\log T < 5$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3678,7 +3678,7 @@ class PropertyTable:
             shape=9,
             dtype=np.float32,
             unit="snap_mass*snap_length/snap_time**2",
-            description="Momentum flow rate of warm gas particles ($5 < \\log T < 7$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R.",
+            description="Momentum flow rate of warm gas particles ($5 < \\log T < 7$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
@@ -3696,7 +3696,7 @@ class PropertyTable:
             shape=9,
             dtype=np.float32,
             unit="snap_mass*snap_length/snap_time**2",
-            description="Momentum flow rate of hot gas particles ($7 < \\log T$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R.",
+            description="Momentum flow rate of hot gas particles ($7 < \\log T$) through spherical shells. Contains 9 entries: inflow rate at 0.1R, 0.3R, R, outflow rate at 0.1R, 0.3R, R, fast outflow rate at 0.1R, 0.3R, R. Pseudo-evolution correction applied to every SO radius.",
             lossy_compression_filter="FMantissa9",
             dmo_property=False,
             particle_properties=[
diff --git a/documentation/footnote_flow_rates.tex b/documentation/footnote_flow_rates.tex
index 3c852f7f..67218e42 100644
--- a/documentation/footnote_flow_rates.tex
+++ b/documentation/footnote_flow_rates.tex
@@ -7,16 +7,35 @@
 \begin{equation}
     v_{r,i} = (\underline{v_i} - \underline{v_{COM}}) \cdot \underline{\hat{r}} - \dot{R}
 \end{equation}
-where $\underline{v_i}$ is the velocity of the particle, $\underline{v_{COM}}$ is the centre of mass velocity of all particles within $R$ (therefore we use a different value for each spherical shell). The final term accounts for the "pseudo-evolution" of the halo radius and is given by
+where $\underline{v_i}$ is the velocity of the particle, $\underline{v_{COM}}$ is the
+centre of mass velocity of all particles within $R$ (therefore we use a different value
+for each spherical shell). The final term accounts for the "pseudo-evolution" of the halo
+radius and is given by
 
 \begin{equation}
-    \dot{R} = f \frac{2}{3} \left(\frac{GHM_{SO}}{100}\right)^\frac{1}{3} \left( 2 \Omega_\gamma + \frac{3}{2} \Omega_m \right)
+    \dot{R} = -\frac{1}{3} f R_{SO} H \frac{\mathrm{d} \ln \rho_{\rm ref}}{\mathrm{d} \ln a},
 \end{equation}
+where $\rho_{\rm ref}$ is the reference density used to define the SO radius. Since
+this depends on the SO definition, we can write $\dot{R} = c f R_{SO} H$ with
 
-This is required since accretion rates are often measured by subtracting the halo mass between consecutive snapshots and dividing by the time interval.
-To be consistent with this method we must consider that the virial radius is defined w.r.t background density (which decreases in time). Hence, the virial radius actually moves
+\begin{equation}
+    c = \begin{cases}
+        1 & \Delta_{m} \; ({\rm e.g.} \; 200_{m}) \\[4pt]
+        \frac{2}{3} \left( 2 \Omega_r + \frac{3}{2} \Omega_m \right) & \Delta_{c} \; ({\rm e.g.} \; 200_{c}, 500_{c}) \\[4pt]
+        \frac{1}{3} \left[ 2 \left( 2 \Omega_r + \frac{3}{2} \Omega_m \right) - \frac{\mathrm{d} \ln \Delta_{BN}}{\mathrm{d} \ln a} \right] & BN98
+    \end{cases}
+\end{equation}
+
+where $\Delta_{BN}$ is the Bryan \& Norman (1998) critical density multiple. The density
+parameters are evaluated at the redshift of the snapshot, not at $z=0$.
+
+This is required since accretion rates are often measured by subtracting
+the halo mass between consecutive snapshots and dividing by the time interval.
+To be consistent with this method we must consider that the virial radius is
+defined w.r.t background density (which decreases in time). Hence, the virial radius actually moves
 outward with a velocity that we can compute analytically. This means that static
-particles at the virial radius actually become inflowing. The expression is derived by taking the partial differential of the analytic expression for $R_{200}$ w.r.t. time.
+particles at the virial radius actually become inflowing. The expression is derived by taking the partial
+differential of the analytic expression for $R_{SO}$ w.r.t. time at fixed $M_{SO}$.
 Note that $v_{r,i}$ does not include the Hubble flow relative to the halo centre.
 
 To calculate the mass inflow (outflow) rate
@@ -40,5 +59,10 @@
     \frac{1}{dR} \sum_{i} m_i \left(v_{r,i}^2 + \frac{c_s^2}{\gamma}\right),
 \end{equation}
 
-where $c_s$ is the sound speed and $\gamma$ = 5/3 (the second term accounts for pressure). For the gas phases we also calculate "fast outflow" rates. These are calculated by using the equations above, but only for particles that satisfy $v_{r,i} > V_{max} / 4$, where $V_{max}$ is the maximum circular velocity of the halo. The flow rates are always positive, so to compute the net rate you must subtract the inflow rate from the outflow rate. Flow rates are only calculated for the
-following SO definitions: $200_{c}$, $200_{m}$, $BN98$. To calculate the total gas flow rate the individual phases should be summed together.
+where $c_s$ is the sound speed and $\gamma$ = 5/3 (the second term accounts for pressure).
+For the gas phases we also calculate "fast outflow" rates.
+These are calculated by using the equations above, but only for particles that satisfy $v_{r,i} > V_{max} / 4$,
+where $V_{max}$ is the maximum circular velocity of the halo. The flow rates are always positive,
+so to compute the net rate you must subtract the inflow rate from the outflow rate.
+Flow rates are calculated for all SO definitions except those with a radius that is a multiple of another SO radius.
+To calculate the total gas flow rate the individual phases should be summed together.
diff --git a/documentation/pseudo_evolution.pdf b/documentation/pseudo_evolution.pdf
new file mode 100644
index 00000000..74f7575b
Binary files /dev/null and b/documentation/pseudo_evolution.pdf differ
diff --git a/misc/convert_eagle.py b/misc/convert_eagle.py
index c5df9ab8..d6b94d7d 100644
--- a/misc/convert_eagle.py
+++ b/misc/convert_eagle.py
@@ -232,8 +232,8 @@
             "Omega_m": header.attrs["Omega0"],
             "Omega_k": 0,
             "Omega_nu_0": 0,
-            "Omega_r": astropy.cosmology.Planck13.Ogamma(z),
-            "Omega_g": astropy.cosmology.Planck13.Ogamma(z),
+            "Omega_r": astropy.cosmology.Planck13.Ogamma0,
+            "Omega_g": astropy.cosmology.Planck13.Ogamma0,
             "Omega_lambda": header.attrs["OmegaLambda"],
             "Redshift": z,
             "H0 [internal units]": h * 100,
diff --git a/tests/dummy_halo_generator.py b/tests/dummy_halo_generator.py
index 4fcf4397..4c912927 100644
--- a/tests/dummy_halo_generator.py
+++ b/tests/dummy_halo_generator.py
@@ -21,6 +21,7 @@
 
 from SOAP.core.swift_units import unit_registry_from_snapshot
 from SOAP.core.snapshot_datasets import SnapshotDatasets
+from SOAP.core.swift_cells import compute_virBN98, compute_dlog_virBN98_dloga
 from SOAP.property_table import PropertyTable
 from SOAP.particle_filter.recently_heated_gas_filter import RecentlyHeatedGasFilter
 from SOAP.particle_filter.cold_dense_gas_filter import ColdDenseGasFilter
@@ -475,16 +476,14 @@ def __init__(self, reg: unyt.UnitRegistry, snap: h5py.File):
         )
         self.mean_density = self.critical_density * self.cosmology["Omega_m"]
         # Compute the BN98 critical density multiple
-        Omega_k = self.cosmology["Omega_k"]
-        Omega_Lambda = self.cosmology["Omega_lambda"]
-        Omega_m = self.cosmology["Omega_m"]
-        bnx = -(Omega_k / self.a**2 + Omega_Lambda) / (
-            Omega_k / self.a**2 + Omega_m / self.a**3 + Omega_Lambda
-        )
-        self.virBN98 = 18.0 * np.pi**2 + 82.0 * bnx - 39.0 * bnx**2
+        self.virBN98 = compute_virBN98(self.cosmology, self.a)
         if self.virBN98 < 50.0 or self.virBN98 > 1000.0:
             raise RuntimeError("Invalid value for virBN98!")
 
+        # The BN98 density multiple is time dependent, so its logarithmic
+        # derivative is required to calculate the pseudo-evolution of R_BN98
+        self.dlog_virBN98_dloga = compute_dlog_virBN98_dloga(self.cosmology, self.a)
+
         # Get the box size. Assume it's comoving with no h factors.
         comoving_length_unit = self.get_unit("snap_length") * self.a_unit
         self.boxsize = unyt.unyt_quantity(100.0, units=comoving_length_unit)