diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index de5a745e..bb0b8388 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -79,10 +79,19 @@ class Discretisation(object):
         variables : iterable of :class:`pybamm.Variables`
             The variables for which to set slices
         """
-        y_slices = {variable.id: None for variable in variables}
+        # Unpack symbols in variables that are concatenations of variables
+        unpacked_variables = []
+        for symbol in variables:
+            if isinstance(symbol, pybamm.Concatenation):
+                unpacked_variables.extend([var for var in symbol.children])
+            else:
+                unpacked_variables.append(symbol)
+        # Set up y_slices
+        y_slices = {variable.id: None for variable in unpacked_variables}
         start = 0
         end = 0
-        for variable in variables:
+        # Iterate through unpacked variables, adding appropriate slices to y_slices
+        for variable in unpacked_variables:
             # If domain is empty then variable has size 1
             if variable.domain == []:
                 end += 1
@@ -333,7 +342,15 @@ class Discretisation(object):
             Discretised right-hand side equations
 
         """
-        ids = {v.id for v in var_eqn_dict.keys()}
+        # Unpack symbols in variables that are concatenations of variables
+        unpacked_variables = []
+        for symbol in var_eqn_dict.keys():
+            if isinstance(symbol, pybamm.Concatenation):
+                unpacked_variables.extend([var for var in symbol.children])
+            else:
+                unpacked_variables.append(symbol)
+        # Check keys from the given var_eqn_dict against self._y_slices
+        ids = {v.id for v in unpacked_variables}
         if ids != self._y_slices.keys():
             given_variable_names = [v.name for v in var_eqn_dict.keys()]
             raise pybamm.ModelError(
@@ -342,7 +359,7 @@ class Discretisation(object):
             )
 
         equations = list(var_eqn_dict.values())
-        slices = [self._y_slices[var.id] for var in var_eqn_dict.keys()]
+        slices = [self._y_slices[var.id] for var in unpacked_variables]
 
         # sort equations according to slices
         sorted_equations = [eq for _, eq in sorted(zip(slices, equations))]
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 7c484f00..60d0df17 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -232,13 +232,20 @@ class BaseModel(object):
         vars_in_algebraic_keys = set()
         vars_in_eqns = set()
         # Get all variables ids from rhs and algebraic keys and equations
+        # For equations we look through the whole expression tree.
+        # "Variables" can be Concatenations so we also have to look in the whole
+        # expression tree
         for var, eqn in self.rhs.items():
-            vars_in_rhs_keys.add(var.id)
+            vars_in_rhs_keys.update(
+                [x.id for x in var.pre_order() if isinstance(x, pybamm.Variable)]
+            )
             vars_in_eqns.update(
                 [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
             )
         for var, eqn in self.algebraic.items():
-            vars_in_algebraic_keys.add(var.id)
+            vars_in_algebraic_keys.update(
+                [x.id for x in var.pre_order() if isinstance(x, pybamm.Variable)]
+            )
             vars_in_eqns.update(
                 [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
             )

Allow Concatenation of Variables as a key in the model
This will make writing the model down easier, but requires some adaptations:
- Set y_slices for the children of the concatenation, not the concatenation itself
- Take gradient and divergence appropriately
- Reading initial and boundary conditions *should* work ok as-is

### Example

```python3
cn = Variable("conc", ["negative electrode"])
cs = Variable("conc", ["separator"])
cp = Variable("conc", ["positive electrode"])
c = Concatenation(cn, cs, cp)
jn = butler_volmer(cn)
jp = butler_volmer(cp)
j = Concatenation(jn, 0, jp)
model.rhs = {c: div(grad(c)) + sj} # avoids having to write this equation down 3 times
```
diff --git a/tests/test_discretisations/test_discretisation.py b/tests/test_discretisations/test_discretisation.py
index 7327f60b..01995f97 100644
--- a/tests/test_discretisations/test_discretisation.py
+++ b/tests/test_discretisations/test_discretisation.py
@@ -71,6 +71,29 @@ class TestDiscretise(unittest.TestCase):
         np.testing.assert_array_equal(y[disc._y_slices[d.id]], d_true)
         np.testing.assert_array_equal(y[disc._y_slices[jn.id]], jn_true)
 
+        # Variables with a concatenation
+        js = pybamm.Variable("js", domain=["separator"])
+        jp = pybamm.Variable("jp", domain=["positive electrode"])
+        j = pybamm.Concatenation(jn, js, jp)
+        variables = [c, d, j]
+        disc.set_variable_slices(variables)
+        self.assertEqual(
+            disc._y_slices,
+            {
+                c.id: slice(0, 100),
+                d.id: slice(100, 200),
+                jn.id: slice(200, 240),
+                js.id: slice(240, 265),
+                jp.id: slice(265, 300),
+            },
+        )
+        d_true = 4 * combined_submesh.nodes
+        jn_true = mesh["negative electrode"].nodes ** 3
+        y = np.concatenate([c_true, d_true, jn_true])
+        np.testing.assert_array_equal(y[disc._y_slices[c.id]], c_true)
+        np.testing.assert_array_equal(y[disc._y_slices[d.id]], d_true)
+        np.testing.assert_array_equal(y[disc._y_slices[jn.id]], jn_true)
+
     def test_process_symbol_base(self):
         # create discretisation
         mesh = get_mesh_for_testing()
@@ -437,6 +460,36 @@ class TestDiscretise(unittest.TestCase):
         with self.assertRaises(pybamm.ModelError):
             disc.process_model(model)
 
+    def test_process_model_concatenation(self):
+        # concatenation of variables as the key
+        cn = pybamm.Variable("c", domain=["negative electrode"])
+        cs = pybamm.Variable("c", domain=["separator"])
+        cp = pybamm.Variable("c", domain=["positive electrode"])
+        c = pybamm.Concatenation(cn, cs, cp)
+        N = pybamm.grad(c)
+        model = pybamm.BaseModel()
+        model.rhs = {c: pybamm.div(N)}
+        model.initial_conditions = {c: pybamm.Scalar(3)}
+
+        model.boundary_conditions = {N: {"left": 0, "right": 0}}
+        model.check_well_posedness()
+
+        # create discretisation
+        disc = get_discretisation_for_testing()
+        mesh = disc.mesh
+
+        combined_submesh = mesh.combine_submeshes(
+            "negative electrode", "separator", "positive electrode"
+        )
+
+        disc.process_model(model)
+        y0 = model.concatenated_initial_conditions
+        np.testing.assert_array_equal(y0, 3 * np.ones_like(combined_submesh.nodes))
+
+        # grad and div are identity operators here
+        np.testing.assert_array_equal(y0, model.concatenated_rhs.evaluate(None, y0))
+        model.check_well_posedness()
+
     def test_broadcast(self):
         whole_cell = ["negative electrode", "separator", "positive electrode"]
 

pip install -e .[dev,docs]
parse_log_pytest
pytest
flake8 pytest
apt-get update apt-get install -y python3-dev gfortran gcc libopenblas-dev liblapack-dev
3.6
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.13
anyio==3.6.2
anytree==2.8.0
argon2-cffi==21.3.0
argon2-cffi-bindings==21.2.0
async-generator==1.10
attrs @ file:///opt/conda/conda-bld/attrs_1642510447205/work
Babel==2.11.0
backcall==0.2.0
bleach==4.1.0
certifi==2021.5.30
cffi==1.15.1
charset-normalizer==2.0.12
comm==0.1.4
contextvars==2.4
cycler==0.11.0
dataclasses==0.8
decorator==5.1.1
defusedxml==0.7.1
docutils==0.17.1
entrypoints==0.4
flake8==5.0.4
guzzle-sphinx-theme==0.7.11
idna==3.10
imagesize==1.4.1
immutables==0.19
importlib-metadata==4.2.0
iniconfig @ file:///home/linux1/recipes/ci/iniconfig_1610983019677/work
ipykernel==5.5.6
ipython==7.16.3
ipython-genutils==0.2.0
ipywidgets==7.8.5
jedi==0.17.2
Jinja2==3.0.3
json5==0.9.16
jsonschema==3.2.0
jupyter==1.1.1
jupyter-client==7.1.2
jupyter-console==6.4.3
jupyter-core==4.9.2
jupyter-server==1.13.1
jupyterlab==3.2.9
jupyterlab-pygments==0.1.2
jupyterlab-server==2.10.3
jupyterlab_widgets==1.1.11
kiwisolver==1.3.1
MarkupSafe==2.0.1
matplotlib==3.3.4
mccabe==0.7.0
mistune==0.8.4
more-itertools @ file:///tmp/build/80754af9/more-itertools_1637733554872/work
nbclassic==0.3.5
nbclient==0.5.9
nbconvert==6.0.7
nbformat==5.1.3
nest-asyncio==1.6.0
notebook==6.4.10
numpy==1.19.5
packaging @ file:///tmp/build/80754af9/packaging_1637314298585/work
pandas==1.1.5
pandocfilters==1.5.1
parso==0.7.1
pexpect==4.9.0
pickleshare==0.7.5
Pillow==8.4.0
pluggy @ file:///tmp/build/80754af9/pluggy_1615976315926/work
prometheus-client==0.17.1
prompt-toolkit==3.0.36
ptyprocess==0.7.0
py @ file:///opt/conda/conda-bld/py_1644396412707/work
-e git+https://github.com/tinosulzer/PyBaMM.git@de25f11d8b1b51dbcb3c4f9e6b8503b6dcafe5e3#egg=pybamm
pycodestyle==2.9.1
pycparser==2.21
pyflakes==2.5.0
Pygments==2.14.0
pyparsing @ file:///tmp/build/80754af9/pyparsing_1635766073266/work
pyrsistent==0.18.0
pytest==6.2.4
python-dateutil==2.9.0.post0
pytz==2025.2
pyzmq==25.1.2
requests==2.27.1
scipy==1.5.4
Send2Trash==1.8.3
six==1.17.0
sniffio==1.2.0
snowballstemmer==2.2.0
Sphinx==4.3.2
sphinxcontrib-applehelp==1.0.2
sphinxcontrib-devhelp==1.0.2
sphinxcontrib-htmlhelp==2.0.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==1.0.3
sphinxcontrib-serializinghtml==1.1.5
terminado==0.12.1
testpath==0.6.0
toml @ file:///tmp/build/80754af9/toml_1616166611790/work
tornado==6.1
traitlets==4.3.3
typing_extensions @ file:///opt/conda/conda-bld/typing_extensions_1647553014482/work
urllib3==1.26.20
wcwidth==0.2.13
webencodings==0.5.1
websocket-client==1.3.1
widgetsnbextension==3.6.10
zipp @ file:///tmp/build/80754af9/zipp_1633618647012/work

name: PyBaMM
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - attrs=21.4.0=pyhd3eb1b0_0
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2021.5.30=py36h06a4308_0
  - importlib_metadata=4.8.1=hd3eb1b0_0
  - iniconfig=1.1.1=pyhd3eb1b0_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - more-itertools=8.12.0=pyhd3eb1b0_0
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - packaging=21.3=pyhd3eb1b0_0
  - pip=21.2.2=py36h06a4308_0
  - pluggy=0.13.1=py36h06a4308_0
  - py=1.11.0=pyhd3eb1b0_0
  - pyparsing=3.0.4=pyhd3eb1b0_0
  - pytest=6.2.4=py36h06a4308_2
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - toml=0.10.2=pyhd3eb1b0_0
  - typing_extensions=4.1.1=pyh06a4308_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zipp=3.6.0=pyhd3eb1b0_0
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.13
      - anyio==3.6.2
      - anytree==2.8.0
      - argon2-cffi==21.3.0
      - argon2-cffi-bindings==21.2.0
      - async-generator==1.10
      - babel==2.11.0
      - backcall==0.2.0
      - bleach==4.1.0
      - cffi==1.15.1
      - charset-normalizer==2.0.12
      - comm==0.1.4
      - contextvars==2.4
      - cycler==0.11.0
      - dataclasses==0.8
      - decorator==5.1.1
      - defusedxml==0.7.1
      - docutils==0.17.1
      - entrypoints==0.4
      - flake8==5.0.4
      - guzzle-sphinx-theme==0.7.11
      - idna==3.10
      - imagesize==1.4.1
      - immutables==0.19
      - importlib-metadata==4.2.0
      - ipykernel==5.5.6
      - ipython==7.16.3
      - ipython-genutils==0.2.0
      - ipywidgets==7.8.5
      - jedi==0.17.2
      - jinja2==3.0.3
      - json5==0.9.16
      - jsonschema==3.2.0
      - jupyter==1.1.1
      - jupyter-client==7.1.2
      - jupyter-console==6.4.3
      - jupyter-core==4.9.2
      - jupyter-server==1.13.1
      - jupyterlab==3.2.9
      - jupyterlab-pygments==0.1.2
      - jupyterlab-server==2.10.3
      - jupyterlab-widgets==1.1.11
      - kiwisolver==1.3.1
      - markupsafe==2.0.1
      - matplotlib==3.3.4
      - mccabe==0.7.0
      - mistune==0.8.4
      - nbclassic==0.3.5
      - nbclient==0.5.9
      - nbconvert==6.0.7
      - nbformat==5.1.3
      - nest-asyncio==1.6.0
      - notebook==6.4.10
      - numpy==1.19.5
      - pandas==1.1.5
      - pandocfilters==1.5.1
      - parso==0.7.1
      - pexpect==4.9.0
      - pickleshare==0.7.5
      - pillow==8.4.0
      - prometheus-client==0.17.1
      - prompt-toolkit==3.0.36
      - ptyprocess==0.7.0
      - pycodestyle==2.9.1
      - pycparser==2.21
      - pyflakes==2.5.0
      - pygments==2.14.0
      - pyrsistent==0.18.0
      - python-dateutil==2.9.0.post0
      - pytz==2025.2
      - pyzmq==25.1.2
      - requests==2.27.1
      - scipy==1.5.4
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.2.0
      - snowballstemmer==2.2.0
      - sphinx==4.3.2
      - sphinxcontrib-applehelp==1.0.2
      - sphinxcontrib-devhelp==1.0.2
      - sphinxcontrib-htmlhelp==2.0.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==1.0.3
      - sphinxcontrib-serializinghtml==1.1.5
      - terminado==0.12.1
      - testpath==0.6.0
      - tornado==6.1
      - traitlets==4.3.3
      - urllib3==1.26.20
      - wcwidth==0.2.13
      - webencodings==0.5.1
      - websocket-client==1.3.1
      - widgetsnbextension==3.6.10
prefix: /opt/conda/envs/PyBaMM

tests/test_discretisations/test_discretisation.py::TestDiscretise::test_discretise_slicing tests/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_concatenation
tests/test_discretisations/test_discretisation.py::TestDiscretise::test_broadcast tests/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenate_init tests/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenation tests/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenation_of_scalars tests/test_discretisations/test_discretisation.py::TestDiscretise::test_core_NotImplementedErrors tests/test_discretisations/test_discretisation.py::TestDiscretise::test_discretise_spatial_operator tests/test_discretisations/test_discretisation.py::TestDiscretise::test_discretise_spatial_variable tests/test_discretisations/test_discretisation.py::TestDiscretise::test_process_complex_expression tests/test_discretisations/test_discretisation.py::TestDiscretise::test_process_dict tests/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_dae tests/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_ode tests/test_discretisations/test_discretisation.py::TestDiscretise::test_process_symbol_base tests/test_discretisations/test_discretisation.py::TestDiscretise::test_process_variables_dict
pybamm/discretisations/discretisation.py pybamm/models/base_model.py
pybamm/discretisations/discretisation.py pybamm/models/base_model.py