diff --git a/gusto/core/__init__.py b/gusto/core/__init__.py
index 163f9ddf..432c66b9 100644
--- a/gusto/core/__init__.py
+++ b/gusto/core/__init__.py
@@ -1,11 +1,12 @@
-from gusto.core.configuration import *           # noqa
-from gusto.core.coordinates import *             # noqa
-from gusto.core.coord_transforms import *        # noqa
-from gusto.core.domain import *                  # noqa
-from gusto.core.fields import *                  # noqa
-from gusto.core.function_spaces import *         # noqa
-from gusto.core.io import *                      # noqa
-from gusto.core.kernels import *                 # noqa
-from gusto.core.labels import *                  # noqa
-from gusto.core.logging import *                 # noqa
-from gusto.core.meshes import *                  # noqa
\ No newline at end of file
+from gusto.core.configuration import *             # noqa
+from gusto.core.conservative_projection import *   # noqa
+from gusto.core.coordinates import *               # noqa
+from gusto.core.coord_transforms import *          # noqa
+from gusto.core.domain import *                    # noqa
+from gusto.core.fields import *                    # noqa
+from gusto.core.function_spaces import *           # noqa
+from gusto.core.io import *                        # noqa
+from gusto.core.kernels import *                   # noqa
+from gusto.core.labels import *                    # noqa
+from gusto.core.logging import *                   # noqa
+from gusto.core.meshes import *                    # noqa
\ No newline at end of file
diff --git a/gusto/core/configuration.py b/gusto/core/configuration.py
index 252f2518..a3796b40 100644
--- a/gusto/core/configuration.py
+++ b/gusto/core/configuration.py
@@ -8,7 +8,8 @@
     "IntegrateByParts", "TransportEquationType", "OutputParameters",
     "BoussinesqParameters", "CompressibleParameters",
     "ShallowWaterParameters",
-    "EmbeddedDGOptions", "RecoveryOptions", "SUPGOptions", "MixedFSOptions",
+    "EmbeddedDGOptions", "ConservativeEmbeddedDGOptions", "RecoveryOptions",
+    "ConservativeRecoveryOptions", "SUPGOptions", "MixedFSOptions",
     "SpongeLayerParameters", "DiffusionParameters", "BoundaryLayerParameters"
 ]
 
@@ -164,6 +165,14 @@ class EmbeddedDGOptions(WrapperOptions):
     embedding_space = None
 
 
+class ConservativeEmbeddedDGOptions(EmbeddedDGOptions):
+    """Specifies options for a conservative embedded DG method."""
+
+    project_back_method = 'conservative_project'
+    rho_name = None
+    orig_rho_space = None
+
+
 class RecoveryOptions(WrapperOptions):
     """Specifies options for a recovery wrapper method."""
 
@@ -177,6 +186,15 @@ class RecoveryOptions(WrapperOptions):
     broken_method = 'interpolate'
 
 
+class ConservativeRecoveryOptions(RecoveryOptions):
+    """Specifies options for a conservative recovery wrapper method."""
+
+    rho_name = None
+    orig_rho_space = None
+    project_high_method = 'conservative_project'
+    project_low_method = 'conservative_project'
+
+
 class SUPGOptions(WrapperOptions):
     """Specifies options for an SUPG scheme."""
 
diff --git a/gusto/core/conservative_projection.py b/gusto/core/conservative_projection.py
new file mode 100644
index 00000000..1ab1455d
--- /dev/null
+++ b/gusto/core/conservative_projection.py
@@ -0,0 +1,93 @@
+"""
+This provides an operator for perform a conservative projection.
+
+The :class:`ConservativeProjector` provided in this module is an operator that
+projects a field such as a mixing ratio from one function space to another,
+weighted by a density field to ensure that mass is conserved by the projection.
+"""
+
+from firedrake import (Function, TestFunction, TrialFunction, lhs, rhs, inner,
+                       dx, LinearVariationalProblem, LinearVariationalSolver,
+                       Constant, assemble)
+import ufl
+
+__all__ = ["ConservativeProjector"]
+
+
+class ConservativeProjector(object):
+    """
+    Projects a field such that mass is conserved.
+
+    This object is designed for projecting fields such as mixing ratios of
+    tracer species from one function space to another, but weighted by density
+    such that mass is conserved by the projection.
+    """
+
+    def __init__(self, rho_source, rho_target, m_source, m_target,
+                 subtract_mean=False):
+        """
+        Args:
+            rho_source (:class:`Function`): the density to use for weighting the
+                source mixing ratio field. Can also be a :class:`ufl.Expr`.
+            rho_target (:class:`Function`): the density to use for weighting the
+                target mixing ratio field. Can also be a :class:`ufl.Expr`.
+            m_source (:class:`Function`): the source mixing ratio field. Can
+                also be a :class:`ufl.Expr`.
+            m_target (:class:`Function`): the target mixing ratio field to
+                compute.
+            subtract_mean (bool, optional): whether to solve the projection by
+                subtracting the mean value of m for both sides. This is more
+                expensive as it involves calculating the mean, but will ensure
+                preservation of a constant when projecting to a continuous
+                space. Default to False.
+
+        Raises:
+            RuntimeError: the geometric shape of the two rho fields must be equal.
+            RuntimeError: the geometric shape of the two m fields must be equal.
+        """
+
+        self.subtract_mean = subtract_mean
+
+        if not isinstance(rho_source, (ufl.core.expr.Expr, Function)):
+            raise ValueError("Can only recover UFL expression or Functions not '%s'" % type(rho_source))
+
+        if not isinstance(rho_target, (ufl.core.expr.Expr, Function)):
+            raise ValueError("Can only recover UFL expression or Functions not '%s'" % type(rho_target))
+
+        if not isinstance(m_source, (ufl.core.expr.Expr, Function)):
+            raise ValueError("Can only recover UFL expression or Functions not '%s'" % type(m_source))
+
+        # Check shape values
+        if m_source.ufl_shape != m_target.ufl_shape:
+            raise RuntimeError('Shape mismatch between source %s and target function spaces %s in project' % (m_source.ufl_shape, m_target.ufl_shape))
+
+        if rho_source.ufl_shape != rho_target.ufl_shape:
+            raise RuntimeError('Shape mismatch between source %s and target function spaces %s in project' % (rho_source.ufl_shape, rho_target.ufl_shape))
+
+        self.m_source = m_source
+        self.m_target = m_target
+
+        V = self.m_target.function_space()
+        mesh = V.mesh()
+
+        self.m_mean = Constant(0.0, domain=mesh)
+        self.volume = assemble(Constant(1.0, domain=mesh)*dx)
+
+        test = TestFunction(V)
+        m_trial = TrialFunction(V)
+        eqn = (rho_source*inner(test, m_source - self.m_mean)*dx
+               - rho_target*inner(test, m_trial - self.m_mean)*dx)
+        problem = LinearVariationalProblem(lhs(eqn), rhs(eqn), self.m_target)
+        self.solver = LinearVariationalSolver(problem)
+
+    def project(self):
+        """Apply the projection."""
+
+        # Compute mean value
+        if self.subtract_mean:
+            self.m_mean.assign(assemble(self.m_source*dx) / self.volume)
+
+        # Solve projection
+        self.solver.solve()
+
+        return self.m_target
diff --git a/gusto/equations/prognostic_equations.py b/gusto/equations/prognostic_equations.py
index bf8dd784..985e0ff5 100644
--- a/gusto/equations/prognostic_equations.py
+++ b/gusto/equations/prognostic_equations.py
@@ -304,6 +304,17 @@ def add_tracers_to_prognostics(self, domain, active_tracers):
                 name of the active tracer.
         """
 
+        # Check if there are any conservatively transported tracers.
+        # If so, ensure that the reference density is indexed before this tracer.
+        for i in range(len(active_tracers) - 1):
+            tracer = active_tracers[i]
+            if tracer.transport_eqn == TransportEquationType.tracer_conservative:
+                ref_density = next(x for x in active_tracers if x.name == tracer.density_name)
+                j = active_tracers.index(ref_density)
+                if j > i:
+                    # Swap the indices of the tracer and the reference density
+                    active_tracers[i], active_tracers[j] = active_tracers[j], active_tracers[i]
+
         # Loop through tracer fields and add field names and spaces
         for tracer in active_tracers:
             if isinstance(tracer, ActiveTracer):
diff --git a/gusto/recovery/reversible_recovery.py b/gusto/recovery/reversible_recovery.py
index d9ad661d..fc8fbd33 100644
--- a/gusto/recovery/reversible_recovery.py
+++ b/gusto/recovery/reversible_recovery.py
@@ -3,9 +3,12 @@
 higher-order function space.
 """
 
+from gusto.core.conservative_projection import ConservativeProjector
 from firedrake import (Projector, Function, Interpolator)
 from .recovery import Recoverer
 
+__all__ = ["ReversibleRecoverer", "ConservativeRecoverer"]
+
 
 class ReversibleRecoverer(object):
     """
@@ -13,10 +16,11 @@ class ReversibleRecoverer(object):
     field into a higher-order discontinuous space. This uses the recovery
     operator, but with further adjustments to ensure reversibility.
 
-    :arg source_field:      the source field.
-    :arg target_field:      the target_field.
-    :arg reconstruct_opts:  an object containing the various options for the
-                               reconstruction.
+    Args:
+        source_field (:class:`Function`): the source field.
+        target_field (:class:`Function`): the target field.
+        reconstruct_opts (:class:`RecoveryOptions`): an object containing the
+            various options for the reconstruction.
     """
     def __init__(self, source_field, target_field, reconstruct_opts):
 
@@ -92,3 +96,63 @@ def project(self):
         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_high.assign(self.q_corr_high + self.q_rec_high)
+
+
+class ConservativeRecoverer(object):
+    """
+    An object for performing a reconstruction of a low-order discontinuous
+    field into a higher-order discontinuous space, but such that mass is
+    conserved. This uses the recovery operator, but with further adjustments to
+    ensure both reversibility and mass conservation.
+
+    Args:
+        source_field (:class:`Function`): the source field.
+        target_field (:class:`Function`): the target field.
+        source_density (:class:`Function`): the source density field.
+        target_density (:class:`Function`): the target density field.
+        reconstruct_opts (:class:`RecoveryOptions`): an object containing the
+            various options for the reconstruction.
+    """
+    def __init__(self, source_field, target_field, source_density,
+                 target_density, reconstruct_opts):
+
+        self.opts = reconstruct_opts
+
+        # Declare the fields used by the reconstructor
+        self.q_low = source_field
+        self.q_high = target_field
+        self.q_recovered = Function(self.opts.recovered_space)
+        self.q_corr_low = Function(source_field.function_space())
+        self.q_corr_high = Function(target_field.function_space())
+        self.q_rec_high = Function(target_field.function_space())
+
+        # -------------------------------------------------------------------- #
+        # Set up the operators for different transformations
+        # -------------------------------------------------------------------- #
+
+        # Does recovery by first projecting into broken space then averaging
+        self.recoverer = Recoverer(self.q_low, self.q_recovered,
+                                   method=self.opts.broken_method,
+                                   boundary_method=self.opts.boundary_method)
+
+        # Obtain the recovered field in the higher order space
+        self.projector_high = Projector(self.q_recovered, self.q_rec_high)
+
+        # Obtain the correction in the lower order space
+        # Swap density arguments!
+        self.projector_low = ConservativeProjector(target_density, source_density,
+                                                   self.q_rec_high, self.q_corr_low,
+                                                   subtract_mean=True)
+
+        # Final injection operator
+        # Should identify low order field in higher order space
+        self.injector = ConservativeProjector(source_density, target_density,
+                                              self.q_corr_low, self.q_corr_high)
+
+    def project(self):
+        self.recoverer.project()
+        self.projector_high.project()
+        self.projector_low.project()
+        self.q_corr_low.assign(self.q_low - self.q_corr_low)
+        self.injector.project()
+        self.q_high.assign(self.q_corr_high + self.q_rec_high)
diff --git a/gusto/time_discretisation/time_discretisation.py b/gusto/time_discretisation/time_discretisation.py
index eddb7245..1aef8ad6 100644
--- a/gusto/time_discretisation/time_discretisation.py
+++ b/gusto/time_discretisation/time_discretisation.py
@@ -94,7 +94,7 @@ def __init__(self, domain, field_name=None, solver_parameters=None,
                             'Time discretisation: suboption SUPG is currently not implemented within MixedOptions')
                     else:
                         raise RuntimeError(
-                            f'Time discretisation: suboption wrapper {wrapper_name} not implemented')
+                            f'Time discretisation: suboption wrapper {suboption.name} not implemented')
             elif self.wrapper_name == "embedded_dg":
                 self.wrapper = EmbeddedDGWrapper(self, options)
             elif self.wrapper_name == "recovered":
diff --git a/gusto/time_discretisation/wrappers.py b/gusto/time_discretisation/wrappers.py
index f37e402c..a1803cad 100644
--- a/gusto/time_discretisation/wrappers.py
+++ b/gusto/time_discretisation/wrappers.py
@@ -11,8 +11,9 @@
 )
 from firedrake.fml import Term
 from gusto.core.configuration import EmbeddedDGOptions, RecoveryOptions, SUPGOptions
-from gusto.recovery import Recoverer, ReversibleRecoverer
+from gusto.recovery import Recoverer, ReversibleRecoverer, ConservativeRecoverer
 from gusto.core.labels import transporting_velocity
+from gusto.core.conservative_projection import ConservativeProjector
 import ufl
 
 __all__ = ["EmbeddedDGWrapper", "RecoveryWrapper", "SUPGWrapper", "MixedFSWrapper"]
@@ -34,6 +35,7 @@ def __init__(self, time_discretisation, wrapper_options):
         self.options = wrapper_options
         self.solver_parameters = None
         self.original_space = None
+        self.is_conservative = False
 
     @abstractmethod
     def setup(self, original_space):
@@ -123,6 +125,7 @@ def setup(self, original_space, post_apply_bcs):
 
         self.x_in = Function(self.function_space)
         self.x_out = Function(self.function_space)
+        self.x_in_orig = Function(original_space)
 
         if self.time_discretisation.idx is None:
             self.x_projected = Function(self.original_space)
@@ -134,6 +137,19 @@ def setup(self, original_space, post_apply_bcs):
                                              bcs=post_apply_bcs)
         elif self.options.project_back_method == 'recover':
             self.x_out_projector = Recoverer(self.x_out, self.x_projected)
+        elif self.options.project_back_method == 'conservative_project':
+            self.is_conservative = True
+            self.rho_name = self.options.rho_name
+            self.rho_in_orig = Function(self.options.orig_rho_space)
+            self.rho_out_orig = Function(self.options.orig_rho_space)
+            self.rho_in_embedded = Function(self.function_space)
+            self.rho_out_embedded = Function(self.function_space)
+            self.x_in_projector = ConservativeProjector(
+                self.rho_in_orig, self.rho_in_embedded,
+                self.x_in_orig, self.x_in)
+            self.x_out_projector = ConservativeProjector(
+                self.rho_out_embedded, self.rho_out_orig,
+                self.x_out, self.x_projected, subtract_mean=True)
         else:
             raise NotImplementedError(
                 'EmbeddedDG Wrapper: project_back_method'
@@ -152,10 +168,15 @@ def pre_apply(self, x_in):
             x_in (:class:`Function`): the original input field.
         """
 
-        try:
-            self.x_in.interpolate(x_in)
-        except NotImplementedError:
-            self.x_in.project(x_in)
+        self.x_in_orig.assign(x_in)
+
+        if self.is_conservative:
+            self.x_in_projector.project()
+        else:
+            try:
+                self.x_in.interpolate(x_in)
+            except NotImplementedError:
+                self.x_in.project(x_in)
 
     def post_apply(self, x_out):
         """
@@ -215,7 +236,7 @@ def setup(self, original_space, post_apply_bcs):
         # Internal variables to be used
         # -------------------------------------------------------------------- #
 
-        self.x_in_tmp = Function(self.original_space)
+        self.x_in_orig = Function(self.original_space)
         self.x_in = Function(self.function_space)
         self.x_out = Function(self.function_space)
 
@@ -225,7 +246,19 @@ def setup(self, original_space, post_apply_bcs):
             self.x_projected = Function(equation.spaces[self.time_discretisation.idx])
 
         # Operator to recover to higher discontinuous space
-        self.x_recoverer = ReversibleRecoverer(self.x_in_tmp, self.x_in, self.options)
+        if self.options.project_low_method == 'conservative_project':
+            self.is_conservative = True
+            self.rho_name = self.options.rho_name
+            self.rho_in_orig = Function(self.options.orig_rho_space)
+            self.rho_out_orig = Function(self.options.orig_rho_space)
+            self.rho_in_embedded = Function(self.function_space)
+            self.rho_out_embedded = Function(self.function_space)
+            self.x_recoverer = ConservativeRecoverer(self.x_in_orig, self.x_in,
+                                                     self.rho_in_orig,
+                                                     self.rho_in_embedded,
+                                                     self.options)
+        else:
+            self.x_recoverer = ReversibleRecoverer(self.x_in_orig, self.x_in, self.options)
 
         # Operators for projecting back
         self.interp_back = (self.options.project_low_method == 'interpolate')
@@ -237,6 +270,10 @@ def setup(self, original_space, post_apply_bcs):
         elif self.options.project_low_method == 'recover':
             self.x_out_projector = Recoverer(self.x_out, self.x_projected,
                                              method=self.options.broken_method)
+        elif self.options.project_low_method == 'conservative_project':
+            self.x_out_projector = ConservativeProjector(
+                self.rho_out_embedded, self.rho_out_orig,
+                self.x_out, self.x_projected, subtract_mean=True)
         else:
             raise NotImplementedError(
                 'Recovery Wrapper: project_back_method'
@@ -251,7 +288,7 @@ def pre_apply(self, x_in):
             x_in (:class:`Function`): the original input field.
         """
 
-        self.x_in_tmp.assign(x_in)
+        self.x_in_orig.assign(x_in)
         self.x_recoverer.project()
 
     def post_apply(self, x_out):
@@ -416,6 +453,7 @@ def setup(self):
         self.function_space = MixedFunctionSpace(self.wrapper_spaces)
         self.x_in = Function(self.function_space)
         self.x_out = Function(self.function_space)
+        self.is_conservative = any([subwrapper.is_conservative for subwrapper in self.subwrappers.values()])
 
     def pre_apply(self, x_in):
         """
@@ -430,6 +468,8 @@ def pre_apply(self, x_in):
 
             if field_name in self.subwrappers:
                 subwrapper = self.subwrappers[field_name]
+                if subwrapper.is_conservative:
+                    self.pre_update_rho(subwrapper)
                 subwrapper.pre_apply(field)
                 x_in_sub.assign(subwrapper.x_in)
             else:
@@ -449,6 +489,34 @@ def post_apply(self, x_out):
             if field_name in self.subwrappers:
                 subwrapper = self.subwrappers[field_name]
                 subwrapper.x_out.assign(field)
+                if subwrapper.is_conservative:
+                    self.post_update_rho(subwrapper)
                 subwrapper.post_apply(x_out_sub)
             else:
                 x_out_sub.assign(field)
+
+    def pre_update_rho(self, subwrapper):
+        """
+        Updates the stored density field for the pre-apply for the subwrapper.
+
+        Args:
+            subwrapper (:class:`Wrapper`): the original input field.
+        """
+
+        rho_subwrapper = self.subwrappers[subwrapper.rho_name]
+
+        subwrapper.rho_in_orig.assign(rho_subwrapper.x_in_orig)
+        subwrapper.rho_in_embedded.assign(rho_subwrapper.x_in)
+
+    def post_update_rho(self, subwrapper):
+        """
+        Updates the stored density field for the post-apply for the subwrapper.
+
+        Args:
+            subwrapper (:class:`Wrapper`): the original input field.
+        """
+
+        rho_subwrapper = self.subwrappers[subwrapper.rho_name]
+
+        subwrapper.rho_out_orig.assign(rho_subwrapper.x_projected)
+        subwrapper.rho_out_embedded.assign(rho_subwrapper.x_out)
diff --git a/integration-tests/transport/test_tracer_conservative_transport.py b/integration-tests/transport/test_tracer_conservative_transport.py
new file mode 100644
index 00000000..8ff8b549
--- /dev/null
+++ b/integration-tests/transport/test_tracer_conservative_transport.py
@@ -0,0 +1,206 @@
+"""
+Tests the conservative transport of a mixing ratio and dry density, both when
+they are defined on the same and different function spaces. This checks
+that there is conservation of the total species mass (dry density times the
+mixing ratio) and that there is consistency (a constant field will remain
+constant).
+"""
+
+from gusto import *
+from firedrake import (
+    PeriodicIntervalMesh, ExtrudedMesh, exp, cos, sin, SpatialCoordinate,
+    assemble, dx, FunctionSpace, pi, min_value, as_vector, BrokenElement,
+    errornorm
+)
+import pytest
+
+
+def setup_conservative_transport(dirname, pair_of_spaces, desirable_property):
+
+    # Domain
+    Lx = 2000.
+    Hz = 2000.
+
+    # Time parameters
+    dt = 2.
+    tmax = 2000.
+
+    nlayers = 10.  # horizontal layers
+    columns = 10.  # number of columns
+
+    # Define the spaces for the tracers
+    if pair_of_spaces == 'same_order_1':
+        rho_d_space = 'DG'
+        m_X_space = 'DG'
+        space_order = 1
+    elif pair_of_spaces == 'diff_order_0':
+        rho_d_space = 'DG'
+        m_X_space = 'theta'
+        space_order = 0
+    elif pair_of_spaces == 'diff_order_1':
+        rho_d_space = 'DG'
+        m_X_space = 'theta'
+        space_order = 1
+
+    period_mesh = PeriodicIntervalMesh(columns, Lx)
+    mesh = ExtrudedMesh(period_mesh, layers=nlayers, layer_height=Hz/nlayers)
+    domain = Domain(mesh, dt, "CG", space_order)
+    x, z = SpatialCoordinate(mesh)
+
+    V_rho = domain.spaces(rho_d_space)
+    V_m_X = domain.spaces(m_X_space)
+
+    m_X = ActiveTracer(name='m_X', space=m_X_space,
+                       variable_type=TracerVariableType.mixing_ratio,
+                       transport_eqn=TransportEquationType.tracer_conservative,
+                       density_name='rho_d')
+
+    rho_d = ActiveTracer(name='rho_d', space=rho_d_space,
+                         variable_type=TracerVariableType.density,
+                         transport_eqn=TransportEquationType.conservative)
+
+    # Define m_X first to test that the tracers will be
+    # automatically re-ordered such that the density field
+    # is indexed before the mixing ratio.
+    tracers = [m_X, rho_d]
+
+    # Equation
+    V = domain.spaces("HDiv")
+    eqn = CoupledTransportEquation(domain, active_tracers=tracers, Vu=V)
+
+    # IO
+    output = OutputParameters(dirname=dirname)
+    io = IO(domain, output)
+
+    if pair_of_spaces == 'diff_order_0':
+        VCG1 = FunctionSpace(mesh, 'CG', 1)
+        VDG1 = domain.spaces('DG1_equispaced')
+
+        suboptions = {
+            'rho_d': RecoveryOptions(
+                embedding_space=VDG1,
+                recovered_space=VCG1,
+                project_low_method='recover',
+                boundary_method=BoundaryMethod.taylor
+            ),
+            'm_X': ConservativeRecoveryOptions(
+                embedding_space=VDG1,
+                recovered_space=VCG1,
+                boundary_method=BoundaryMethod.taylor,
+                rho_name='rho_d',
+                orig_rho_space=V_rho
+            )
+        }
+    elif pair_of_spaces == 'diff_order_1':
+        Vt_brok = FunctionSpace(mesh, BrokenElement(V_m_X.ufl_element()))
+        suboptions = {
+            'rho_d': EmbeddedDGOptions(embedding_space=Vt_brok),
+            'm_X': ConservativeEmbeddedDGOptions(
+                rho_name='rho_d',
+                orig_rho_space=V_rho
+            )
+        }
+    else:
+        suboptions = {}
+
+    opts = MixedFSOptions(suboptions=suboptions)
+
+    transport_scheme = SSPRK3(domain, options=opts, increment_form=False)
+    transport_methods = [DGUpwind(eqn, "m_X"), DGUpwind(eqn, "rho_d")]
+
+    # Timestepper
+    time_varying = True
+    stepper = PrescribedTransport(
+        eqn, transport_scheme, io, time_varying, transport_methods
+    )
+
+    # Initial Conditions
+    # Specify locations of the two Gaussians
+    xc1 = 5.*Lx/8.
+    zc1 = Hz/2.
+
+    xc2 = 3.*Lx/8.
+    zc2 = Hz/2.
+
+    def l2_dist(xc, zc):
+        return min_value(abs(x-xc), Lx-abs(x-xc))**2 + (z-zc)**2
+
+    lc = 2.*Lx/25.
+    m0 = 0.02
+
+    # Set the initial state from the configuration choice
+    if desirable_property == 'conservation':
+        f0 = 0.05
+
+        rho_t = 0.5
+        rho_b = 1.
+
+        rho_d_0 = rho_b + z*(rho_t-rho_b)/Hz
+
+        g1 = f0*exp(-l2_dist(xc1, zc1)/(lc**2))
+        g2 = f0*exp(-l2_dist(xc2, zc2)/(lc**2))
+
+        m_X_0 = m0 + g1 + g2
+
+    else:
+        f0 = 0.5
+        rho_b = 0.5
+
+        g1 = f0*exp(-l2_dist(xc1, zc1)/(lc**2))
+        g2 = f0*exp(-l2_dist(xc2, zc2)/(lc**2))
+
+        rho_d_0 = rho_b + g1 + g2
+
+        # Constant mass field
+        m_X_0 = m0 + 0*x
+
+    # Set up the divergent, time-varying, velocity field
+    U = Lx/tmax
+    W = U/10.
+
+    def u_t(t):
+        xd = x - U*t
+        u = U - (W*pi*Lx/Hz)*cos(pi*t/tmax)*cos(2*pi*xd/Lx)*cos(pi*z/Hz)
+        w = 2*pi*W*cos(pi*t/tmax)*sin(2*pi*xd/Lx)*sin(pi*z/Hz)
+
+        u_expr = as_vector((u, w))
+
+        return u_expr
+
+    stepper.setup_prescribed_expr(u_t)
+
+    stepper.fields("m_X").interpolate(m_X_0)
+    stepper.fields("rho_d").interpolate(rho_d_0)
+    stepper.fields("u").project(u_t(0))
+
+    m_X_init = Function(V_m_X)
+    rho_d_init = Function(V_rho)
+
+    m_X_init.assign(stepper.fields("m_X"))
+    rho_d_init.assign(stepper.fields("rho_d"))
+
+    return stepper, m_X_init, rho_d_init
+
+
+@pytest.mark.parametrize("pair_of_spaces", ["same_order_1", "diff_order_0", "diff_order_1"])
+@pytest.mark.parametrize("desirable_property", ["consistency", "conservation"])
+def test_conservative_transport(tmpdir, pair_of_spaces, desirable_property):
+
+    # Setup and run
+    dirname = str(tmpdir)
+
+    stepper, m_X_0, rho_d_0 = \
+        setup_conservative_transport(dirname, pair_of_spaces, desirable_property)
+
+    # Run for five timesteps
+    stepper.run(t=0, tmax=10)
+    m_X = stepper.fields("m_X")
+    rho_d = stepper.fields("rho_d")
+
+    # Perform the check
+    if desirable_property == 'consistency':
+        assert errornorm(m_X_0, m_X) < 2e-13, "conservative transport is not consistent"
+    else:
+        rho_X_init = assemble(m_X_0*rho_d_0*dx)
+        rho_X_final = assemble(m_X*rho_d*dx)
+        assert abs((rho_X_init - rho_X_final)/rho_X_init) < 1e-14, "conservative transport is not conservative"
diff --git a/unit-tests/recovery_tests/test_conservative_recovery.py b/unit-tests/recovery_tests/test_conservative_recovery.py
new file mode 100644
index 00000000..796e7539
--- /dev/null
+++ b/unit-tests/recovery_tests/test_conservative_recovery.py
@@ -0,0 +1,112 @@
+"""
+Test whether the conservative recovery process is working appropriately.
+"""
+
+from firedrake import (PeriodicIntervalMesh, IntervalMesh, ExtrudedMesh,
+                       SpatialCoordinate, FiniteElement, FunctionSpace,
+                       TensorProductElement, Function, interval, norm, errornorm,
+                       assemble)
+from gusto import *
+import numpy as np
+import pytest
+
+np.random.seed(0)
+
+
+@pytest.fixture
+def mesh(geometry):
+
+    L = 100.
+    H = 100.
+
+    deltax = L / 5.
+    deltaz = H / 5.
+    nlayers = int(H/deltaz)
+    ncolumns = int(L/deltax)
+
+    if geometry == "periodic":
+        m = PeriodicIntervalMesh(ncolumns, L)
+    elif geometry == "non-periodic":
+        m = IntervalMesh(ncolumns, L)
+
+    extruded_mesh = ExtrudedMesh(m, layers=nlayers, layer_height=deltaz)
+
+    return extruded_mesh
+
+
+def expr(geometry, mesh, configuration):
+
+    x, z = SpatialCoordinate(mesh)
+
+    if configuration == 'rho_constant':
+        rho_expr = Constant(2.0)
+        if geometry == "periodic":
+            m_expr = np.random.randn() + np.random.randn() * z
+        elif geometry == "non-periodic":
+            m_expr = np.random.randn() + np.random.randn() * x + np.random.randn() * z
+
+    elif configuration == 'm_constant':
+        m_expr = Constant(0.01)
+        if geometry == "periodic":
+            rho_expr = np.random.randn() + np.random.randn() * z
+        elif geometry == "non-periodic":
+            rho_expr = np.random.randn() + np.random.randn() * x + np.random.randn() * z
+
+    return rho_expr, m_expr
+
+
+@pytest.mark.parametrize("configuration", ["m_constant", "rho_constant"])
+@pytest.mark.parametrize("geometry", ["periodic", "non-periodic"])
+def test_conservative_recovery(geometry, mesh, configuration):
+
+    rho_expr, m_expr = expr(geometry, mesh, configuration)
+
+    # construct theta elemnt
+    cell = mesh._base_mesh.ufl_cell().cellname()
+    w_hori = FiniteElement("DG", cell, 0)
+    w_vert = FiniteElement("CG", interval, 1)
+    theta_element = TensorProductElement(w_hori, w_vert)
+
+    # spaces
+    DG0 = FunctionSpace(mesh, "DG", 0)
+    CG1 = FunctionSpace(mesh, "CG", 1)
+    DG1 = FunctionSpace(mesh, "DG", 1)
+    Vt = FunctionSpace(mesh, theta_element)
+
+    # set up density
+    rho_DG1 = Function(DG1).interpolate(rho_expr)
+    rho_DG0 = Function(DG0).project(rho_DG1)
+
+    # mixing ratio fields
+    m_Vt = Function(Vt).interpolate(m_expr)
+    m_DG1_approx = Function(DG1).interpolate(m_expr)
+    m_Vt_back = Function(Vt)
+    m_DG1 = Function(DG1)
+
+    options = ConservativeRecoveryOptions(embedding_space=DG1,
+                                          recovered_space=CG1,
+                                          boundary_method=BoundaryMethod.taylor)
+
+    # make the recoverers and do the recovery
+    conservative_recoverer = ConservativeRecoverer(m_Vt, m_DG1,
+                                                   rho_DG0, rho_DG1, options)
+    back_projector = ConservativeProjector(rho_DG1, rho_DG0, m_DG1, m_Vt_back,
+                                           subtract_mean=True)
+
+    conservative_recoverer.project()
+    back_projector.project()
+
+    # check various aspects of the process
+    m_high_diff = errornorm(m_DG1, m_DG1_approx) / norm(m_DG1_approx)
+    m_low_diff = errornorm(m_Vt_back, m_Vt) / norm(m_Vt)
+    mass_low = assemble(rho_DG0*m_Vt*dx)
+    mass_high = assemble(rho_DG1*m_DG1*dx)
+
+    assert (mass_low - mass_high) / mass_high < 5e-14, \
+        f'Conservative recovery on {geometry} vertical slice not conservative for {configuration} configuration'
+    assert m_low_diff < 2e-14, \
+        f'Irreversible conservative recovery on {geometry} vertical slice for {configuration} configuration'
+
+    if configuration in ['m_constant', 'rho_constant']:
+        assert m_high_diff < 2e-14, \
+            f'Inaccurate conservative recovery on {geometry} vertical slice for {configuration} configuration'
diff --git a/unit-tests/recovery_tests/test_reversible_recovery.py b/unit-tests/recovery_tests/test_reversible_recovery.py
index 721b0158..474e9ba0 100644
--- a/unit-tests/recovery_tests/test_reversible_recovery.py
+++ b/unit-tests/recovery_tests/test_reversible_recovery.py
@@ -9,8 +9,8 @@
 """
 
 from firedrake import (IntervalMesh, CubedSphereMesh, IcosahedralSphereMesh,
-                       SpatialCoordinate, FunctionSpace,
-                       Function, norm, errornorm, as_vector)
+                       SpatialCoordinate, FunctionSpace, Interpolator,
+                       Projector, Function, norm, errornorm, as_vector)
 from gusto import *
 import numpy as np
 import pytest
diff --git a/unit-tests/test_conservative_projection.py b/unit-tests/test_conservative_projection.py
new file mode 100644
index 00000000..bb47c891
--- /dev/null
+++ b/unit-tests/test_conservative_projection.py
@@ -0,0 +1,64 @@
+"""
+This tests the ConservativeProjector object, by projecting a mixing ratio from
+DG1 to DG0, relative to different density fields, and checking that the tracer
+mass is conserved.
+"""
+
+from firedrake import (UnitSquareMesh, FunctionSpace, Constant,
+                       Function, assemble, dx, sin, SpatialCoordinate)
+from gusto import ConservativeProjector
+import pytest
+
+
+@pytest.mark.parametrize("projection", ["discontinuous", "continuous"])
+def test_conservative_projection(projection):
+
+    # Set up mesh on plane
+    mesh = UnitSquareMesh(3, 3)
+
+    # Function spaces and functions
+    DG0 = FunctionSpace(mesh, "DG", 0)
+    DG1 = FunctionSpace(mesh, "DG", 1)
+
+    rho_DG0 = Function(DG0)
+    rho_DG1 = Function(DG1)
+    m_DG1 = Function(DG1)
+
+    if projection == "continuous":
+        CG1 = FunctionSpace(mesh, "CG", 1)
+        m_CG1 = Function(CG1)
+    else:
+        m_DG0 = Function(DG0)
+
+    # Projector object
+    if projection == "continuous":
+        projector = ConservativeProjector(rho_DG1, rho_DG0, m_DG1, m_CG1,
+                                          subtract_mean=True)
+    else:
+        projector = ConservativeProjector(rho_DG1, rho_DG0, m_DG1, m_DG0)
+
+    # Initial conditions
+    x, y = SpatialCoordinate(mesh)
+
+    rho_expr = Constant(1.0) + 0.5*x*y**2
+    m_expr = Constant(2.0) + 0.6*sin(x)
+
+    rho_DG1.interpolate(rho_expr)
+    m_DG1.interpolate(m_expr)
+    rho_DG0.project(rho_DG1)
+
+    # Test projection
+    projector.project()
+
+    tol = 1e-14
+    mass_DG1 = assemble(rho_DG1*m_DG1*dx)
+
+    if projection == "continuous":
+        mass_CG1 = assemble(rho_DG0*m_CG1*dx)
+
+        assert abs(mass_CG1 - mass_DG1) < tol, "continuous projection is not conservative"
+
+    else:
+        mass_DG0 = assemble(rho_DG0*m_DG0*dx)
+
+        assert abs(mass_DG0 - mass_DG1) < tol, "discontinuous projection is not conservative"