diff --git a/gusto/core/configuration.py b/gusto/core/configuration.py
index 4f82a35f7..40852033a 100644
--- a/gusto/core/configuration.py
+++ b/gusto/core/configuration.py
@@ -269,6 +269,19 @@ class BoundaryLayerParameters(Configuration):
     mu = 100.                   # Interior penalty coefficient for vertical diffusion
 
 
+class HeldSuarezParameters(Configuration):
+    """
+    Parameters used in the default configuration for the Held Suarez test case.
+    """
+    T0stra = 200               # Stratosphere temp
+    T0surf = 315               # Surface temperature at equator
+    T0horiz = 60               # Equator to pole temperature difference
+    T0vert = 10                # Stability parameter
+    sigmab = 0.7               # Height of the boundary layer
+    tau_d = 40 * 24 * 60 * 60  # 40 day time scale
+    tau_u = 4 * 24 * 60 * 60   # 4 day timescale
+
+
 class SubcyclingOptions(Configuration):
     """
     Describes the process of subcycling a time discretisation, by dividing the
diff --git a/gusto/diagnostics/diagnostics.py b/gusto/diagnostics/diagnostics.py
index 90192f022..00272edac 100644
--- a/gusto/diagnostics/diagnostics.py
+++ b/gusto/diagnostics/diagnostics.py
@@ -1,14 +1,15 @@
 """Common diagnostic fields."""
 
 
-from firedrake import (assemble, dot, dx, Function, sqrt, TestFunction,
+from firedrake import (dot, dx, Function, sqrt, TestFunction,
                        TrialFunction, Constant, grad, inner, FacetNormal,
                        LinearVariationalProblem, LinearVariationalSolver,
                        ds_b, ds_v, ds_t, dS_h, dS_v, ds, dS, div, avg, pi,
                        TensorFunctionSpace, SpatialCoordinate, as_vector,
-                       Projector, Interpolator, FunctionSpace, FiniteElement,
+                       Projector, assemble, FunctionSpace, FiniteElement,
                        TensorProductElement)
 from firedrake.assign import Assigner
+from firedrake.__future__ import interpolate
 from ufl.domain import extract_unique_domain
 
 from abc import ABCMeta, abstractmethod, abstractproperty
@@ -193,7 +194,7 @@ def setup(self, domain, state_fields, space=None):
 
             # Solve method must be declared in diagnostic's own setup routine
             if self.method == 'interpolate':
-                self.evaluator = Interpolator(self.expr, self.field)
+                self.evaluator = interpolate(self.expr, self.space)
             elif self.method == 'project':
                 self.evaluator = Projector(self.expr, self.field)
             elif self.method == 'assign':
@@ -207,7 +208,7 @@ def compute(self):
         logger.debug(f'Computing diagnostic {self.name} with {self.method} method')
 
         if self.method == 'interpolate':
-            self.evaluator.interpolate()
+            self.field.assign(assemble(self.evaluator))
         elif self.method == 'assign':
             self.evaluator.assign()
         elif self.method == 'project':
@@ -294,7 +295,7 @@ def setup(self, domain, state_fields, space=None):
 
             # Solve method must be declared in diagnostic's own setup routine
             if self.method == 'interpolate':
-                self.evaluator = Interpolator(self.expr, self.field)
+                self.evaluator = interpolate(self.expr, self.space)
             elif self.method == 'project':
                 self.evaluator = Projector(self.expr, self.field)
             elif self.method == 'assign':
diff --git a/gusto/diagnostics/shallow_water_diagnostics.py b/gusto/diagnostics/shallow_water_diagnostics.py
index eb2dc4a32..024e58ecf 100644
--- a/gusto/diagnostics/shallow_water_diagnostics.py
+++ b/gusto/diagnostics/shallow_water_diagnostics.py
@@ -2,9 +2,10 @@
 
 
 from firedrake import (
-    dx, TestFunction, TrialFunction, grad, inner, curl, Function, Interpolator,
+    dx, TestFunction, TrialFunction, grad, inner, curl, Function, assemble,
     LinearVariationalProblem, LinearVariationalSolver, conditional
 )
+from firedrake.__future__ import interpolate
 from gusto.diagnostics.diagnostics import DiagnosticField, Energy
 
 __all__ = ["ShallowWaterKineticEnergy", "ShallowWaterPotentialEnergy",
@@ -321,14 +322,14 @@ def setup(self, domain, state_fields):
 
         qsat_expr = self.equation.compute_saturation(state_fields.X(
             self.equation.field_name))
-        self.qsat_interpolator = Interpolator(qsat_expr, self.qsat_func)
+        self.qsat_interpolate = interpolate(qsat_expr, space)
         self.expr = conditional(q_t < self.qsat_func, q_t, self.qsat_func)
 
         super().setup(domain, state_fields, space=space)
 
     def compute(self):
         """Performs the computation of the diagnostic field."""
-        self.qsat_interpolator.interpolate()
+        self.qsat_func.assign(assemble(self.qsat_interpolate))
         super().compute()
 
 
@@ -371,7 +372,7 @@ def setup(self, domain, state_fields):
 
         qsat_expr = self.equation.compute_saturation(state_fields.X(
             self.equation.field_name))
-        self.qsat_interpolator = Interpolator(qsat_expr, self.qsat_func)
+        self.qsat_interpolate = interpolate(qsat_expr, space)
         vapour = conditional(q_t < self.qsat_func, q_t, self.qsat_func)
         self.expr = q_t - vapour
 
@@ -379,5 +380,5 @@ def setup(self, domain, state_fields):
 
     def compute(self):
         """Performs the computation of the diagnostic field."""
-        self.qsat_interpolator.interpolate()
+        self.qsat_func.assign(assemble(self.qsat_interpolate))
         super().compute()
diff --git a/gusto/physics/__init__.py b/gusto/physics/__init__.py
index 91f09c2ec..b1407a299 100644
--- a/gusto/physics/__init__.py
+++ b/gusto/physics/__init__.py
@@ -2,4 +2,5 @@
 from gusto.physics.chemistry import *                       # noqa
 from gusto.physics.boundary_and_turbulence import *         # noqa
 from gusto.physics.microphysics import *                    # noqa
-from gusto.physics.shallow_water_microphysics import *      # noqa
\ No newline at end of file
+from gusto.physics.shallow_water_microphysics import *      # noqa
+from gusto.physics.held_suarez_forcing import *             # noqa
diff --git a/gusto/physics/held_suarez_forcing.py b/gusto/physics/held_suarez_forcing.py
new file mode 100644
index 000000000..d7573d009
--- /dev/null
+++ b/gusto/physics/held_suarez_forcing.py
@@ -0,0 +1,188 @@
+import numpy as np
+from firedrake import (Interpolator, Function, dx, pi, SpatialCoordinate,
+                       split, conditional, ge, sin, dot, ln, cos, inner, Projector)
+from firedrake.fml import subject
+from gusto.core.coord_transforms import lonlatr_from_xyz
+from gusto.recovery import Recoverer, BoundaryMethod
+from gusto.physics.physics_parametrisation import PhysicsParametrisation
+from gusto.core.labels import prognostic
+from gusto.equations import thermodynamics
+from gusto.core.configuration import HeldSuarezParameters
+from gusto.core import logger
+
+
+class Relaxation(PhysicsParametrisation):
+    """
+    Relaxation term for Held Suarez
+    """
+
+    def __init__(self, equation, variable_name, parameters, hs_parameters=None):
+        """
+        Args:
+            equation (:class:`PrognosticEquationSet`): the model's equation.
+            variable_name (str): the name of the variable
+            hs_parameters (:class'Configuration'): contains the parameters for the Held-suariez test case
+
+        """
+        label_name = f'relaxation_{variable_name}'
+        if hs_parameters is None:
+            hs_parameters = HeldSuarezParameters()
+            logger.warning('Using default Held-Suarez parameters')
+        super().__init__(equation, label_name, hs_parameters)
+
+        if equation.domain.on_sphere:
+            x, y, z = SpatialCoordinate(equation.domain.mesh)
+            _, lat, _ = lonlatr_from_xyz(x, y, z)
+        else:
+            # TODO: this could be determined some other way
+            # Take a mid-latitude
+            lat = pi / 4
+
+        self.X = Function(equation.X.function_space())
+        X = self.X
+        self.domain = equation.domain
+        theta_idx = equation.field_names.index('theta')
+        self.theta = X.subfunctions[theta_idx]
+        Vt = equation.domain.spaces('theta')
+        rho_idx = equation.field_names.index('rho')
+        rho = split(X)[rho_idx]
+
+        boundary_method = BoundaryMethod.extruded if equation.domain.vertical_degree == 0 else None
+        self.rho_averaged = Function(Vt)
+        self.rho_recoverer = Recoverer(rho, self.rho_averaged, boundary_method=boundary_method)
+        self.exner = Function(Vt)
+        self.exner_interpolator = Interpolator(
+            thermodynamics.exner_pressure(equation.parameters,
+                                          self.rho_averaged, self.theta), self.exner)
+        self.sigma = Function(Vt)
+        kappa = equation.parameters.kappa
+
+        T0surf = hs_parameters.T0surf
+        T0horiz = hs_parameters.T0horiz
+        T0vert = hs_parameters.T0vert
+        T0stra = hs_parameters.T0stra
+
+        sigma_b = hs_parameters.sigmab
+        tau_d = hs_parameters.tau_d
+        tau_u = hs_parameters.tau_u
+
+        theta_condition = (T0surf - T0horiz * sin(lat)**2 - (T0vert * ln(self.exner) * cos(lat)**2)/kappa)
+        Theta_eq = conditional(T0stra/self.exner >= theta_condition, T0stra/self.exner, theta_condition)
+
+        # timescale of temperature forcing
+        tau_cond = (self.sigma**(1/kappa) - sigma_b) / (1 - sigma_b)
+        newton_freq = 1 / tau_d + (1/tau_u - 1/tau_d) * conditional(0 >= tau_cond, 0, tau_cond) * cos(lat)**4
+        forcing_expr = newton_freq * (self.theta - Theta_eq)
+
+        # Create source for forcing
+        self.source_relaxation = Function(Vt)
+        self.source_interpolator = Interpolator(forcing_expr, self.source_relaxation)
+
+        # Add relaxation term to residual
+        test = equation.tests[theta_idx]
+        dx_reduced = dx(degree=equation.domain.max_quad_degree)
+        forcing_form = test * self.source_relaxation * dx_reduced
+        equation.residual += self.label(subject(prognostic(forcing_form, 'theta'), X), self.evaluate)
+
+    def evaluate(self, x_in, dt):
+        """
+        Evalutes the source term generated by the physics.
+
+        Args:
+            x_in: (:class:`Function`): the (mixed) field to be evolved.
+            dt: (:class:`Constant`): the timestep, which can be the time
+                interval for the scheme.
+        """
+        self.X.assign(x_in)
+        self.rho_recoverer.project()
+        self.exner_interpolator.interpolate()
+
+        # Determine sigma:= exner / exner_surf
+        exner_columnwise, index_data = self.domain.coords.get_column_data(self.exner, self.domain)
+        sigma_columnwise = np.zeros_like(exner_columnwise)
+        for col in range(len(exner_columnwise[:, 0])):
+            sigma_columnwise[col, :] = exner_columnwise[col, :] / exner_columnwise[col, 0]
+        self.domain.coords.set_field_from_column_data(self.sigma, sigma_columnwise, index_data)
+
+        self.source_interpolator.interpolate()
+
+
+class RayleighFriction(PhysicsParametrisation):
+    """
+    Forcing term on the velocity of the form
+    F_u = -u / a,
+    where a is some friction factor
+    """
+    def __init__(self, equation, hs_parameters=None):
+        """
+         Args:
+            equation (:class:`PrognosticEquationSet`): the model's equation.
+            hs_parameters (:class'Configuration'): contains the parameters for the Held-suariez test case
+        """
+        label_name = 'rayleigh_friction'
+        if hs_parameters is None:
+            hs_parameters = HeldSuarezParameters()
+            logger.warning('Using default Held-Suarez parameters')
+        super().__init__(equation, label_name, hs_parameters)
+
+        self.domain = equation.domain
+        self.X = Function(equation.X.function_space())
+        X = self.X
+        k = equation.domain.k
+        u_idx = equation.field_names.index('u')
+        u = split(X)[u_idx]
+        theta_idx = equation.field_names.index('theta')
+        self.theta = X.subfunctions[theta_idx]
+        rho_idx = equation.field_names.index('rho')
+        rho = split(X)[rho_idx]
+        Vt = equation.domain.spaces('theta')
+        Vu = equation.domain.spaces('HDiv')
+        u_hori = u - k*dot(u, k)
+
+        boundary_method = BoundaryMethod.extruded if self.domain == 0 else None
+        self.rho_averaged = Function(Vt)
+        self.exner = Function(Vt)
+        self.rho_recoverer = Recoverer(rho, self.rho_averaged, boundary_method=boundary_method)
+        self.exner_interpolator = Interpolator(
+            thermodynamics.exner_pressure(equation.parameters,
+                                          self.rho_averaged, self.theta), self.exner)
+
+        self.sigma = Function(Vt)
+        sigmab = hs_parameters.sigmab
+        kappa = equation.parameters.kappa
+        tau_fric = 24 * 60 * 60
+
+        tau_cond = (self.sigma**(1/kappa) - sigmab) / (1 - sigmab)
+        wind_timescale = conditional(ge(0, tau_cond), 0, tau_cond) / tau_fric
+        forcing_expr = u_hori * wind_timescale
+
+        self.source_friction = Function(Vu)
+        self.source_projector = Projector(forcing_expr, self.source_friction)
+
+        tests = equation.tests
+        test = tests[u_idx]
+        dx_reduced = dx(degree=equation.domain.max_quad_degree)
+        source_form = inner(test, self.source_friction) * dx_reduced
+        equation.residual += self.label(subject(prognostic(source_form, 'u'), X), self.evaluate)
+
+    def evaluate(self, x_in, dt):
+        """
+        Evaluates the source term generated by the physics. This does nothing if
+        the implicit formulation is not used.
+
+        Args:
+            x_in: (:class: 'Function'): the (mixed) field to be evolved.
+            dt: (:class: 'Constant'): the timestep, which can be the time
+                interval for the scheme.
+        """
+        self.X.assign(x_in)
+        self.rho_recoverer.project()
+        self.exner_interpolator.interpolate()
+        # Determine sigma:= exner / exner_surf
+        exner_columnwise, index_data = self.domain.coords.get_column_data(self.exner, self.domain)
+        sigma_columnwise = np.zeros_like(exner_columnwise)
+        for col in range(len(exner_columnwise[:, 0])):
+            sigma_columnwise[col, :] = exner_columnwise[col, :] / exner_columnwise[col, 0]
+        self.domain.coords.set_field_from_column_data(self.sigma, sigma_columnwise, index_data)
+
+        self.source_projector.project()
diff --git a/gusto/physics/microphysics.py b/gusto/physics/microphysics.py
index e8313ab11..a5cbe3da5 100644
--- a/gusto/physics/microphysics.py
+++ b/gusto/physics/microphysics.py
@@ -4,9 +4,10 @@
 """
 
 from firedrake import (
-    Interpolator, conditional, Function, dx, min_value, max_value, Constant, pi,
-    Projector
+    conditional, Function, dx, min_value, max_value, Constant, pi,
+    Projector, assemble
 )
+from firedrake.__future__ import interpolate
 from firedrake.fml import identity, Term, subject
 from gusto.equations import Phases, TracerVariableType
 from gusto.recovery import Recoverer, BoundaryMethod
@@ -170,8 +171,7 @@ def __init__(self, equation, vapour_name='water_vapour',
         # Add terms to equations and make interpolators
         # -------------------------------------------------------------------- #
         self.source = [Function(V) for factor in factors]
-        self.source_interpolators = [Interpolator(sat_adj_expr*factor, source)
-                                     for factor, source in zip(factors, self.source)]
+        self.source_interpolate = [interpolate(sat_adj_expr*factor, V) for factor in factors]
 
         tests = [equation.tests[idx] for idx in V_idxs]
 
@@ -195,8 +195,8 @@ def evaluate(self, x_in, dt):
         if isinstance(self.equation, CompressibleEulerEquations):
             self.rho_recoverer.project()
         # Evaluate the source
-        for interpolator in self.source_interpolators:
-            interpolator.interpolate()
+        for interpolator, src in zip(self.source_interpolate, self.source):
+            src.assign(assemble(interpolator))
 
 
 class AdvectedMoments(Enum):
@@ -440,7 +440,7 @@ def __init__(self, equation, cloud_name='cloud_water', rain_name='rain',
                                                         min_value(accu_rate, self.cloud_water / self.dt),
                                                         min_value(accr_rate + accu_rate, self.cloud_water / self.dt))))
 
-        self.source_interpolator = Interpolator(rain_expr, self.source)
+        self.source_interpolate = interpolate(rain_expr, Vt)
 
         # Add term to equation's residual
         test_cl = equation.tests[self.cloud_idx]
@@ -464,7 +464,7 @@ def evaluate(self, x_in, dt):
         self.rain.assign(x_in.subfunctions[self.rain_idx])
         self.cloud_water.assign(x_in.subfunctions[self.cloud_idx])
         # Evaluate the source
-        self.source.assign(self.source_interpolator.interpolate())
+        self.source.assign(assemble(self.source_interpolate))
 
 
 class EvaporationOfRain(PhysicsParametrisation):
@@ -609,8 +609,7 @@ def __init__(self, equation, rain_name='rain', vapour_name='water_vapour',
         # Add terms to equations and make interpolators
         # -------------------------------------------------------------------- #
         self.source = [Function(V) for factor in factors]
-        self.source_interpolators = [Interpolator(evap_rate*factor, source)
-                                     for factor, source in zip(factors, self.source)]
+        self.source_interpolate = [interpolate(evap_rate*factor, V) for factor in factors]
 
         tests = [equation.tests[idx] for idx in V_idxs]
 
@@ -634,5 +633,5 @@ def evaluate(self, x_in, dt):
         if isinstance(self.equation, CompressibleEulerEquations):
             self.rho_recoverer.project()
         # Evaluate the source
-        for interpolator in self.source_interpolators:
-            interpolator.interpolate()
+        for interpolator, src in zip(self.source_interpolate, self.source):
+            src.assign(assemble(interpolator))
diff --git a/gusto/physics/physics_parametrisation.py b/gusto/physics/physics_parametrisation.py
index 2d43dae48..c304c23f6 100644
--- a/gusto/physics/physics_parametrisation.py
+++ b/gusto/physics/physics_parametrisation.py
@@ -8,7 +8,8 @@
 """
 
 from abc import ABCMeta, abstractmethod
-from firedrake import Interpolator, Function, dx, Projector
+from firedrake import Function, dx, Projector, assemble
+from firedrake.__future__ import interpolate
 from firedrake.fml import subject
 from gusto.core.labels import PhysicsLabel
 from gusto.core.logging import logger
@@ -117,14 +118,14 @@ def __init__(self, equation, variable_name, rate_expression,
 
         # Handle method of evaluating source/sink
         if self.method == 'interpolate':
-            self.source_interpolator = Interpolator(expression, V)
+            self.source_interpolate = interpolate(expression, V)
         else:
             self.source_projector = Projector(expression, V)
 
         # If not time-varying, evaluate for the first time here
         if not self.time_varying:
             if self.method == 'interpolate':
-                self.source.assign(self.source_interpolator.interpolate())
+                self.source.assign(assemble(self.source_interpolate))
             else:
                 self.source.assign(self.source_projector.project())
 
@@ -140,7 +141,7 @@ def evaluate(self, x_in, dt):
         if self.time_varying:
             logger.info(f'Evaluating physics parametrisation {self.label.label}')
             if self.method == 'interpolate':
-                self.source.assign(self.source_interpolator.interpolate())
+                self.source.assign(assemble(self.source_interpolate))
             else:
                 self.source.assign(self.source_projector.project())
         else:
diff --git a/gusto/physics/shallow_water_microphysics.py b/gusto/physics/shallow_water_microphysics.py
index bd3f8a02e..7f1b91227 100644
--- a/gusto/physics/shallow_water_microphysics.py
+++ b/gusto/physics/shallow_water_microphysics.py
@@ -3,8 +3,9 @@
 """
 
 from firedrake import (
-    Interpolator, conditional, Function, dx, min_value, max_value, Constant
+    conditional, Function, dx, min_value, max_value, Constant, assemble
 )
+from firedrake.__future__ import interpolate
 from firedrake.fml import subject
 from gusto.core.logging import logger
 from gusto.physics.physics_parametrisation import PhysicsParametrisation
@@ -134,7 +135,7 @@ def __init__(self, equation, saturation_curve,
                                             self.evaluate)
 
         # interpolator does the conversion of vapour to rain
-        self.source_interpolator = Interpolator(conditional(
+        self.source_interpolate = interpolate(conditional(
             self.water_v > self.saturation_curve,
             (1/self.tau)*gamma_r*(self.water_v - self.saturation_curve),
             0), Vv)
@@ -159,7 +160,7 @@ def evaluate(self, x_in, dt):
         if self.set_tau_to_dt:
             self.tau.assign(dt)
         self.water_v.assign(x_in.subfunctions[self.Vv_idx])
-        self.source.assign(self.source_interpolator.interpolate())
+        self.source.assign(assemble(self.source_interpolate))
 
 
 class SWSaturationAdjustment(PhysicsParametrisation):
@@ -321,8 +322,8 @@ def __init__(self, equation, saturation_curve,
         # Add terms to equations and make interpolators
         # sources have the same order as V_idxs and factors
         self.source = [Function(Vc) for factor in factors]
-        self.source_interpolators = [Interpolator(sat_adj_expr*factor, source)
-                                     for factor, source in zip(factors, self.source)]
+        self.source_interpolate = [interpolate(sat_adj_expr*factor, Vc)
+                                   for factor in factors]
 
         # test functions have the same order as factors and sources (vapour,
         # cloud, depth, buoyancy) so that the correct test function multiplies
@@ -359,5 +360,5 @@ def evaluate(self, x_in, dt):
         self.cloud.assign(x_in.subfunctions[self.Vc_idx])
         if self.time_varying_gamma_v:
             self.gamma_v.interpolate(self.gamma_v_computation(x_in))
-        for interpolator in self.source_interpolators:
-            interpolator.interpolate()
+        for interpolator, src in zip(self.source_interpolate, self.source):
+            src.assign(assemble(interpolator))
diff --git a/gusto/recovery/recovery.py b/gusto/recovery/recovery.py
index b88c6ad6a..9d186dd48 100644
--- a/gusto/recovery/recovery.py
+++ b/gusto/recovery/recovery.py
@@ -13,7 +13,8 @@
                        Function, FunctionSpace, Interpolator, Projector,
                        SpatialCoordinate, TensorProductElement,
                        VectorFunctionSpace, as_vector, function, interval,
-                       VectorElement)
+                       VectorElement, assemble)
+from firedrake.__future__ import interpolate
 from gusto.recovery import Averager
 from .recovery_kernels import (BoundaryRecoveryExtruded, BoundaryRecoveryHCurl,
                                BoundaryGaussianElimination)
@@ -144,14 +145,14 @@ def __init__(self, x_inout, method=BoundaryMethod.extruded, eff_coords=None):
             V_broken = FunctionSpace(mesh, BrokenElement(V_inout.ufl_element()))
             self.x_DG1_wrong = Function(V_broken)
             self.x_DG1_correct = Function(V_broken)
-            self.interpolator = Interpolator(self.x_inout, self.x_DG1_wrong)
+            self.interpolate = interpolate(self.x_inout, V_broken)
             self.averager = Averager(self.x_DG1_correct, self.x_inout)
             self.kernel = BoundaryGaussianElimination(V_broken)
 
     def apply(self):
         """Applies the boundary recovery process."""
         if self.method == BoundaryMethod.taylor:
-            self.interpolator.interpolate()
+            self.x_DG1_wrong.assign(assemble(self.interpolate))
             self.kernel.apply(self.x_DG1_wrong, self.x_DG1_correct,
                               self.act_coords, self.eff_coords, self.num_ext)
             self.averager.project()
@@ -275,7 +276,7 @@ def __init__(self, x_in, x_out, method='interpolate', boundary_method=None):
                         self.boundary_recoverers.append(BoundaryRecoverer(x_out_scalars[i],
                                                                           method=BoundaryMethod.taylor,
                                                                           eff_coords=eff_coords[i]))
-                    self.interpolate_to_vector = Interpolator(as_vector(x_out_scalars), self.x_out)
+                    self.interpolate_to_vector = interpolate(as_vector(x_out_scalars), V_out)
 
     def project(self):
         """Perform the whole recovery step."""
@@ -294,7 +295,7 @@ def project(self):
                     # Correct at boundaries
                     boundary_recoverer.apply()
                 # Combine the components to obtain the vector field
-                self.interpolate_to_vector.interpolate()
+                self.x_out.assign(assemble(self.interpolate_to_vector))
             else:
                 # Extrapolate at boundaries
                 self.boundary_recoverer.apply()
diff --git a/gusto/recovery/reversible_recovery.py b/gusto/recovery/reversible_recovery.py
index fc8fbd332..d87449817 100644
--- a/gusto/recovery/reversible_recovery.py
+++ b/gusto/recovery/reversible_recovery.py
@@ -4,7 +4,8 @@
 """
 
 from gusto.core.conservative_projection import ConservativeProjector
-from firedrake import (Projector, Function, Interpolator)
+from firedrake import (Projector, Function, assemble)
+from firedrake.__future__ import interpolate
 from .recovery import Recoverer
 
 __all__ = ["ReversibleRecoverer", "ConservativeRecoverer"]
@@ -52,7 +53,7 @@ def __init__(self, source_field, target_field, reconstruct_opts):
         elif self.opts.project_high_method == 'project':
             self.projector_high = Projector(self.q_recovered, self.q_rec_high)
         elif self.opts.project_high_method == 'interpolate':
-            self.projector_high = Interpolator(self.q_recovered, self.q_rec_high)
+            self.projector_high = interpolate(self.q_recovered, target_field.function_space())
             self.interp_high = True
         else:
             raise ValueError(f'Method {self.opts.project_high_method} '
@@ -68,7 +69,7 @@ def __init__(self, source_field, target_field, reconstruct_opts):
         elif self.opts.project_low_method == 'project':
             self.projector_low = Projector(self.q_rec_high, self.q_corr_low)
         elif self.opts.project_low_method == 'interpolate':
-            self.projector_low = Interpolator(self.q_rec_high, self.q_corr_low)
+            self.projector_low = interpolate(self.q_rec_high, source_field.function_space())
             self.interp_low = True
         else:
             raise ValueError(f'Method {self.opts.project_low_method} '
@@ -84,17 +85,17 @@ def __init__(self, source_field, target_field, reconstruct_opts):
         elif self.opts.injection_method == 'project':
             self.injector = Projector(self.q_corr_low, self.q_corr_high)
         elif self.opts.injection_method == 'interpolate':
-            self.injector = Interpolator(self.q_corr_low, self.q_corr_high)
+            self.injector = interpolate(self.q_corr_low, target_field.function_space())
             self.interp_inj = True
         else:
             raise ValueError(f'Method {self.opts.injection_method} for injection not valid')
 
     def project(self):
         self.recoverer.project()
-        self.projector_high.interpolate() if self.interp_high else self.projector_high.project()
-        self.projector_low.interpolate() if self.interp_low else self.projector_low.project()
+        self.q_rec_high.assign(assemble(self.projector_high)) if self.interp_high else self.projector_high.project()
+        self.q_corr_low.assign(assemble(self.projector_low)) if self.interp_low else self.projector_low.project()
         self.q_corr_low.assign(self.q_low - self.q_corr_low)
-        self.injector.interpolate() if self.interp_inj else self.injector.project()
+        self.q_corr_high.assign(assemble(self.injector)) if self.interp_inj else self.injector.project()
         self.q_high.assign(self.q_corr_high + self.q_rec_high)
 
 
diff --git a/gusto/solvers/linear_solvers.py b/gusto/solvers/linear_solvers.py
index 9aac1df3e..213d09df2 100644
--- a/gusto/solvers/linear_solvers.py
+++ b/gusto/solvers/linear_solvers.py
@@ -11,10 +11,11 @@
     rhs, FacetNormal, div, dx, jump, avg, dS, dS_v, dS_h, ds_v, ds_t, ds_b,
     ds_tb, inner, action, dot, grad, Function, VectorSpaceBasis, cross,
     BrokenElement, FunctionSpace, MixedFunctionSpace, DirichletBC, as_vector,
-    Interpolator, conditional
+    assemble, conditional
 )
 from firedrake.fml import Term, drop
 from firedrake.petsc import flatten_parameters
+from firedrake.__future__ import interpolate
 from pyop2.profiling import timed_function, timed_region
 
 from gusto.equations.active_tracers import TracerVariableType
@@ -660,11 +661,11 @@ def _setup_solver(self):
             qtbar = split(equation.X_ref)[3]
 
             # set up interpolators that use the X_ref values for D and b_e
-            self.q_sat_expr_interpolator = Interpolator(
-                equation.compute_saturation(equation.X_ref), self.q_sat_func)
-            self.q_v_interpolator = Interpolator(
+            self.q_sat_expr_interpolate = interpolate(
+                equation.compute_saturation(equation.X_ref), VD)
+            self.q_v_interpolate = interpolate(
                 conditional(qtbar < self.q_sat_func, qtbar, self.q_sat_func),
-                self.qvbar)
+                VD)
 
             # bbar was be_bar and here we correct to become bbar
             bbar += equation.parameters.beta2 * self.qvbar
@@ -729,8 +730,8 @@ def trace_nullsp(T):
     @timed_function("Gusto:UpdateReferenceProfiles")
     def update_reference_profiles(self):
         if self.equations.equivalent_buoyancy:
-            self.q_sat_expr_interpolator.interpolate()
-            self.q_v_interpolator.interpolate()
+            self.q_sat_func.assign(assemble(self.q_sat_expr_interpolate))
+            self.qvbar.assign(assemble(self.q_v_interpolate))
 
     @timed_function("Gusto:LinearSolve")
     def solve(self, xrhs, dy):
diff --git a/gusto/time_discretisation/explicit_runge_kutta.py b/gusto/time_discretisation/explicit_runge_kutta.py
index b167cf80a..16732bd50 100644
--- a/gusto/time_discretisation/explicit_runge_kutta.py
+++ b/gusto/time_discretisation/explicit_runge_kutta.py
@@ -15,8 +15,8 @@
 
 
 __all__ = [
-    "ForwardEuler", "ExplicitRungeKutta", "SSPRK3", "RK4", "Heun",
-    "RungeKuttaFormulation"
+    "ForwardEuler", "ExplicitRungeKutta", "SSPRK2", "SSPRK3", "SSPRK4",
+    "RK4", "Heun", "RungeKuttaFormulation"
 ]
 
 
@@ -530,21 +530,132 @@ def __init__(
                          augmentation=augmentation)
 
 
+class SSPRK2(ExplicitRungeKutta):
+    u"""
+    Implements 2nd order Strong-Stability-Preserving Runge-Kutta methods
+    for solving ∂y/∂t = F(y).                                                 \n
+
+    Spiteri & Ruuth, 2002, SIAM J. Numer. Anal.                               \n
+    "A new class of optimal high-order strong-stability-preserving time       \n
+    discretisation methods".                                                  \n
+
+    The 2-stage method (Heun's method) can be written as:                     \n
+
+    k0 = F[y^n]                                                               \n
+    k1 = F{y^n + dt*k0]                                                       \n
+    y^(n+1) = y^n + (1/2)*dt*(k0+k1)                                          \n
+
+    The 3-stage method can be written as:                                     \n
+
+    k0 = F[y^n]                                                               \n
+    k1 = F[y^n + (1/2*dt*k0]                                                  \n
+    k2 = F[y^n + (1/2)*dt*(k0+k1)]                                            \n
+    y^(n+1) = y^n + (1/3)*dt*(k0 + k1 + k2)                                   \n
+
+    The 4-stage method can be written as:                                     \n
+
+    k0 = F[y^n]                                                               \n
+    k1 = F[y^n + (1/3)*dt*k1]                                                 \n
+    k2 = F[y^n + (1/3)*dt*(k0+k1)]                                            \n
+    k3 = F[y^n + (1/3)*dt*(k0+k1+k2)]                                         \n
+    y^(n+1) = y^n + (1/4)*dt*(k0 + k1 + k2 + k3)                              \n
+    """
+    def __init__(
+            self, domain, field_name=None, subcycling_options=None,
+            rk_formulation=RungeKuttaFormulation.increment,
+            solver_parameters=None, limiter=None, options=None,
+            augmentation=None, stages=2
+    ):
+        """
+        Args:
+            domain (:class:`Domain`): the model's domain object, containing the
+                mesh and the compatible function spaces.
+            field_name (str, optional): name of the field to be evolved.
+                Defaults to None.
+            subcycling_options(:class:`SubcyclingOptions`, optional): an object
+                containing options for subcycling the time discretisation.
+                Defaults to None.
+            rk_formulation (:class:`RungeKuttaFormulation`, optional):
+                an enumerator object, describing the formulation of the Runge-
+                Kutta scheme. Defaults to the increment form.
+            solver_parameters (dict, optional): dictionary of parameters to
+                pass to the underlying solver. Defaults to None.
+            limiter (:class:`Limiter` object, optional): a limiter to apply to
+                the evolving field to enforce monotonicity. Defaults to None.
+            options (:class:`AdvectionOptions`, optional): an object containing
+                options to either be passed to the spatial discretisation, or
+                to control the "wrapper" methods, such as Embedded DG or a
+                recovery method. Defaults to None.
+            augmentation (:class:`Augmentation`): allows the equation solved in
+                this time discretisation to be augmented, for instances with
+                extra terms of another auxiliary variable. Defaults to None.
+            stages (int, optional): number of stages: (2, 3, 4). Defaults to 2.
+        """
+
+        if stages == 2:
+            butcher_matrix = np.array([
+                [1., 0.],
+                [0.5, 0.5]
+            ])
+            self.cfl_limit = 1
+
+        elif stages == 3:
+            butcher_matrix = np.array([
+                [1./2., 0., 0.],
+                [1./2., 1./2., 0.],
+                [1./3., 1./3., 1./3.]
+            ])
+            self.cfl_limit = 2
+        elif stages == 4:
+            butcher_matrix = np.array([
+                [1./3., 0., 0., 0.],
+                [1./3., 1./3., 0., 0.],
+                [1./3., 1./3., 1./3., 0.],
+                [1./4., 1./4., 1./4., 1./4.]
+            ])
+            self.cfl_limit = 3
+        else:
+            raise ValueError(f"{stages} stage 2rd order SSPRK not implemented")
+
+        super().__init__(domain, butcher_matrix, field_name=field_name,
+                         subcycling_options=subcycling_options,
+                         rk_formulation=rk_formulation,
+                         solver_parameters=solver_parameters,
+                         limiter=limiter, options=options,
+                         augmentation=augmentation)
+
+
 class SSPRK3(ExplicitRungeKutta):
     u"""
-    Implements the 3-stage Strong-Stability-Preserving Runge-Kutta method
-    for solving ∂y/∂t = F(y). It can be written as:                           \n
+    Implements 3rd order Strong-Stability-Preserving Runge-Kutta methods
+    for solving ∂y/∂t = F(y).                                                 \n
+
+    Spiteri & Ruuth, 2002, SIAM J. Numer. Anal.                               \n
+    "A new class of optimal high-order strong-stability-preserving time       \n
+    discretisation methods".                                                  \n
+
+    The 3-stage method can be written as:                                     \n
 
     k0 = F[y^n]                                                               \n
     k1 = F[y^n + dt*k1]                                                       \n
     k2 = F[y^n + (1/4)*dt*(k0+k1)]                                            \n
     y^(n+1) = y^n + (1/6)*dt*(k0 + k1 + 4*k2)                                 \n
+
+    The 4-stage method can be written as:                                     \n
+
+    k0 = F[y^n]                                                               \n
+    k1 = F[y^n + (1/2)*dt*k1]                                                 \n
+    k2 = F[y^n + (1/2)*dt*(k0+k1)]                                            \n
+    k3 = F[y^n + (1/6)*dt*(k0+k1+k2)]                                         \n
+    y^(n+1) = y^n + (1/6)*dt*(k0 + k1 + k2 + 3*k3)                            \n
+
+    The 5-stage method has numerically optimised coefficients.                \n
     """
     def __init__(
             self, domain, field_name=None, subcycling_options=None,
             rk_formulation=RungeKuttaFormulation.increment,
             solver_parameters=None, limiter=None, options=None,
-            augmentation=None
+            augmentation=None, stages=3
     ):
         """
         Args:
@@ -569,13 +680,36 @@ def __init__(
             augmentation (:class:`Augmentation`): allows the equation solved in
                 this time discretisation to be augmented, for instances with
                 extra terms of another auxiliary variable. Defaults to None.
+            stages (int, optional): number of stages: (3, 4, 5). Defaults to 3.
         """
 
-        butcher_matrix = np.array([
-            [1., 0., 0.],
-            [1./4., 1./4., 0.],
-            [1./6., 1./6., 2./3.]
-        ])
+        if stages == 3:
+            butcher_matrix = np.array([
+                [1., 0., 0.],
+                [1./4., 1./4., 0.],
+                [1./6., 1./6., 2./3.]
+            ])
+            self.cfl_limit = 1
+        elif stages == 4:
+            butcher_matrix = np.array([
+                [1./2., 0., 0., 0.],
+                [1./2., 1./2., 0., 0.],
+                [1./6., 1./6., 1./6., 0.],
+                [1./6., 1./6., 1./6., 1./2.]
+            ])
+            self.cfl_limit = 2
+        elif stages == 5:
+            self.cfl_limit = 2.65062919294483
+            butcher_matrix = np.array([
+                [0.37726891511710, 0., 0., 0., 0.],
+                [0.37726891511710, 0.37726891511710, 0., 0., 0.],
+                [0.16352294089771, 0.16352294089771, 0.16352294089771, 0., 0.],
+                [0.14904059394856, 0.14831273384724, 0.14831273384724, 0.34217696850008, 0.],
+                [0.19707596384481, 0.11780316509765, 0.11709725193772, 0.27015874934251, 0.29786487010104]
+            ])
+        else:
+            raise ValueError(f"{stages} stage 3rd order SSPRK not implemented")
+
         super().__init__(domain, butcher_matrix, field_name=field_name,
                          subcycling_options=subcycling_options,
                          rk_formulation=rk_formulation,
@@ -584,26 +718,22 @@ def __init__(
                          augmentation=augmentation)
 
 
-class RK4(ExplicitRungeKutta):
+class SSPRK4(ExplicitRungeKutta):
     u"""
-    Implements the classic 4-stage Runge-Kutta method.
+    Implements 4th order Strong-Stability-Preserving Runge-Kutta methods
+    for solving ∂y/∂t = F(y).                                                 \n
 
-    The classic 4-stage Runge-Kutta method for solving ∂y/∂t = F(y). It can be
-    written as:                                                               \n
+    Spiteri & Ruuth, 2002, SIAM J. Numer. Anal.                               \n
+    "A new class of optimal high-order strong-stability-preserving time       \n
+    discretisation methods".                                                  \n
 
-    k0 = F[y^n]                                                               \n
-    k1 = F[y^n + 1/2*dt*k1]                                                   \n
-    k2 = F[y^n + 1/2*dt*k2]                                                   \n
-    k3 = F[y^n + dt*k3]                                                       \n
-    y^(n+1) = y^n + (1/6) * dt * (k0 + 2*k1 + 2*k2 + k3)                      \n
-
-    where superscripts indicate the time-level.                               \n
+    The 5-stage method has numerically optimised coefficients.                \n
     """
     def __init__(
             self, domain, field_name=None, subcycling_options=None,
             rk_formulation=RungeKuttaFormulation.increment,
             solver_parameters=None, limiter=None, options=None,
-            augmentation=None
+            augmentation=None, stages=5
     ):
         """
         Args:
@@ -628,13 +758,21 @@ def __init__(
             augmentation (:class:`Augmentation`): allows the equation solved in
                 this time discretisation to be augmented, for instances with
                 extra terms of another auxiliary variable. Defaults to None.
+            stages (int, optional): number of stages: (5,). Defaults to 5.
         """
-        butcher_matrix = np.array([
-            [0.5, 0., 0., 0.],
-            [0., 0.5, 0., 0.],
-            [0., 0., 1., 0.],
-            [1./6., 1./3., 1./3., 1./6.]
-        ])
+
+        if stages == 5:
+            self.cfl_limit = 1.50818004975927
+            butcher_matrix = np.array([
+                [0.39175222700392, 0., 0., 0., 0.],
+                [0.21766909633821, 0.36841059262959, 0., 0., 0.],
+                [0.08269208670950, 0.13995850206999, 0.25189177424738, 0., 0.],
+                [0.06796628370320, 0.11503469844438, 0.20703489864929, 0.54497475021237, 0.],
+                [0.14681187618661, 0.24848290924556, 0.10425883036650, 0.27443890091960, 0.22600748319395]
+            ])
+        else:
+            raise ValueError(f"{stages} stage 4rd order SSPRK not implemented")
+
         super().__init__(domain, butcher_matrix, field_name=field_name,
                          subcycling_options=subcycling_options,
                          rk_formulation=rk_formulation,
@@ -643,18 +781,20 @@ def __init__(
                          augmentation=augmentation)
 
 
-class Heun(ExplicitRungeKutta):
+class RK4(ExplicitRungeKutta):
     u"""
-    Implements Heun's method.
+    Implements the classic 4-stage Runge-Kutta method.
 
-    The 2-stage Runge-Kutta scheme known as Heun's method,for solving
-    ∂y/∂t = F(y). It can be written as:                                       \n
+    The classic 4-stage Runge-Kutta method for solving ∂y/∂t = F(y). It can be
+    written as:                                                               \n
 
-    y_1 = F[y^n]                                                              \n
-    y^(n+1) = (1/2)y^n + (1/2)F[y_1]                                          \n
+    k0 = F[y^n]                                                               \n
+    k1 = F[y^n + 1/2*dt*k1]                                                   \n
+    k2 = F[y^n + 1/2*dt*k2]                                                   \n
+    k3 = F[y^n + dt*k3]                                                       \n
+    y^(n+1) = y^n + (1/6) * dt * (k0 + 2*k1 + 2*k2 + k3)                      \n
 
-    where superscripts indicate the time-level and subscripts indicate the stage
-    number.
+    where superscripts indicate the time-level.                               \n
     """
     def __init__(
             self, domain, field_name=None, subcycling_options=None,
@@ -686,10 +826,11 @@ def __init__(
                 this time discretisation to be augmented, for instances with
                 extra terms of another auxiliary variable. Defaults to None.
         """
-
         butcher_matrix = np.array([
-            [1., 0.],
-            [0.5, 0.5]
+            [0.5, 0., 0., 0.],
+            [0., 0.5, 0., 0.],
+            [0., 0., 1., 0.],
+            [1./6., 1./3., 1./3., 1./6.]
         ])
         super().__init__(domain, butcher_matrix, field_name=field_name,
                          subcycling_options=subcycling_options,
@@ -697,3 +838,24 @@ def __init__(
                          solver_parameters=solver_parameters,
                          limiter=limiter, options=options,
                          augmentation=augmentation)
+
+
+def Heun(domain, field_name=None, subcycling_options=None,
+         rk_formulation=RungeKuttaFormulation.increment,
+         solver_parameters=None, limiter=None, options=None,
+         augmentation=None):
+    u"""
+    Implements Heun's method.
+
+    The 2-stage Runge-Kutta scheme known as Heun's method,for solving
+    ∂y/∂t = F(y). It can be written as:                                       \n
+
+    y_1 = F[y^n]                                                              \n
+    y^(n+1) = (1/2)y^n + (1/2)F[y_1]                                          \n
+
+    where superscripts indicate the time-level and subscripts indicate the stage
+    number.
+    """
+    return SSPRK2(domain, field_name=field_name, subcycling_options=subcycling_options,
+                  rk_formulation=rk_formulation, solver_parameters=solver_parameters,
+                  limiter=limiter, options=options, augmentation=augmentation, stages=2)
diff --git a/gusto/time_discretisation/wrappers.py b/gusto/time_discretisation/wrappers.py
index f0aae843e..18381bc1f 100644
--- a/gusto/time_discretisation/wrappers.py
+++ b/gusto/time_discretisation/wrappers.py
@@ -6,9 +6,10 @@
 
 from abc import ABCMeta, abstractmethod
 from firedrake import (
-    FunctionSpace, Function, BrokenElement, Projector, Interpolator,
-    Constant, as_ufl, dot, grad, TestFunction, MixedFunctionSpace
+    FunctionSpace, Function, BrokenElement, Projector, Constant, as_ufl, dot, grad,
+    TestFunction, MixedFunctionSpace, assemble
 )
+from firedrake.__future__ import interpolate
 from firedrake.fml import Term
 from gusto.core.configuration import EmbeddedDGOptions, RecoveryOptions, SUPGOptions
 from gusto.recovery import Recoverer, ReversibleRecoverer, ConservativeRecoverer
@@ -264,7 +265,7 @@ def setup(self, original_space, post_apply_bcs):
         # Operators for projecting back
         self.interp_back = (self.options.project_low_method == 'interpolate')
         if self.options.project_low_method == 'interpolate':
-            self.x_out_projector = Interpolator(self.x_out, self.x_projected)
+            self.x_out_projector = interpolate(self.x_out, self.original_space)
         elif self.options.project_low_method == 'project':
             self.x_out_projector = Projector(self.x_out, self.x_projected,
                                              bcs=post_apply_bcs)
@@ -302,7 +303,7 @@ def post_apply(self, x_out):
         """
 
         if self.interp_back:
-            self.x_out_projector.interpolate()
+            self.x_projected.assign(assemble(self.x_out_projector))
         else:
             self.x_out_projector.project()
         x_out.assign(self.x_projected)
diff --git a/gusto/timestepping/semi_implicit_quasi_newton.py b/gusto/timestepping/semi_implicit_quasi_newton.py
index 7920cb65c..deacb6af2 100644
--- a/gusto/timestepping/semi_implicit_quasi_newton.py
+++ b/gusto/timestepping/semi_implicit_quasi_newton.py
@@ -4,10 +4,11 @@
 """
 
 from firedrake import (
-    Function, Constant, TrialFunctions, DirichletBC, div, Interpolator,
+    Function, Constant, TrialFunctions, DirichletBC, div, assemble,
     LinearVariationalProblem, LinearVariationalSolver
 )
 from firedrake.fml import drop, replace_subject
+from firedrake.__future__ import interpolate
 from pyop2.profiling import timed_stage
 from gusto.core import TimeLevelFields, StateFields
 from gusto.core.labels import (transport, diffusion, time_derivative,
@@ -244,9 +245,8 @@ def __init__(self, equation_set, io, transport_schemes, spatial_methods,
             V_DG = equation_set.domain.spaces('DG')
             self.predictor_field_in = Function(V_DG)
             div_factor = Constant(1.0) - (Constant(1.0) - self.alpha)*self.dt*div(self.x.n('u'))
-            self.predictor_interpolator = Interpolator(
-                self.x.star(predictor)*div_factor, self.predictor_field_in
-            )
+            self.predictor_interpolate = interpolate(
+                self.x.star(predictor)*div_factor, V_DG)
 
     def _apply_bcs(self):
         """
@@ -344,7 +344,7 @@ def transport_fields(self, outer, xstar, xp):
                     # Pre-multiply this variable by (1 - dt*beta*div(u))
                     V = xstar(name).function_space()
                     field_out = Function(V)
-                    self.predictor_interpolator.interpolate()
+                    self.predictor_field_in.assign(assemble(self.predictor_interpolate))
                     scheme.apply(field_out, self.predictor_field_in)
 
                     # xp is xstar plus the increment from the transported predictor
diff --git a/integration-tests/model/test_time_discretisation.py b/integration-tests/model/test_time_discretisation.py
index 6d484bfd1..0c259e221 100644
--- a/integration-tests/model/test_time_discretisation.py
+++ b/integration-tests/model/test_time_discretisation.py
@@ -10,9 +10,19 @@ def run(timestepper, tmax, f_end):
 
 @pytest.mark.parametrize(
     "scheme", [
-        "ssprk3_increment", "TrapeziumRule", "ImplicitMidpoint", "QinZhang",
+        "ssprk2_increment_2", "ssprk2_predictor_2", "ssprk2_linear_2",
+        "ssprk2_increment_3", "ssprk2_predictor_3", "ssprk2_linear_3",
+        "ssprk2_increment_4", "ssprk2_predictor_4", "ssprk2_linear_4",
+
+        "ssprk3_increment_3", "ssprk3_predictor_3", "ssprk3_linear_3",
+        "ssprk3_increment_4", "ssprk3_predictor_4", "ssprk3_linear_4",
+        "ssprk3_increment_5", "ssprk3_predictor_5", "ssprk3_linear_5",
+
+        "ssprk4_increment_5", "ssprk4_predictor_5", "ssprk4_linear_5",
+
+        "TrapeziumRule", "ImplicitMidpoint", "QinZhang",
         "RK4", "Heun", "BDF2", "TR_BDF2", "AdamsBashforth", "Leapfrog",
-        "AdamsMoulton", "AdamsMoulton", "ssprk3_predictor", "ssprk3_linear"
+        "AdamsMoulton", "AdamsMoulton"
     ]
 )
 def test_time_discretisation(tmpdir, scheme, tracer_setup):
@@ -30,12 +40,52 @@ def test_time_discretisation(tmpdir, scheme, tracer_setup):
         V = domain.spaces("DG")
         eqn = AdvectionEquation(domain, V, "f")
 
-    if scheme == "ssprk3_increment":
+    if scheme == "ssprk2_increment_2":
+        transport_scheme = SSPRK2(domain, rk_formulation=RungeKuttaFormulation.increment)
+    elif scheme == "ssprk2_predictor_2":
+        transport_scheme = SSPRK2(domain, rk_formulation=RungeKuttaFormulation.predictor)
+    elif scheme == "ssprk2_linear_2":
+        transport_scheme = SSPRK2(domain, rk_formulation=RungeKuttaFormulation.linear)
+    elif scheme == "ssprk2_increment_3":
+        transport_scheme = SSPRK2(domain, rk_formulation=RungeKuttaFormulation.increment, stages=3)
+    elif scheme == "ssprk2_predictor_3":
+        transport_scheme = SSPRK2(domain, rk_formulation=RungeKuttaFormulation.predictor, stages=3)
+    elif scheme == "ssprk2_linear_3":
+        transport_scheme = SSPRK2(domain, rk_formulation=RungeKuttaFormulation.linear, stages=3)
+    elif scheme == "ssprk2_increment_4":
+        transport_scheme = SSPRK2(domain, rk_formulation=RungeKuttaFormulation.increment, stages=4)
+    elif scheme == "ssprk2_predictor_4":
+        transport_scheme = SSPRK2(domain, rk_formulation=RungeKuttaFormulation.predictor, stages=4)
+    elif scheme == "ssprk2_linear_4":
+        transport_scheme = SSPRK2(domain, rk_formulation=RungeKuttaFormulation.linear, stages=4)
+
+    elif scheme == "ssprk3_increment_3":
         transport_scheme = SSPRK3(domain, rk_formulation=RungeKuttaFormulation.increment)
-    elif scheme == "ssprk3_predictor":
+    elif scheme == "ssprk3_predictor_3":
         transport_scheme = SSPRK3(domain, rk_formulation=RungeKuttaFormulation.predictor)
-    elif scheme == "ssprk3_linear":
+    elif scheme == "ssprk3_linear_3":
         transport_scheme = SSPRK3(domain, rk_formulation=RungeKuttaFormulation.linear)
+    elif scheme == "ssprk3_increment_4":
+        transport_scheme = SSPRK3(domain, rk_formulation=RungeKuttaFormulation.increment, stages=4)
+    elif scheme == "ssprk3_predictor_4":
+        transport_scheme = SSPRK3(domain, rk_formulation=RungeKuttaFormulation.predictor, stages=4)
+    elif scheme == "ssprk3_linear_4":
+        transport_scheme = SSPRK3(domain, rk_formulation=RungeKuttaFormulation.linear, stages=4)
+
+    elif scheme == "ssprk3_increment_5":
+        transport_scheme = SSPRK3(domain, rk_formulation=RungeKuttaFormulation.increment, stages=5)
+    elif scheme == "ssprk3_predictor_5":
+        transport_scheme = SSPRK3(domain, rk_formulation=RungeKuttaFormulation.predictor, stages=5)
+    elif scheme == "ssprk3_linear_5":
+        transport_scheme = SSPRK3(domain, rk_formulation=RungeKuttaFormulation.linear, stages=5)
+
+    elif scheme == "ssprk4_increment_5":
+        transport_scheme = SSPRK4(domain, rk_formulation=RungeKuttaFormulation.increment)
+    elif scheme == "ssprk4_predictor_5":
+        transport_scheme = SSPRK4(domain, rk_formulation=RungeKuttaFormulation.predictor)
+    elif scheme == "ssprk4_linear_5":
+        transport_scheme = SSPRK4(domain, rk_formulation=RungeKuttaFormulation.linear)
+
     elif scheme == "TrapeziumRule":
         transport_scheme = TrapeziumRule(domain)
     elif scheme == "ImplicitMidpoint":
diff --git a/integration-tests/physics/test_held_suarez_friction.py b/integration-tests/physics/test_held_suarez_friction.py
new file mode 100644
index 000000000..dfe7e34f9
--- /dev/null
+++ b/integration-tests/physics/test_held_suarez_friction.py
@@ -0,0 +1,114 @@
+"""
+This tests the Rayleigh friction term used in the Held Suarez test case.
+"""
+
+from gusto import *
+import gusto.equations.thermodynamics as td
+from gusto.core.labels import physics_label
+from firedrake import (Constant, PeriodicIntervalMesh, as_vector, norm,
+                       ExtrudedMesh, Function, dot)
+from firedrake.fml import identity, drop
+
+
+def run_apply_rayleigh_friction(dirname):
+    # ------------------------------------------------------------------------ #
+    # Set up model objects
+    # ------------------------------------------------------------------------ #
+
+    dt = 3600.0
+
+    # declare grid shape, with length L and height H
+    L = 500.
+    H = 500.
+    nlayers = int(H / 5.)
+    ncolumns = int(L / 5.)
+
+    # make mesh and domain
+    m = PeriodicIntervalMesh(ncolumns, L)
+    mesh = ExtrudedMesh(m, layers=nlayers, layer_height=(H / nlayers))
+    domain = Domain(mesh, dt, "CG", 0)
+
+    # Set up equation
+    parameters = CompressibleParameters()
+    eqn = CompressibleEulerEquations(domain, parameters)
+
+    # I/O
+    output = OutputParameters(dirname=dirname+"/held_suarez_friction",
+                              dumpfreq=1,
+                              dumplist=['u'])
+    io = IO(domain, output)
+
+    # Physics scheme
+    physics_parametrisation = RayleighFriction(eqn)
+
+    time_discretisation = BackwardEuler(domain)
+
+    # time_discretisation = ForwardEuler(domain)
+    physics_schemes = [(physics_parametrisation, time_discretisation)]
+
+    # Only want time derivatives and physics terms in equation, so drop the rest
+    eqn.residual = eqn.residual.label_map(lambda t: any(t.has_label(time_derivative, physics_label)),
+                                          map_if_true=identity, map_if_false=drop)
+
+    # Time stepper
+    scheme = ForwardEuler(domain)
+    stepper = SplitPhysicsTimestepper(eqn, scheme, io,
+                                      physics_schemes=physics_schemes)
+
+    # ------------------------------------------------------------------------ #
+    # Initial conditions
+    # ------------------------------------------------------------------------ #
+
+    Vu = domain.spaces("HDiv")
+    Vt = domain.spaces("theta")
+    Vr = domain.spaces("DG")
+
+    # Declare prognostic fields
+    u0 = stepper.fields("u")
+    rho0 = stepper.fields("rho")
+    theta0 = stepper.fields("theta")
+
+    # Set a background state with constant pressure and temperature
+    pressure = Function(Vr).interpolate(Constant(100000.))
+    temperature = Function(Vt).interpolate(Constant(295.))
+    theta_d = td.theta(parameters, temperature, pressure)
+
+    theta0.project(theta_d)
+    rho0.interpolate(pressure / (temperature*parameters.R_d))
+
+    # Constant horizontal wind
+    u0.project(as_vector([864, 0.0]))
+
+    # Answer: slower winds than initially
+    u_true = Function(Vu)
+    u_true.project(as_vector([828, 0.0]))
+
+    # ------------------------------------------------------------------------ #
+    # Run
+    # ------------------------------------------------------------------------ #
+
+    stepper.run(t=0, tmax=dt)
+
+    return mesh, stepper, u_true
+
+
+def test_rayleigh_friction(tmpdir):
+
+    dirname = str(tmpdir)
+    mesh, stepper, u_true = run_apply_rayleigh_friction(dirname)
+
+    u_final = stepper.fields('u')
+
+    # Project into CG1 to get sensible values
+    e_x = as_vector([1.0, 0.0])
+    e_z = as_vector([0.0, 1.0])
+
+    DG0 = FunctionSpace(mesh, "DG", 0)
+    u_x_final = Function(DG0).project(dot(u_final, e_x))
+    u_x_true = Function(DG0).project(dot(u_true, e_x))
+    u_z_final = Function(DG0).project(dot(u_final, e_z))
+    u_z_true = Function(DG0).project(dot(u_true, e_z))
+
+    denom = norm(u_x_true)
+    assert norm(u_x_final - u_x_true) / denom < 0.0001, 'Final horizontal wind is incorrect'
+    assert norm(u_z_final - u_z_true) < 1e-12, 'Final vertical wind is incorrect'
diff --git a/integration-tests/physics/test_held_suarez_relaxation.py b/integration-tests/physics/test_held_suarez_relaxation.py
new file mode 100644
index 000000000..615ba8d29
--- /dev/null
+++ b/integration-tests/physics/test_held_suarez_relaxation.py
@@ -0,0 +1,107 @@
+"""
+This tests the Held-Suarez physics routine to apply Rayleigh friction.
+"""
+from gusto import *
+import gusto.equations.thermodynamics as td
+from gusto.core.labels import physics_label
+from firedrake import (Constant, PeriodicIntervalMesh, as_vector,
+                       ExtrudedMesh, Function)
+from firedrake.fml import identity, drop
+import pytest
+
+
+def run_held_suarez_relaxation(dirname, temp):
+
+    # ------------------------------------------------------------------------ #
+    # Set up model objects
+    # ------------------------------------------------------------------------ #
+
+    dt = 3600.0
+
+    # declare grid shape, with length L and height H
+    L = 500.
+    H = 500.
+    nlayers = int(H / 5.)
+    ncolumns = int(L / 5.)
+
+    # make mesh and domain
+    m = PeriodicIntervalMesh(ncolumns, L)
+    mesh = ExtrudedMesh(m, layers=nlayers, layer_height=(H / nlayers))
+    domain = Domain(mesh, dt, "CG", 0)
+
+    # Set up equation
+    parameters = CompressibleParameters()
+    eqn = CompressibleEulerEquations(domain, parameters)
+
+    # I/O
+    output = OutputParameters(dirname=dirname+"/held_suarez_friction",
+                              dumpfreq=1,
+                              dumplist=['u'])
+    io = IO(domain, output)
+
+    # Physics scheme
+    physics_parametrisation = Relaxation(eqn, 'theta', parameters)
+
+    time_discretisation = ForwardEuler(domain)
+
+    # time_discretisation = ForwardEuler(domain)
+    physics_schemes = [(physics_parametrisation, time_discretisation)]
+
+    # Only want time derivatives and physics terms in equation, so drop the rest
+    eqn.residual = eqn.residual.label_map(lambda t: any(t.has_label(time_derivative, physics_label)),
+                                          map_if_true=identity, map_if_false=drop)
+
+    # Time stepper
+    scheme = ForwardEuler(domain)
+    stepper = SplitPhysicsTimestepper(eqn, scheme, io,
+                                      physics_schemes=physics_schemes)
+
+    # ------------------------------------------------------------------------ #
+    # Initial conditions
+    # ------------------------------------------------------------------------ #
+
+    Vt = domain.spaces("theta")
+    Vr = domain.spaces("DG")
+
+    # Declare prognostic fields
+    u0 = stepper.fields("u")
+    rho0 = stepper.fields("rho")
+    theta0 = stepper.fields("theta")
+
+    # Set a background state with constant pressure and temperature
+    p0 = 100000.
+    pressure = Function(Vr).interpolate(Constant(p0))
+    temperature = Function(Vt).interpolate(Constant(temp))
+    theta_d = td.theta(parameters, temperature, pressure)
+
+    theta0.project(theta_d)
+    rho0.interpolate(p0 / (theta_d*parameters.R_d))  # This ensures that exner = 1
+
+    # Constant horizontal wind
+    u0.project(as_vector([1.0, 1.0]))
+    theta_initial = Function(Vt)
+    theta_initial.interpolate(theta_d)
+
+    # ------------------------------------------------------------------------ #
+    # Run
+    # ------------------------------------------------------------------------ #
+
+    stepper.run(t=0, tmax=dt)
+
+    return stepper, theta_initial
+
+
+@pytest.mark.parametrize('temp', [280, 290])
+def test_held_suarez_relaxation(tmpdir, temp):
+    # By configuring the fields we have set the equilibrium temperature to 285K
+    # We test a temperature value eith side to check it moves in the right direction
+    dirname = str(tmpdir)
+    stepper, theta_initial = run_held_suarez_relaxation(dirname, temp)
+
+    theta_final = stepper.fields('theta')
+    final_data = theta_final.dat.data
+    initial_data = theta_initial.dat.data
+    if temp == 280:
+        assert np.mean(final_data) > np.mean(initial_data)
+    if temp == 290:
+        assert np.mean(final_data) < np.mean(initial_data)