diff --git a/gusto/active_tracers.py b/gusto/active_tracers.py
index ab0d59b28..60329e5fb 100644
--- a/gusto/active_tracers.py
+++ b/gusto/active_tracers.py
@@ -53,6 +53,7 @@ class ActiveTracer(object):
     """
     def __init__(self, name, space, variable_type,
                  transport_eqn=TransportEquationType.advective,
+                 density_name=None,
                  phase=Phases.gas, chemical=None):
         """
         Args:
@@ -63,6 +64,10 @@ def __init__(self, name, space, variable_type,
             transport_eqn (:class:`TransportEquationType`, optional): enumerator
                 indicating the type of transport equation to be solved (e.g.
                 advective). Defaults to `TransportEquationType.advective`.
+            density_name (str): the name of the associated density for a mixing
+                ratio when using the tracer_conservative transport. Defaults to None,
+                but raises an error if tracer_conservative transport is used without
+                a specified density.
             phase (:class:`Phases`, optional): enumerator indicating the phase
                 of the tracer variable. Defaults to `Phases.gas`.
             chemical (str, optional): string to describe the chemical that this
@@ -75,14 +80,21 @@ def __init__(self, name, space, variable_type,
         self.name = name
         self.space = space
         self.transport_eqn = transport_eqn
+        self.density_name = density_name
         self.variable_type = variable_type
         self.phase = phase
         self.chemical = chemical
 
-        if (variable_type == TracerVariableType.density and transport_eqn == TransportEquationType.advective):
+        if (variable_type == TracerVariableType.density
+                and transport_eqn == TransportEquationType.advective):
             logger.warning('Active tracer initialised which describes a '
                            + 'density but solving the advective transport eqn')
 
+        if (transport_eqn == TransportEquationType.tracer_conservative
+                and density_name is None):
+            raise ValueError(f'Active tracer {name} using tracer conservative '
+                             + 'transport needs an associated density.')
+
 
 class WaterVapour(ActiveTracer):
     """An object encoding the details of water vapour as a tracer."""
diff --git a/gusto/common_forms.py b/gusto/common_forms.py
index 739485cd8..beedf20bd 100644
--- a/gusto/common_forms.py
+++ b/gusto/common_forms.py
@@ -12,7 +12,7 @@
 __all__ = ["advection_form", "continuity_form", "vector_invariant_form",
            "kinetic_energy_form", "advection_equation_circulation_form",
            "diffusion_form", "linear_advection_form", "linear_continuity_form",
-           "split_continuity_form"]
+           "split_continuity_form", "tracer_conservative_form"]
 
 
 def advection_form(test, q, ubar):
@@ -262,3 +262,28 @@ def split_continuity_form(equation):
                 map_if_true=drop)
 
     return equation
+
+
+def tracer_conservative_form(test, q, rho, ubar):
+    u"""
+    The form corresponding to the continuity transport operator.
+
+    This describes ∇.(u*q*rho) for transporting velocity u and a
+    transported tracer (mixing ratio), q, with an associated density, rho.
+
+    Args:
+        test (:class:`TestFunction`): the test function.
+        q (:class:`ufl.Expr`): the tracer to be transported.
+        rho (:class:`ufl.Expr`): the reference density that will
+        mulitply with q before taking the divergence.
+        ubar (:class:`ufl.Expr`): the transporting velocity.
+
+    Returns:
+        class:`LabelledForm`: a labelled transport form.
+    """
+
+    q_rho = q*rho
+    L = inner(test, div(outer(q_rho, ubar)))*dx
+    form = transporting_velocity(L, ubar)
+
+    return transport(form, TransportEquationType.tracer_conservative)
diff --git a/gusto/configuration.py b/gusto/configuration.py
index 5c909aa59..a2b9b6dd2 100644
--- a/gusto/configuration.py
+++ b/gusto/configuration.py
@@ -31,6 +31,8 @@ class TransportEquationType(Enum):
     conservative: ∂q/∂t + ∇.(u*q) = 0                                         \n
     vector_invariant: ∂q/∂t + (∇×q)×u + (1/2)∇(q.u) + (1/2)[(∇q).u -(∇u).q)] = 0
     circulation: ∂q/∂t + (∇×q)×u + non-transport terms = 0
+    tracer_conservative: ∂(q*rho)/∂t + ∇.(u*q*rho) = 0, for a reference density of rho
+    for the tracer, q.
     """
 
     no_transport = 702
@@ -38,6 +40,7 @@ class TransportEquationType(Enum):
     conservative = 291
     vector_invariant = 9081
     circulation = 512
+    tracer_conservative = 296
 
 
 class Configuration(object):
diff --git a/gusto/equations.py b/gusto/equations.py
index 73953053e..b98181d30 100644
--- a/gusto/equations.py
+++ b/gusto/equations.py
@@ -20,7 +20,8 @@
 from gusto.common_forms import (
     advection_form, continuity_form, vector_invariant_form,
     kinetic_energy_form, advection_equation_circulation_form,
-    diffusion_form, linear_continuity_form, linear_advection_form
+    diffusion_form, linear_continuity_form, linear_advection_form,
+    tracer_conservative_form
 )
 from gusto.active_tracers import ActiveTracer, Phases, TracerVariableType
 from gusto.configuration import TransportEquationType
@@ -430,6 +431,43 @@ def add_tracers_to_prognostics(self, domain, active_tracers):
             else:
                 raise TypeError(f'Tracers must be ActiveTracer objects, not {type(tracer)}')
 
+    def generate_tracer_mass_terms(self, active_tracers):
+        """
+        Adds the mass forms for the active tracers to the equation set.
+
+        Args:
+            active_tracers (list): A list of :class:`ActiveTracer` objects that
+                encode the metadata for the active tracers.
+
+        Returns:
+            :class:`LabelledForm`: a labelled form containing the mass
+                terms for the active tracers. This is the usual mass form
+                unless using tracer_conservative, where it is multiplied
+                by the reference density.
+        """
+
+        for i, tracer in enumerate(active_tracers):
+            idx = self.field_names.index(tracer.name)
+            tracer_prog = split(self.X)[idx]
+            tracer_test = self.tests[idx]
+
+            if tracer.transport_eqn == TransportEquationType.tracer_conservative:
+                ref_density_idx = self.field_names.index(tracer.density_name)
+                ref_density = split(self.X)[ref_density_idx]
+                q = tracer_prog*ref_density
+                mass = subject(prognostic(inner(q, tracer_test)*dx,
+                                          self.field_names[idx]), self.X)
+            else:
+                mass = subject(prognostic(inner(tracer_prog, tracer_test)*dx,
+                                          self.field_names[idx]), self.X)
+
+            if i == 0:
+                mass_form = time_derivative(mass)
+            else:
+                mass_form += time_derivative(mass)
+
+        return mass_form
+
     def generate_tracer_transport_terms(self, active_tracers):
         """
         Adds the transport forms for the active tracers to the equation set.
@@ -473,6 +511,13 @@ def generate_tracer_transport_terms(self, active_tracers):
                     tracer_adv = prognostic(
                         continuity_form(tracer_test, tracer_prog, u),
                         tracer.name)
+                elif tracer.transport_eqn == TransportEquationType.tracer_conservative:
+                    ref_density_idx = self.field_names.index(tracer.density_name)
+                    ref_density = split(self.X)[ref_density_idx]
+                    tracer_adv = prognostic(
+                        tracer_conservative_form(tracer_test, tracer_prog,
+                                                 ref_density, u), tracer.name)
+
                 else:
                     raise ValueError(f'Transport eqn {tracer.transport_eqn} not recognised')
 
@@ -562,9 +607,9 @@ def __init__(self, domain, active_tracers, Vu=None):
         self.tests = TestFunctions(W)
         self.X = Function(W)
 
-        mass_form = self.generate_mass_terms()
-
-        self.residual = subject(mass_form, self.X)
+        # Add mass forms for the tracers, which will use
+        # mass*density for any tracer_conservative terms
+        self.residual = self.generate_tracer_mass_terms(active_tracers)
 
         # Add transport of tracers
         self.residual += self.generate_tracer_transport_terms(active_tracers)
@@ -574,7 +619,6 @@ def __init__(self, domain, active_tracers, Vu=None):
 # Specified Equation Sets
 # ============================================================================ #
 
-
 class ShallowWaterEquations(PrognosticEquationSet):
     u"""
     Class for the (rotating) shallow-water equations, which evolve the velocity
diff --git a/gusto/time_discretisation.py b/gusto/time_discretisation.py
index 7efec29e7..a1221c6b9 100644
--- a/gusto/time_discretisation.py
+++ b/gusto/time_discretisation.py
@@ -14,7 +14,7 @@
     NonlinearVariationalSolver, DirichletBC, split, Constant
 )
 from firedrake.fml import (
-    replace_subject, replace_test_function, Term, all_terms, drop
+    replace_subject, replace_test_function, Term, all_terms, drop, keep
 )
 from firedrake.formmanipulation import split_form
 from firedrake.utils import cached_property
@@ -794,8 +794,9 @@ class ExplicitMultistage(ExplicitTimeDiscretisation):
     #  [ b_0  b_1  .       b_s]
     # ---------------------------------------------------------------------------
 
-    def __init__(self, domain, butcher_matrix, field_name=None, fixed_subcycles=None,
-                 subcycle_by_courant=None, solver_parameters=None,
+    def __init__(self, domain, butcher_matrix, field_name=None,
+                 fixed_subcycles=None, subcycle_by_courant=None,
+                 increment_form=True, solver_parameters=None,
                  limiter=None, options=None):
         """
         Args:
@@ -814,6 +815,9 @@ def __init__(self, domain, butcher_matrix, field_name=None, fixed_subcycles=None
                 for one sub-cycle. Defaults to None, in which case adaptive
                 sub-cycling is not used. This option cannot be specified with the
                 `fixed_subcycles` argument.
+            increment_form (bool, optional): whether to write the RK scheme in
+                "increment form", solving for increments rather than updated
+                fields. Defaults to True.
             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
@@ -830,6 +834,7 @@ def __init__(self, domain, butcher_matrix, field_name=None, fixed_subcycles=None
                          limiter=limiter, options=options)
         self.butcher_matrix = butcher_matrix
         self.nbutcher = int(np.shape(self.butcher_matrix)[0])
+        self.increment_form = increment_form
 
     @property
     def nStages(self):
@@ -846,58 +851,161 @@ def setup(self, equation, apply_bcs=True, *active_labels):
         """
         super().setup(equation, apply_bcs, *active_labels)
 
-        self.k = [Function(self.fs) for i in range(self.nStages)]
+        if not self.increment_form:
+            self.field_i = [Function(self.fs) for i in range(self.nStages+1)]
+        else:
+            self.k = [Function(self.fs) for i in range(self.nStages)]
+
+    @cached_property
+    def solver(self):
+        if self.increment_form:
+            return super().solver
+        else:
+            # In this case, don't set snes_type to ksp only, as we do want the
+            # outer Newton iteration
+            solver_list = []
+
+            for stage in range(self.nStages):
+                # setup linear solver using lhs and rhs defined in derived class
+                problem = NonlinearVariationalProblem(
+                    self.lhs[stage].form - self.rhs[stage].form,
+                    self.field_i[stage+1], bcs=self.bcs)
+                solver_name = self.field_name+self.__class__.__name__+str(stage)
+                solver = NonlinearVariationalSolver(
+                    problem, solver_parameters=self.solver_parameters,
+                    options_prefix=solver_name)
+                solver_list.append(solver)
+            return solver_list
 
     @cached_property
     def lhs(self):
         """Set up the discretisation's left hand side (the time derivative)."""
-        return super(ExplicitMultistage, self).lhs
+
+        if self.increment_form:
+            l = self.residual.label_map(
+                lambda t: t.has_label(time_derivative),
+                map_if_true=replace_subject(self.x_out, self.idx),
+                map_if_false=drop)
+
+            return l.form
+
+        else:
+            lhs_list = []
+            for stage in range(self.nStages):
+                l = self.residual.label_map(
+                    lambda t: t.has_label(time_derivative),
+                    map_if_true=replace_subject(self.field_i[stage+1], self.idx),
+                    map_if_false=drop)
+                lhs_list.append(l)
+
+            return lhs_list
 
     @cached_property
     def rhs(self):
         """Set up the time discretisation's right hand side."""
-        r = self.residual.label_map(
-            all_terms,
-            map_if_true=replace_subject(self.x1, old_idx=self.idx))
 
-        r = r.label_map(
-            lambda t: t.has_label(time_derivative),
-            map_if_true=drop,
-            map_if_false=lambda t: -1*t)
-
-        # If there are no active labels, we may have no terms at this point
-        # So that we can still do xnp1 = xn, put in a zero term here
-        if len(r.terms) == 0:
-            logger.warning('No terms detected for RHS of explicit problem. '
-                           + 'Adding a zero term to avoid failure.')
-            null_term = Constant(0.0)*self.residual.label_map(
+        if self.increment_form:
+            r = self.residual.label_map(
+                all_terms,
+                map_if_true=replace_subject(self.x1, old_idx=self.idx))
+
+            r = r.label_map(
                 lambda t: t.has_label(time_derivative),
-                # Drop label from this
-                map_if_true=lambda t: time_derivative.remove(t),
-                map_if_false=drop)
-            r += null_term
+                map_if_true=drop,
+                map_if_false=lambda t: -1*t)
+
+            # If there are no active labels, we may have no terms at this point
+            # So that we can still do xnp1 = xn, put in a zero term here
+            if self.increment_form and len(r.terms) == 0:
+                logger.warning('No terms detected for RHS of explicit problem. '
+                               + 'Adding a zero term to avoid failure.')
+                null_term = Constant(0.0)*self.residual.label_map(
+                    lambda t: t.has_label(time_derivative),
+                    # Drop label from this
+                    map_if_true=lambda t: time_derivative.remove(t),
+                    map_if_false=drop)
+                r += null_term
+
+            return r.form
 
-        return r.form
+        else:
+            rhs_list = []
+
+            for stage in range(self.nStages):
+                r = self.residual.label_map(
+                    all_terms,
+                    map_if_true=replace_subject(self.field_i[0], old_idx=self.idx))
+
+                r = r.label_map(
+                    lambda t: t.has_label(time_derivative),
+                    map_if_true=keep,
+                    map_if_false=lambda t: -self.butcher_matrix[stage, 0]*self.dt*t)
+
+                for i in range(1, stage+1):
+                    r_i = self.residual.label_map(
+                        lambda t: t.has_label(time_derivative),
+                        map_if_true=drop,
+                        map_if_false=replace_subject(self.field_i[i], old_idx=self.idx)
+                    )
+
+                    r -= self.butcher_matrix[stage, i]*self.dt*r_i
+
+                rhs_list.append(r)
+
+            return rhs_list
 
     def solve_stage(self, x0, stage):
-        self.x1.assign(x0)
-        for i in range(stage):
-            self.x1.assign(self.x1 + self.dt*self.butcher_matrix[stage-1, i]*self.k[i])
-        for evaluate in self.evaluate_source:
-            evaluate(self.x1, self.dt)
-        if self.limiter is not None:
-            self.limiter.apply(self.x1)
-        self.solver.solve()
-        self.k[stage].assign(self.x_out)
 
-        if (stage == self.nStages - 1):
+        if self.increment_form:
             self.x1.assign(x0)
-            for i in range(self.nStages):
-                self.x1.assign(self.x1 + self.dt*self.butcher_matrix[stage, i]*self.k[i])
-            self.x1.assign(self.x1)
 
+            # Use x0 as a first guess (otherwise may not converge)
+            self.x_out.assign(x0)
+
+            for i in range(stage):
+                self.x1.assign(self.x1 + self.dt*self.butcher_matrix[stage-1, i]*self.k[i])
+            for evaluate in self.evaluate_source:
+                evaluate(self.x1, self.dt)
             if self.limiter is not None:
                 self.limiter.apply(self.x1)
+            self.solver.solve()
+
+            self.k[stage].assign(self.x_out)
+
+            if (stage == self.nStages - 1):
+                self.x1.assign(x0)
+                for i in range(self.nStages):
+                    self.x1.assign(self.x1 + self.dt*self.butcher_matrix[stage, i]*self.k[i])
+                self.x1.assign(self.x1)
+
+                if self.limiter is not None:
+                    self.limiter.apply(self.x1)
+
+        else:
+            # Set initial field
+            if stage == 0:
+                self.field_i[0].assign(x0)
+
+            # Use x0 as a first guess (otherwise may not converge)
+            self.field_i[stage+1].assign(x0)
+
+            # Update field_i for physics / limiters
+            for evaluate in self.evaluate_source:
+                # TODO: not implemented! Here we need to evaluate the m-th term
+                # in the i-th RHS with field_m
+                raise NotImplementedError(
+                    'Physics not implemented with RK schemes that do not use '
+                    + 'the increment form')
+            if self.limiter is not None:
+                self.limiter.apply(self.field_i[stage])
+
+            # Obtain field_ip1 = field_n - dt* sum_m{a_im*F[field_m]}
+            self.solver[stage].solve()
+
+            if (stage == self.nStages - 1):
+                self.x1.assign(self.field_i[stage+1])
+                if self.limiter is not None:
+                    self.limiter.apply(self.x1)
 
     def apply_cycle(self, x_out, x_in):
         """
@@ -926,8 +1034,8 @@ class ForwardEuler(ExplicitMultistage):
     y^(n+1) = y^n + dt*k0                                                       \n
     """
     def __init__(self, domain, field_name=None, fixed_subcycles=None,
-                 subcycle_by_courant=None, solver_parameters=None,
-                 limiter=None, options=None):
+                 subcycle_by_courant=None, increment_form=True,
+                 solver_parameters=None, limiter=None, options=None):
         """
         Args:
             domain (:class:`Domain`): the model's domain object, containing the
@@ -943,6 +1051,9 @@ def __init__(self, domain, field_name=None, fixed_subcycles=None,
                 for one sub-cycle. Defaults to None, in which case adaptive
                 sub-cycling is not used. This option cannot be specified with the
                 `fixed_subcycles` argument.
+            increment_form (bool, optional): whether to write the RK scheme in
+                "increment form", solving for increments rather than updated
+                fields. Defaults to True.
             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
@@ -956,6 +1067,7 @@ def __init__(self, domain, field_name=None, fixed_subcycles=None,
         super().__init__(domain, butcher_matrix, field_name=field_name,
                          fixed_subcycles=fixed_subcycles,
                          subcycle_by_courant=subcycle_by_courant,
+                         increment_form=increment_form,
                          solver_parameters=solver_parameters,
                          limiter=limiter, options=options)
 
@@ -971,8 +1083,8 @@ class SSPRK3(ExplicitMultistage):
     y^(n+1) = y^n + (1/6)*dt*(k0 + k1 + 4*k2)                                 \n
     """
     def __init__(self, domain, field_name=None, fixed_subcycles=None,
-                 subcycle_by_courant=None, solver_parameters=None,
-                 limiter=None, options=None):
+                 subcycle_by_courant=None, increment_form=True,
+                 solver_parameters=None, limiter=None, options=None):
         """
         Args:
             domain (:class:`Domain`): the model's domain object, containing the
@@ -988,6 +1100,9 @@ def __init__(self, domain, field_name=None, fixed_subcycles=None,
                 for one sub-cycle. Defaults to None, in which case adaptive
                 sub-cycling is not used. This option cannot be specified with the
                 `fixed_subcycles` argument.
+            increment_form (bool, optional): whether to write the RK scheme in
+                "increment form", solving for increments rather than updated
+                fields. Defaults to True.
             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
@@ -998,9 +1113,11 @@ def __init__(self, domain, field_name=None, fixed_subcycles=None,
                 recovery method. Defaults to None.
         """
         butcher_matrix = np.array([[1., 0., 0.], [1./4., 1./4., 0.], [1./6., 1./6., 2./3.]])
+
         super().__init__(domain, butcher_matrix, field_name=field_name,
                          fixed_subcycles=fixed_subcycles,
                          subcycle_by_courant=subcycle_by_courant,
+                         increment_form=increment_form,
                          solver_parameters=solver_parameters,
                          limiter=limiter, options=options)
 
@@ -1021,7 +1138,8 @@ class RK4(ExplicitMultistage):
     where superscripts indicate the time-level.                               \n
     """
     def __init__(self, domain, field_name=None, fixed_subcycles=None,
-                 subcycle_by_courant=None, solver_parameters=None,
+                 subcycle_by_courant=None, increment_form=True,
+                 solver_parameters=None,
                  limiter=None, options=None):
         """
         Args:
@@ -1038,6 +1156,9 @@ def __init__(self, domain, field_name=None, fixed_subcycles=None,
                 for one sub-cycle. Defaults to None, in which case adaptive
                 sub-cycling is not used. This option cannot be specified with the
                 `fixed_subcycles` argument.
+            increment_form (bool, optional): whether to write the RK scheme in
+                "increment form", solving for increments rather than updated
+                fields. Defaults to True.
             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
@@ -1051,6 +1172,7 @@ def __init__(self, domain, field_name=None, fixed_subcycles=None,
         super().__init__(domain, butcher_matrix, field_name=field_name,
                          fixed_subcycles=fixed_subcycles,
                          subcycle_by_courant=subcycle_by_courant,
+                         increment_form=increment_form,
                          solver_parameters=solver_parameters,
                          limiter=limiter, options=options)
 
@@ -1069,8 +1191,8 @@ class Heun(ExplicitMultistage):
     number.
     """
     def __init__(self, domain, field_name=None, fixed_subcycles=None,
-                 subcycle_by_courant=None, solver_parameters=None,
-                 limiter=None, options=None):
+                 subcycle_by_courant=None, increment_form=True,
+                 solver_parameters=None, limiter=None, options=None):
         """
         Args:
             domain (:class:`Domain`): the model's domain object, containing the
@@ -1086,6 +1208,9 @@ def __init__(self, domain, field_name=None, fixed_subcycles=None,
                 for one sub-cycle. Defaults to None, in which case adaptive
                 sub-cycling is not used. This option cannot be specified with the
                 `fixed_subcycles` argument.
+            increment_form (bool, optional): whether to write the RK scheme in
+                "increment form", solving for increments rather than updated
+                fields. Defaults to True.
             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
@@ -1099,6 +1224,7 @@ def __init__(self, domain, field_name=None, fixed_subcycles=None,
         super().__init__(domain, butcher_matrix, field_name=field_name,
                          fixed_subcycles=fixed_subcycles,
                          subcycle_by_courant=subcycle_by_courant,
+                         increment_form=increment_form,
                          solver_parameters=solver_parameters,
                          limiter=limiter, options=options)
 
@@ -1901,7 +2027,6 @@ def apply(self, x_out, *x_in):
         """
         if self.initial_timesteps < self.nlevels-1:
             self.initial_timesteps += 1
-            print(self.initial_timesteps)
             solver = self.solver0
         else:
             solver = self.solver
diff --git a/gusto/transport_methods.py b/gusto/transport_methods.py
index 4fc951b6d..2329015ef 100644
--- a/gusto/transport_methods.py
+++ b/gusto/transport_methods.py
@@ -4,7 +4,7 @@
 
 from firedrake import (
     dx, dS, dS_v, dS_h, ds_t, ds_b, ds_v, dot, inner, outer, jump,
-    grad, div, FacetNormal, Function, sign, avg, cross, curl
+    grad, div, FacetNormal, Function, sign, avg, cross, curl, split
 )
 from firedrake.fml import Term, keep, drop
 from gusto.configuration import IntegrateByParts, TransportEquationType
@@ -36,6 +36,12 @@ def __init__(self, equation, variable):
 
         self.transport_equation_type = self.original_form.terms[0].get(transport)
 
+        if self.transport_equation_type == TransportEquationType.tracer_conservative:
+            # Extract associated density of the variable
+            tracer = next(x for x in self.equation.active_tracers if x.name == variable)
+            density_idx = self.equation.field_names.index(tracer.density_name)
+            self.conservative_density = split(self.equation.X)[density_idx]
+
     def replace_form(self, equation):
         """
         Replaces the form for the transport term in the equation with the
@@ -174,6 +180,8 @@ def __init__(self, equation, variable, ibp=IntegrateByParts.ONCE,
                 raise NotImplementedError('Outflow not implemented for upwind vector invariant')
             form = upwind_vector_invariant_form(self.domain, self.test, self.field, ibp=ibp)
 
+        elif self.transport_equation_type == TransportEquationType.tracer_conservative:
+            form = upwind_tracer_conservative_form(self.domain, self.test, self.field, self.conservative_density, ibp=ibp)
         else:
             raise NotImplementedError('Upwind transport scheme has not been '
                                       + 'implemented for this transport equation type')
@@ -300,6 +308,66 @@ def upwind_continuity_form(domain, test, q, ibp=IntegrateByParts.ONCE, outflow=F
     return ibp_label(transport(form, TransportEquationType.conservative), ibp)
 
 
+def upwind_tracer_conservative_form(domain, test, q, rho,
+                                    ibp=IntegrateByParts.ONCE, outflow=False):
+    u"""
+    The form corresponding to the DG upwind continuity transport operator.
+
+    This discretises ∇.(u*q*rho), for transporting velocity u, transported
+    variable q, and its reference density, rho. Although the tracer q obeys an advection
+    equation, the transport term is in a conservative form.
+
+    Args:
+        domain (:class:`Domain`): the model's domain object, containing the
+            mesh and the compatible function spaces.
+        test (:class:`TestFunction`): the test function.
+        q (:class:`ufl.Expr`): the variable to be transported.
+        rho (:class:`ufl.Expr`): the reference density for the tracer.
+        ibp (:class:`IntegrateByParts`, optional): an enumerator representing
+            the number of times to integrate by parts. Defaults to
+            `IntegrateByParts.ONCE`.
+        outflow (bool, optional): whether to include outflow at the domain
+            boundaries, through exterior facet terms. Defaults to False.
+
+    Raises:
+        ValueError: Can only use outflow option when the integration by parts
+            option is not "never".
+
+    Returns:
+        class:`LabelledForm`: a labelled transport form.
+    """
+
+    if outflow and ibp == IntegrateByParts.NEVER:
+        raise ValueError("outflow is True and ibp is None are incompatible options")
+    Vu = domain.spaces("HDiv")
+    dS_ = (dS_v + dS_h) if Vu.extruded else dS
+    ubar = Function(Vu)
+
+    if ibp == IntegrateByParts.ONCE:
+        L = -inner(grad(test), outer(inner(q, rho), ubar))*dx
+    else:
+        L = inner(test, div(outer(inner(q, rho), ubar)))*dx
+
+    if ibp != IntegrateByParts.NEVER:
+        n = FacetNormal(domain.mesh)
+        un = 0.5*(dot(ubar, n) + abs(dot(ubar, n)))
+
+        L += dot(jump(test), (un('+')*q('+')*rho('+') - un('-')*q('-')*rho('-')))*dS_
+
+        if ibp == IntegrateByParts.TWICE:
+            L -= (inner(test('+'), dot(ubar('+'), n('+'))*q('+')*rho('+'))
+                  + inner(test('-'), dot(ubar('-'), n('-'))*q('-')*rho('-')))*dS_
+
+    if outflow:
+        n = FacetNormal(domain.mesh)
+        un = 0.5*(dot(ubar, n) + abs(dot(ubar, n)))
+        L += test*un*q*rho*(ds_v + ds_t + ds_b)
+
+    form = transporting_velocity(L, ubar)
+
+    return ibp_label(transport(form, TransportEquationType.tracer_conservative), ibp)
+
+
 def vector_manifold_advection_form(domain, test, q, ibp=IntegrateByParts.ONCE, outflow=False):
     """
     Form for advective transport operator including vector manifold correction.
diff --git a/integration-tests/data/linear_sw_wave_chkpt.h5 b/integration-tests/data/linear_sw_wave_chkpt.h5
index 11420f7c7..03400f12c 100644
Binary files a/integration-tests/data/linear_sw_wave_chkpt.h5 and b/integration-tests/data/linear_sw_wave_chkpt.h5 differ
diff --git a/integration-tests/equations/test_coupled_transport.py b/integration-tests/equations/test_coupled_transport.py
index e22bbd81d..da474bb4f 100644
--- a/integration-tests/equations/test_coupled_transport.py
+++ b/integration-tests/equations/test_coupled_transport.py
@@ -6,7 +6,7 @@
 
 """
 
-from firedrake import norm
+from firedrake import norm, BrokenElement
 from gusto import *
 import pytest
 
@@ -63,3 +63,53 @@ def test_coupled_transport_scalar(tmpdir, geometry, equation_form1, equation_for
         'The transport error for tracer 1 is greater than the permitted tolerance'
     assert error2 < setup.tol, \
         'The transport error for tracer 2 is greater than the permitted tolerance'
+
+
+@pytest.mark.parametrize("m_X_space", ['DG', 'theta'])
+def test_conservative_coupled_transport(tmpdir, m_X_space, tracer_setup):
+
+    setup = tracer_setup(tmpdir, "slice")
+    domain = setup.domain
+    mesh = domain.mesh
+
+    rho_d = ActiveTracer(name='f1', space='DG',
+                         variable_type=TracerVariableType.density,
+                         transport_eqn=TransportEquationType.conservative)
+
+    m_X = ActiveTracer(name='f2', space=m_X_space,
+                       variable_type=TracerVariableType.mixing_ratio,
+                       transport_eqn=TransportEquationType.tracer_conservative,
+                       density_name='f1')
+
+    tracers = [rho_d, m_X]
+
+    V = domain.spaces("HDiv")
+    eqn = CoupledTransportEquation(domain, active_tracers=tracers, Vu=V)
+
+    if m_X_space == 'theta':
+        V_m_X = domain.spaces(m_X_space)
+        Vt_brok = FunctionSpace(mesh, BrokenElement(V_m_X.ufl_element()))
+        suboptions = {'f1': RecoveryOptions(embedding_space=Vt_brok,
+                                            recovered_space=V_m_X,
+                                            project_low_method='recover'),
+                      'f2': EmbeddedDGOptions()}
+        opts = MixedFSOptions(suboptions=suboptions)
+
+        transport_scheme = SSPRK3(domain, options=opts, increment_form=False)
+    else:
+        transport_scheme = SSPRK3(domain, increment_form=False)
+
+    transport_method = [DGUpwind(eqn, 'f1'), DGUpwind(eqn, 'f2')]
+
+    timestepper = PrescribedTransport(eqn, transport_scheme, setup.io, transport_method)
+
+    # Initial conditions
+    timestepper.fields("f1").interpolate(setup.f_init)
+    timestepper.fields("f2").interpolate(setup.f_init)
+    timestepper.fields("u").project(setup.uexpr)
+
+    error1, error2 = run(timestepper, setup.tmax, setup.f_end)
+    assert error1 < setup.tol, \
+        'The transport error for rho_d is greater than the permitted tolerance'
+    assert error2 < setup.tol, \
+        'The transport error for m_X is greater than the permitted tolerance'
diff --git a/integration-tests/model/test_time_discretisation.py b/integration-tests/model/test_time_discretisation.py
index 6f114fe90..9523ef75e 100644
--- a/integration-tests/model/test_time_discretisation.py
+++ b/integration-tests/model/test_time_discretisation.py
@@ -8,8 +8,10 @@ def run(timestepper, tmax, f_end):
     return norm(timestepper.fields("f") - f_end) / norm(f_end)
 
 
-@pytest.mark.parametrize("scheme", ["ssprk", "TrapeziumRule", "ImplicitMidpoint", "QinZhang",
-                                    "RK4", "Heun", "BDF2", "TR_BDF2", "AdamsBashforth", "Leapfrog", "AdamsMoulton"])
+@pytest.mark.parametrize(
+    "scheme", ["ssprk3_increment", "TrapeziumRule", "ImplicitMidpoint",
+               "QinZhang", "RK4", "Heun", "BDF2", "TR_BDF2", "AdamsBashforth",
+               "Leapfrog", "AdamsMoulton", "AdamsMoulton", "ssprk3_predictor"])
 def test_time_discretisation(tmpdir, scheme, tracer_setup):
     if (scheme == "AdamsBashforth"):
         # Tighter stability constraints
@@ -25,8 +27,10 @@ def test_time_discretisation(tmpdir, scheme, tracer_setup):
         V = domain.spaces("DG")
         eqn = AdvectionEquation(domain, V, "f")
 
-    if scheme == "ssprk":
-        transport_scheme = SSPRK3(domain)
+    if scheme == "ssprk3_increment":
+        transport_scheme = SSPRK3(domain, increment_form=True)
+    elif scheme == "ssprk3_predictor":
+        transport_scheme = SSPRK3(domain, increment_form=False)
     elif scheme == "TrapeziumRule":
         transport_scheme = TrapeziumRule(domain)
     elif scheme == "ImplicitMidpoint":
diff --git a/integration-tests/physics/test_suppress_vertical_wind.py b/integration-tests/physics/test_suppress_vertical_wind.py
index 44f1e05ae..cab6efc76 100644
--- a/integration-tests/physics/test_suppress_vertical_wind.py
+++ b/integration-tests/physics/test_suppress_vertical_wind.py
@@ -92,5 +92,5 @@ def test_suppress_vertical_wind(tmpdir):
     vertical_wind = Function(domain.spaces('theta'))
     vertical_wind.interpolate(dot(u, domain.k))
 
-    tol = 1e-12
+    tol = 1e-10
     assert norm(vertical_wind) < tol
diff --git a/unit-tests/test_active_tracer.py b/unit-tests/test_active_tracer.py
index 8c7907c43..65c8ff56b 100644
--- a/unit-tests/test_active_tracer.py
+++ b/unit-tests/test_active_tracer.py
@@ -11,7 +11,8 @@ def test_tracer_classes():
     names = ['mr_v', 'big_blob']
     spaces = ['V', 'U']
     transport_eqns = [TransportEquationType.advective,
-                      TransportEquationType.no_transport]
+                      TransportEquationType.no_transport,
+                      TransportEquationType.tracer_conservative]
     variable_types = [TracerVariableType.mixing_ratio]
 
     for name, space, transport_eqn in zip(names, spaces, transport_eqns):