diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 78a08181a0..52274519ea 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -437,6 +437,34 @@ def __init__(self, child, integration_variable):
         if not isinstance(integration_variable, list):
             integration_variable = [integration_variable]
 
+        name = "integral"
+        for var in integration_variable:
+            if isinstance(var, pybamm.SpatialVariable):
+                # Check that child and integration_variable domains agree
+                if var.domain == child.domain:
+                    self._integration_domain = "primary"
+                elif (
+                    "secondary" in child.auxiliary_domains
+                    and var.domain == child.auxiliary_domains["secondary"]
+                ):
+                    self._integration_domain = "secondary"
+                elif (
+                    "tertiary" in child.auxiliary_domains
+                    and var.domain == child.auxiliary_domains["tertiary"]
+                ):
+                    self._integration_domain = "tertiary"
+                else:
+                    raise pybamm.DomainError(
+                        "integration_variable must be the same as child domain or "
+                        "auxiliary domain"
+                    )
+            elif not isinstance(var, pybamm.IndependentVariable):
+                raise ValueError(
+                    "integration_variable must be of type pybamm.IndependentVariable, "
+                    "not {}".format(type(var))
+                )
+            name += " d{}".format(var.name)
+
         # integral of a child takes the domain from auxiliary domain of the child
         if child.auxiliary_domains != {}:
             domain = child.auxiliary_domains["secondary"]
@@ -448,22 +476,6 @@ def __init__(self, child, integration_variable):
         else:
             domain = []
             auxiliary_domains = {}
-        name = "integral"
-        for var in integration_variable:
-            if isinstance(var, pybamm.SpatialVariable):
-                # Check that child and integration_variable domains agree
-                if child.domain != var.domain:
-                    raise pybamm.DomainError(
-                        "child and integration_variable must have the same domain"
-                    )
-            elif not isinstance(var, pybamm.IndependentVariable):
-                raise ValueError(
-                    """integration_variable must be of type pybamm.IndependentVariable,
-                           not {}""".format(
-                        type(var)
-                    )
-                )
-            name += " d{}".format(var.name)
 
         if any(isinstance(var, pybamm.SpatialVariable) for var in integration_variable):
             name += " {}".format(child.domain)
diff --git a/pybamm/models/submodels/particle/fickian_many_particles.py b/pybamm/models/submodels/particle/fickian_many_particles.py
index 3cb2a639ac..a2515c1491 100644
--- a/pybamm/models/submodels/particle/fickian_many_particles.py
+++ b/pybamm/models/submodels/particle/fickian_many_particles.py
@@ -30,9 +30,8 @@ def get_fundamental_variables(self):
         elif self.domain == "Positive":
             c_s = pybamm.standard_variables.c_s_p
 
-        # TODO: implement c_s_xav for Fickian many particles (tricky because this
-        # requires averaging a secondary domain)
-        variables = self._get_standard_concentration_variables(c_s, c_s)
+        c_s_xav = pybamm.x_average(c_s)
+        variables = self._get_standard_concentration_variables(c_s, c_s_xav)
 
         return variables
 
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 030deaa73d..4551eeb544 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -135,7 +135,7 @@ def test_integral(self):
         inta_sec = pybamm.Integral(a_sec, x)
         self.assertEqual(inta_sec.domain, ["current collector"])
         self.assertEqual(inta_sec.auxiliary_domains, {})
-        # space integral with secondary domain
+        # space integral with tertiary domain
         a_tert = pybamm.Symbol(
             "a",
             domain=["negative electrode"],
@@ -151,6 +151,22 @@ def test_integral(self):
             inta_tert.auxiliary_domains, {"secondary": ["some extra domain"]}
         )
 
+        # space integral *in* secondary domain
+        y = pybamm.SpatialVariable("y", ["current collector"])
+        inta_tert_y = pybamm.Integral(a_tert, y)
+        self.assertEqual(inta_tert_y.domain, ["negative electrode"])
+        self.assertEqual(
+            inta_tert_y.auxiliary_domains, {"secondary": ["some extra domain"]}
+        )
+
+        # space integral *in* tertiary domain
+        z = pybamm.SpatialVariable("z", ["some extra  domain"])
+        inta_tert_z = pybamm.Integral(a_tert, z)
+        self.assertEqual(inta_tert_z.domain, ["negative electrode"])
+        self.assertEqual(
+            inta_tert_z.auxiliary_domains, {"secondary": ["current collector"]}
+        )
+
         # space integral over two variables
         b = pybamm.Symbol("b", domain=["current collector"])
         y = pybamm.SpatialVariable("y", ["current collector"])