diff --git a/docs/source/usersguide/decay_sources.rst b/docs/source/usersguide/decay_sources.rst
index 398680e7464..21981fdaa82 100644
--- a/docs/source/usersguide/decay_sources.rst
+++ b/docs/source/usersguide/decay_sources.rst
@@ -190,6 +190,25 @@ we would run::
 
     r2s.run(timesteps, source_rates, mat_vol_kwargs={'n_samples': 10_000_000})
 
+It is also possible to use multiple meshes by passing a list of meshes instead
+of a single mesh. This can be useful, for example, when different regions of the
+model require different mesh resolutions. The meshes are assumed to be
+**non-overlapping**; each element--material combination across all meshes is
+treated as an independent activation region, and all meshes are handled in a
+single neutron transport solve. For example::
+
+    # Fine mesh near the activation target
+    mesh_fine = openmc.RegularMesh()
+    mesh_fine.dimension = (10, 10, 10)
+    ...
+
+    # Coarse mesh for the surrounding region
+    mesh_coarse = openmc.RegularMesh()
+    mesh_coarse.dimension = (5, 5, 5)
+    ...
+
+    r2s = openmc.deplete.R2SManager(model, [mesh_fine, mesh_coarse])
+
 Direct 1-Step (D1S) Calculations
 ================================
 
diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py
index e6e2dbce24c..6e82e5ebe98 100644
--- a/openmc/deplete/microxs.py
+++ b/openmc/deplete/microxs.py
@@ -36,7 +36,8 @@
     Sequence[openmc.Cell],
     Sequence[openmc.Universe],
     openmc.MeshBase,
-    openmc.Filter
+    openmc.Filter,
+    Sequence[openmc.Filter]
 ]
 
 
@@ -69,8 +70,12 @@ def get_microxs_and_flux(
     ----------
     model : openmc.Model
         OpenMC model object. Must contain geometry, materials, and settings.
-    domains : list of openmc.Material or openmc.Cell or openmc.Universe, or openmc.MeshBase, or openmc.Filter
+    domains : list of openmc.Material or openmc.Cell or openmc.Universe, or openmc.MeshBase, or openmc.Filter, or list of openmc.Filter
         Domains in which to tally reaction rates, or a spatial tally filter.
+        A list of filters can be provided to create one set of tallies per
+        filter (e.g., one :class:`~openmc.MeshMaterialFilter` per mesh) that
+        are all evaluated in a single transport solve. Results are
+        concatenated across all filters in order.
     nuclides : list of str
         Nuclides to get cross sections for. If not specified, all burnable
         nuclides from the depletion chain file are used.
@@ -142,26 +147,24 @@ def get_microxs_and_flux(
     else:
         energy_filter = openmc.EnergyFilter(energies)
 
+    # Build list of domain filters
     if isinstance(domains, openmc.Filter):
-        domain_filter = domains
+        domain_filters = [domains]
     elif isinstance(domains, openmc.MeshBase):
-        domain_filter = openmc.MeshFilter(domains)
+        domain_filters = [openmc.MeshFilter(domains)]
+    elif isinstance(domains, Sequence) and len(domains) > 0 and \
+            isinstance(domains[0], openmc.Filter):
+        domain_filters = list(domains)
     elif isinstance(domains[0], openmc.Material):
-        domain_filter = openmc.MaterialFilter(domains)
+        domain_filters = [openmc.MaterialFilter(domains)]
     elif isinstance(domains[0], openmc.Cell):
-        domain_filter = openmc.CellFilter(domains)
+        domain_filters = [openmc.CellFilter(domains)]
     elif isinstance(domains[0], openmc.Universe):
-        domain_filter = openmc.UniverseFilter(domains)
+        domain_filters = [openmc.UniverseFilter(domains)]
     else:
         raise ValueError(f"Unsupported domain type: {type(domains[0])}")
 
-    flux_tally = openmc.Tally(name='MicroXS flux')
-    flux_tally.filters = [domain_filter, energy_filter]
-    flux_tally.scores = ['flux']
-    model.tallies = [flux_tally]
-
-    # Prepare reaction-rate tally for 'direct' or subset for 'flux' with opts
-    rr_tally = None
+    # Prepare reaction-rate nuclides/reactions
     rr_nuclides: list[str] = []
     rr_reactions: list[str] = []
     if reaction_rate_mode == 'direct':
@@ -177,20 +180,33 @@ def get_microxs_and_flux(
         if rr_reactions:
             rr_reactions = [r for r in rr_reactions if r in set(reactions)]
 
-    # Only construct tally if both lists are non-empty
-    if rr_nuclides and rr_reactions:
-        rr_tally = openmc.Tally(name='MicroXS RR')
-        # Use 1-group energy filter for RR in flux mode
-        if reaction_rate_mode == 'flux':
-            rr_energy_filter = openmc.EnergyFilter(
-                [energy_filter.values[0], energy_filter.values[-1]])
-        else:
-            rr_energy_filter = energy_filter
-        rr_tally.filters = [domain_filter, rr_energy_filter]
-        rr_tally.nuclides = rr_nuclides
-        rr_tally.multiply_density = False
-        rr_tally.scores = rr_reactions
-        model.tallies.append(rr_tally)
+    # Use 1-group energy filter for RR in flux mode
+    has_rr = bool(rr_nuclides and rr_reactions)
+    if has_rr and reaction_rate_mode == 'flux':
+        rr_energy_filter = openmc.EnergyFilter(
+            [energy_filter.values[0], energy_filter.values[-1]])
+    else:
+        rr_energy_filter = energy_filter
+
+    # Create one flux tally (and optionally one RR tally) per domain filter.
+    flux_tallies = []
+    rr_tallies = []
+    model.tallies = []
+    for i, domain_filter in enumerate(domain_filters):
+        flux_tally = openmc.Tally(name=f'MicroXS flux {i}')
+        flux_tally.filters = [domain_filter, energy_filter]
+        flux_tally.scores = ['flux']
+        model.tallies.append(flux_tally)
+        flux_tallies.append(flux_tally)
+
+        if has_rr:
+            rr_tally = openmc.Tally(name=f'MicroXS RR {i}')
+            rr_tally.filters = [domain_filter, rr_energy_filter]
+            rr_tally.nuclides = rr_nuclides
+            rr_tally.multiply_density = False
+            rr_tally.scores = rr_reactions
+            model.tallies.append(rr_tally)
+            rr_tallies.append(rr_tally)
 
     if openmc.lib.is_initialized:
         openmc.lib.finalize()
@@ -227,40 +243,40 @@ def get_microxs_and_flux(
 
         # Read in tally results (on all ranks)
         with StatePoint(statepoint_path) as sp:
-            if rr_tally is not None:
-                rr_tally = sp.tallies[rr_tally.id]
-                rr_tally._read_results()
-            flux_tally = sp.tallies[flux_tally.id]
-            flux_tally._read_results()
-
-    # Get flux values and make energy groups last dimension
-    flux = flux_tally.get_reshaped_data()  # (domains, groups, 1, 1)
-    flux = np.moveaxis(flux, 1, -1)  # (domains, 1, 1, groups)
-
-    # Create list where each item corresponds to one domain
-    fluxes = list(flux.squeeze((1, 2)))
-
-    # If we built a reaction-rate tally, compute microscopic cross sections
-    if rr_tally is not None:
-        # Get reaction rates
-        reaction_rates = rr_tally.get_reshaped_data()  # (domains, groups, nuclides, reactions)
-
-        # Make energy groups last dimension
-        reaction_rates = np.moveaxis(reaction_rates, 1, -1)  # (domains, nuclides, reactions, groups)
-
-        # If RR is 1-group, sum flux over groups
-        if reaction_rate_mode == "flux":
-            flux = flux.sum(axis=-1, keepdims=True)  # (domains, 1, 1, 1)
-
-        # Divide RR by flux to get microscopic cross sections. The indexing
-        # ensures that only non-zero flux values are used, and broadcasting is
-        # applied to align the shapes of reaction_rates and flux for division.
-        xs = np.zeros_like(reaction_rates)  # (domains, nuclides, reactions, groups)
-        d, _, _, g = np.nonzero(flux)
-        xs[d, ..., g] = reaction_rates[d, ..., g] / flux[d, :, :, g]
-
-        # Create lists where each item corresponds to one domain
-        direct_micros = [MicroXS(xs_i, rr_nuclides, rr_reactions) for xs_i in xs]
+            for i in range(len(flux_tallies)):
+                flux_tallies[i] = sp.tallies[flux_tallies[i].id]
+                flux_tallies[i]._read_results()
+                if rr_tallies:
+                    rr_tallies[i] = sp.tallies[rr_tallies[i].id]
+                    rr_tallies[i]._read_results()
+
+    # Concatenate results across all domain filters
+    fluxes = []
+    all_flux_arrays = []
+    for flux_tally in flux_tallies:
+        # Get flux values and make energy groups last dimension
+        flux = flux_tally.get_reshaped_data()  # (domains, groups, 1, 1)
+        flux = np.moveaxis(flux, 1, -1)  # (domains, 1, 1, groups)
+        all_flux_arrays.append(flux)
+        fluxes.extend(flux.squeeze((1, 2)))
+
+    # If we built reaction-rate tallies, compute microscopic cross sections
+    if rr_tallies:
+        direct_micros = []
+        for flux_arr, rr_tally in zip(all_flux_arrays, rr_tallies):
+            flux = flux_arr
+            reaction_rates = rr_tally.get_reshaped_data()
+            reaction_rates = np.moveaxis(reaction_rates, 1, -1)
+
+            # If RR is 1-group, sum flux over groups
+            if reaction_rate_mode == "flux":
+                flux = flux.sum(axis=-1, keepdims=True)
+
+            xs = np.zeros_like(reaction_rates)
+            d, _, _, g = np.nonzero(flux)
+            xs[d, ..., g] = reaction_rates[d, ..., g] / flux[d, :, :, g]
+            direct_micros.extend(
+                MicroXS(xs_i, rr_nuclides, rr_reactions) for xs_i in xs)
 
     # If using flux mode, compute flux-collapsed microscopic XS
     if reaction_rate_mode == 'flux':
@@ -273,9 +289,9 @@ def get_microxs_and_flux(
         ) for flux_i in fluxes]
 
     # Decide which micros to use and merge if needed
-    if reaction_rate_mode == 'flux' and rr_tally is not None:
+    if reaction_rate_mode == 'flux' and rr_tallies:
         micros = [m1.merge(m2) for m1, m2 in zip(flux_micros, direct_micros)]
-    elif rr_tally is not None:
+    elif rr_tallies:
         micros = direct_micros
     else:
         micros = flux_micros
diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py
index 57bbe437ffb..10d7fdd502a 100644
--- a/openmc/deplete/r2s.py
+++ b/openmc/deplete/r2s.py
@@ -13,11 +13,11 @@
 from ..checkvalue import PathLike
 from ..mpi import comm
 from openmc.lib import TemporarySession
-from openmc.utility_funcs import change_directory
 
 
 def get_activation_materials(
-    model: openmc.Model, mmv: openmc.MeshMaterialVolumes
+    model: openmc.Model,
+    mmv_list: list[openmc.MeshMaterialVolumes]
 ) -> openmc.Materials:
     """Get a list of activation materials for each mesh element/material.
 
@@ -31,35 +31,35 @@ def get_activation_materials(
     ----------
     model : openmc.Model
         The full model containing the geometry and materials.
-    mmv : openmc.MeshMaterialVolumes
-        The mesh material volumes object containing the materials and their
-        volumes for each mesh element.
+    mmv_list : list of openmc.MeshMaterialVolumes
+        List of mesh material volumes objects, one per mesh, containing the
+        materials and their volumes for each mesh element.
 
     Returns
     -------
     openmc.Materials
         A list of materials, each corresponding to a unique mesh element and
-        material combination.
+        material combination across all meshes.
 
     """
-    # Get the material ID, volume, and element index for each element-material
-    # combination
-    mat_ids = mmv._materials[mmv._materials > -1]
-    volumes = mmv._volumes[mmv._materials > -1]
-    elems, _ = np.where(mmv._materials > -1)
-
     # Get all materials in the model
     material_dict = model._get_all_materials()
 
     # Create a new activation material for each element-material combination
+    # across all meshes
     materials = openmc.Materials()
-    for elem, mat_id, vol in zip(elems, mat_ids, volumes):
-        mat = material_dict[mat_id]
-        new_mat = mat.clone()
-        new_mat.depletable = True
-        new_mat.name = f'Element {elem}, Material {mat_id}'
-        new_mat.volume = vol
-        materials.append(new_mat)
+    for mesh_idx, mmv in enumerate(mmv_list):
+        mat_ids = mmv._materials[mmv._materials > -1]
+        volumes = mmv._volumes[mmv._materials > -1]
+        elems, _ = np.where(mmv._materials > -1)
+
+        for elem, mat_id, vol in zip(elems, mat_ids, volumes):
+            mat = material_dict[mat_id]
+            new_mat = mat.clone()
+            new_mat.depletable = True
+            new_mat.name = f'Mesh {mesh_idx}, Element {elem}, Material {mat_id}'
+            new_mat.volume = vol
+            materials.append(new_mat)
 
     return materials
 
@@ -70,7 +70,9 @@ class R2SManager:
     This class is responsible for managing the materials and sources needed for
     mesh-based or cell-based R2S calculations. It provides methods to get
     activation materials and decay photon sources based on the mesh/cells and
-    materials in the OpenMC model.
+    materials in the OpenMC model. Multiple meshes can be specified as domains,
+    in which case each element--material combination of each mesh is treated as
+    an activation region (meshes are assumed to be non-overlapping).
 
     This class supports the use of a different models for the neutron and photon
     transport calculation. However, for cell-based calculations, it assumes that
@@ -83,17 +85,20 @@ class R2SManager:
     ----------
     neutron_model : openmc.Model
         The OpenMC model to use for neutron transport.
-    domains : openmc.MeshBase or Sequence[openmc.Cell]
-        The mesh or a sequence of cells that represent the spatial units over
-        which the R2S calculation will be performed.
+    domains : openmc.MeshBase or Sequence[openmc.MeshBase] or Sequence[openmc.Cell]
+        The mesh(es) or a sequence of cells that represent the spatial units
+        over which the R2S calculation will be performed. When a single
+        :class:`~openmc.MeshBase` or a sequence of meshes is given, each
+        element--material combination across all meshes is treated as an
+        activation region.
     photon_model : openmc.Model, optional
         The OpenMC model to use for photon transport calculations. If None, a
         shallow copy of the neutron_model will be created and used.
 
     Attributes
     ----------
-    domains : openmc.MeshBase or Sequence[openmc.Cell]
-        The mesh or a sequence of cells that represent the spatial units over
+    domains : list of openmc.MeshBase or Sequence[openmc.Cell]
+        The meshes or a sequence of cells that represent the spatial units over
         which the R2S calculation will be performed.
     neutron_model : openmc.Model
         The OpenMC model used for neutron transport.
@@ -101,7 +106,7 @@ class R2SManager:
         The OpenMC model used for photon transport calculations.
     method : {'mesh-based', 'cell-based'}
         Indicates whether the R2S calculation uses mesh elements ('mesh-based')
-        as the spatial discetization or a list of a cells ('cell-based').
+        as the spatial discretization or a list of cells ('cell-based').
     results : dict
         A dictionary that stores results from the R2S calculation.
 
@@ -109,7 +114,7 @@ class R2SManager:
     def __init__(
         self,
         neutron_model: openmc.Model,
-        domains: openmc.MeshBase | Sequence[openmc.Cell],
+        domains: openmc.MeshBase | Sequence[openmc.MeshBase] | Sequence[openmc.Cell],
         photon_model: openmc.Model | None = None,
     ):
         self.neutron_model = neutron_model
@@ -126,9 +131,14 @@ def __init__(
             self.photon_model = photon_model
         if isinstance(domains, openmc.MeshBase):
             self.method = 'mesh-based'
+            self.domains = [domains]
+        elif isinstance(domains, Sequence) and len(domains) > 0 and \
+                isinstance(domains[0], openmc.MeshBase):
+            self.method = 'mesh-based'
+            self.domains = list(domains)
         else:
             self.method = 'cell-based'
-        self.domains = domains
+            self.domains = list(domains)
         self.results = {}
 
     def run(
@@ -243,11 +253,13 @@ def step1_neutron_transport(
     ):
         """Run the neutron transport step.
 
-        This step computes the material volume fractions on the mesh, creates a
-        mesh-material filter, and retrieves the fluxes and microscopic cross
-        sections for each mesh/material combination. This step will populate the
-        'fluxes' and 'micros' keys in the results dictionary. For a mesh-based
-        calculation, it will also populate the 'mesh_material_volumes' key.
+        This step computes the material volume fractions on each mesh, creates
+        mesh-material filters, and retrieves the fluxes and microscopic cross
+        sections for each mesh/material combination via a single transport
+        solve. This step will populate the 'fluxes' and 'micros' keys in the
+        results dictionary. For a mesh-based calculation, it will also populate
+        the 'mesh_material_volumes' key (a list of
+        :class:`~openmc.MeshMaterialVolumes`, one per mesh).
 
         Parameters
         ----------
@@ -266,19 +278,28 @@ def step1_neutron_transport(
         output_dir.mkdir(parents=True, exist_ok=True)
 
         if self.method == 'mesh-based':
-            # Compute material volume fractions on the mesh
+            # Compute material volume fractions on each mesh
             if mat_vol_kwargs is None:
                 mat_vol_kwargs = {}
             mat_vol_kwargs.setdefault('bounding_boxes', True)
-            self.results['mesh_material_volumes'] = mmv = comm.bcast(
-                self.domains.material_volumes(self.neutron_model, **mat_vol_kwargs))
 
-            # Save results to file
-            if comm.rank == 0:
-                mmv.save(output_dir / 'mesh_material_volumes.npz')
+            mmv_list = []
+            domain_filters = []
+            for i, mesh in enumerate(self.domains):
+                mmv = comm.bcast(
+                    mesh.material_volumes(self.neutron_model, **mat_vol_kwargs))
+                mmv_list.append(mmv)
+
+                # Save results to file
+                if comm.rank == 0:
+                    mmv.save(output_dir / f'mesh_material_volumes_{i}.npz')
+
+                # Create mesh-material filter for this mesh
+                domain_filters.append(
+                    openmc.MeshMaterialFilter.from_volumes(mesh, mmv))
 
-            # Create mesh-material filter based on what combos were found
-            domains = openmc.MeshMaterialFilter.from_volumes(self.domains, mmv)
+            self.results['mesh_material_volumes'] = mmv_list
+            domains = domain_filters
         else:
             domains: Sequence[openmc.Cell] = self.domains
 
@@ -357,8 +378,9 @@ def step2_activation(
 
         if self.method == 'mesh-based':
             # Get unique material for each (mesh, material) combination
-            mmv = self.results['mesh_material_volumes']
-            self.results['activation_materials'] = get_activation_materials(self.neutron_model, mmv)
+            mmv_list = self.results['mesh_material_volumes']
+            self.results['activation_materials'] = get_activation_materials(
+                self.neutron_model, mmv_list)
         else:
             # Create unique material for each cell
             activation_mats = openmc.Materials()
@@ -468,12 +490,20 @@ def step3_photon_transport(
         # photon model if it is different from the neutron model to account for
         # potential material changes
         if self.method == 'mesh-based' and different_photon_model:
-            self.results['mesh_material_volumes_photon'] = photon_mmv = comm.bcast(
-                self.domains.material_volumes(self.photon_model, **mat_vol_kwargs))
+            if mat_vol_kwargs is None:
+                mat_vol_kwargs = {}
+            photon_mmv_list = []
+            for i, mesh in enumerate(self.domains):
+                photon_mmv = comm.bcast(
+                    mesh.material_volumes(self.photon_model, **mat_vol_kwargs))
+                photon_mmv_list.append(photon_mmv)
+
+                # Save photon MMV results to file
+                if comm.rank == 0:
+                    photon_mmv.save(
+                        output_dir / f'mesh_material_volumes_{i}.npz')
 
-            # Save photon MMV results to file
-            if comm.rank == 0:
-                photon_mmv.save(output_dir / 'mesh_material_volumes.npz')
+            self.results['mesh_material_volumes_photon'] = photon_mmv_list
 
         if comm.rank == 0:
             tally_ids = [tally.id for tally in self.photon_model.tallies]
@@ -543,7 +573,7 @@ def get_decay_photon_source_mesh(
     ) -> list[openmc.IndependentSource]:
         """Create decay photon source for a mesh-based calculation.
 
-        For each mesh element-material combination, an
+        For each mesh element-material combination across all meshes, an
         :class:`~openmc.IndependentSource` is created with a
         :class:`~openmc.stats.Box` spatial distribution based on the bounding
         box of the material within the mesh element. A material constraint is
@@ -575,52 +605,56 @@ def get_decay_photon_source_mesh(
         index_mat = 0
 
         # Get various results from previous steps
-        mat_vols = self.results['mesh_material_volumes']
+        mmv_list = self.results['mesh_material_volumes']
         materials = self.results['activation_materials']
         results = self.results['depletion_results']
-        photon_mat_vols = self.results.get('mesh_material_volumes_photon')
-
-        # Total number of mesh elements
-        n_elements = mat_vols.num_elements
-
-        for index_elem in range(n_elements):
-            # Determine which materials exist in the photon model for this element
-            if photon_mat_vols is not None:
-                photon_materials = {
-                    mat_id
-                    for mat_id, _ in photon_mat_vols.by_element(index_elem)
-                    if mat_id is not None
-                }
-
-            for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True):
-                # Skip void volume
-                if mat_id is None:
-                    continue
-
-                # Skip if this material doesn't exist in photon model
-                if photon_mat_vols is not None and mat_id not in photon_materials:
-                    index_mat += 1
-                    continue
-
-                # Get activated material composition
-                original_mat = materials[index_mat]
-                activated_mat = results[time_index].get_material(str(original_mat.id))
-
-                # Create decay photon source
-                energy = activated_mat.get_decay_photon_energy()
-                if energy is not None:
-                    strength = energy.integral()
-                    space = openmc.stats.Box(*bbox)
-                    sources.append(openmc.IndependentSource(
-                        space=space,
-                        energy=energy,
-                        particle='photon',
-                        strength=strength,
-                        constraints={'domains': [mat_dict[mat_id]]}
-                    ))
+        photon_mmv_list = self.results.get('mesh_material_volumes_photon')
+
+        for mesh_idx, mat_vols in enumerate(mmv_list):
+            photon_mat_vols = photon_mmv_list[mesh_idx] \
+                if photon_mmv_list is not None else None
+
+            # Total number of mesh elements for this mesh
+            n_elements = mat_vols.num_elements
+
+            for index_elem in range(n_elements):
+                # Determine which materials exist in the photon model for this element
+                if photon_mat_vols is not None:
+                    photon_materials = {
+                        mat_id
+                        for mat_id, _ in photon_mat_vols.by_element(index_elem)
+                        if mat_id is not None
+                    }
+
+                for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True):
+                    # Skip void volume
+                    if mat_id is None:
+                        continue
+
+                    # Skip if this material doesn't exist in photon model
+                    if photon_mat_vols is not None and mat_id not in photon_materials:
+                        index_mat += 1
+                        continue
+
+                    # Get activated material composition
+                    original_mat = materials[index_mat]
+                    activated_mat = results[time_index].get_material(str(original_mat.id))
 
-                # Increment index of activated material
-                index_mat += 1
+                    # Create decay photon source
+                    energy = activated_mat.get_decay_photon_energy()
+                    if energy is not None:
+                        strength = energy.integral()
+                        space = openmc.stats.Box(*bbox)
+                        sources.append(openmc.IndependentSource(
+                            space=space,
+                            energy=energy,
+                            particle='photon',
+                            strength=strength,
+                            constraints={'domains': [mat_dict[mat_id]]}
+                        ))
+
+                    # Increment index of activated material
+                    index_mat += 1
 
         return sources
 
@@ -638,10 +672,13 @@ def load_results(self, path: PathLike):
         # Load neutron transport results
         neutron_dir = path / 'neutron_transport'
         if self.method == 'mesh-based':
-            mmv_file = neutron_dir / 'mesh_material_volumes.npz'
-            if mmv_file.exists():
-                self.results['mesh_material_volumes'] = \
-                    openmc.MeshMaterialVolumes.from_npz(mmv_file)
+            mmv_files = sorted(neutron_dir.glob('mesh_material_volumes*.npz'),
+                               key=lambda p: int(p.stem.split('_')[-1])
+                               if p.stem[-1].isdigit() else 0)
+            if mmv_files:
+                self.results['mesh_material_volumes'] = [
+                    openmc.MeshMaterialVolumes.from_npz(f) for f in mmv_files
+                ]
         fluxes_file = neutron_dir / 'fluxes.npy'
         if fluxes_file.exists():
             self.results['fluxes'] = list(np.load(fluxes_file, allow_pickle=True))
@@ -665,10 +702,15 @@ def load_results(self, path: PathLike):
 
         # Load photon mesh material volumes if they exist (for mesh-based calculations)
         if self.method == 'mesh-based':
-            photon_mmv_file = photon_dir / 'mesh_material_volumes.npz'
-            if photon_mmv_file.exists():
-                self.results['mesh_material_volumes_photon'] = \
-                    openmc.MeshMaterialVolumes.from_npz(photon_mmv_file)
+            photon_mmv_files = sorted(
+                photon_dir.glob('mesh_material_volumes*.npz'),
+                key=lambda p: int(p.stem.split('_')[-1])
+                if p.stem[-1].isdigit() else 0)
+            if photon_mmv_files:
+                self.results['mesh_material_volumes_photon'] = [
+                    openmc.MeshMaterialVolumes.from_npz(f)
+                    for f in photon_mmv_files
+                ]
 
         # Load tally IDs from JSON file
         tally_ids_path = photon_dir / 'tally_ids.json'
diff --git a/tests/unit_tests/test_deplete_microxs.py b/tests/unit_tests/test_deplete_microxs.py
index 340ba616323..26529e6ce96 100644
--- a/tests/unit_tests/test_deplete_microxs.py
+++ b/tests/unit_tests/test_deplete_microxs.py
@@ -165,11 +165,11 @@ def capture_run(**kwargs):
 
     # Check that both tallies were created with the expected properties
     tally_names = [t.name for t in captured['tallies']]
-    assert 'MicroXS flux' in tally_names
-    assert 'MicroXS RR' in tally_names
+    assert 'MicroXS flux 0' in tally_names
+    assert 'MicroXS RR 0' in tally_names
 
     # Check that the RR tally has the expected nuclides and reactions
-    rr = next(t for t in captured['tallies'] if t.name == 'MicroXS RR')
+    rr = next(t for t in captured['tallies'] if t.name == 'MicroXS RR 0')
     assert rr.nuclides == ['U235']
     assert rr.scores == ['fission']
 
diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py
index a94f85c8c08..266bd49763a 100644
--- a/tests/unit_tests/test_r2s.py
+++ b/tests/unit_tests/test_r2s.py
@@ -6,7 +6,7 @@
 
 
 @pytest.fixture
-def simple_model_and_mesh(tmp_path):
+def simple_model_and_mesh():
     # Define two materials: water and Ni
     h2o = openmc.Material()
     h2o.add_nuclide("H1", 2.0)
@@ -68,7 +68,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
     nt = Path(outdir) / 'neutron_transport'
     assert (nt / 'fluxes.npy').exists()
     assert (nt / 'micros.h5').exists()
-    assert (nt / 'mesh_material_volumes.npz').exists()
+    assert (nt / 'mesh_material_volumes_0.npz').exists()
     act = Path(outdir) / 'activation'
     assert (act / 'depletion_results.h5').exists()
     pt = Path(outdir) / 'photon_transport'
@@ -78,7 +78,8 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
     # Basic results structure checks
     assert len(r2s.results['fluxes']) == 2
     assert len(r2s.results['micros']) == 2
-    assert len(r2s.results['mesh_material_volumes']) == 2
+    assert len(r2s.results['mesh_material_volumes']) == 1
+    assert len(r2s.results['mesh_material_volumes'][0]) == 2
     assert len(r2s.results['activation_materials']) == 2
     assert len(r2s.results['depletion_results']) == 2
 
@@ -93,11 +94,77 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
     r2s_loaded.load_results(outdir)
     assert len(r2s_loaded.results['fluxes']) == 2
     assert len(r2s_loaded.results['micros']) == 2
-    assert len(r2s_loaded.results['mesh_material_volumes']) == 2
+    assert len(r2s_loaded.results['mesh_material_volumes']) == 1
+    assert len(r2s_loaded.results['mesh_material_volumes'][0]) == 2
     assert len(r2s_loaded.results['activation_materials']) == 2
     assert len(r2s_loaded.results['depletion_results']) == 2
 
 
+def test_r2s_multi_mesh(simple_model_and_mesh, tmp_path):
+    model, _, _ = simple_model_and_mesh
+
+    # Two 1x1x1 meshes that together cover the full domain, split along y.
+    # Each mesh element spans the full x range [-10, 10], crossing the x=0
+    # material boundary, so both meshes contain both materials within their
+    # single element.
+    mesh1 = openmc.RegularMesh()
+    mesh1.lower_left = (-10.0, -10.0, -10.0)
+    mesh1.upper_right = (10.0, 0.0, 10.0)
+    mesh1.dimension = (1, 1, 1)
+    mesh2 = openmc.RegularMesh()
+    mesh2.lower_left = (-10.0, 0.0, -10.0)
+    mesh2.upper_right = (10.0, 10.0, 10.0)
+    mesh2.dimension = (1, 1, 1)
+
+    r2s = R2SManager(model, [mesh1, mesh2])
+    chain = Chain.from_xml(Path(__file__).parents[1] / "chain_ni.xml")
+
+    outdir = r2s.run(
+        timesteps=[(1.0, 'd')],
+        source_rates=[1.0],
+        photon_time_indices=[1],
+        output_dir=tmp_path,
+        chain_file=chain,
+    )
+
+    # Check that per-mesh MMV files were written
+    nt = Path(outdir) / 'neutron_transport'
+    assert (nt / 'fluxes.npy').exists()
+    assert (nt / 'micros.h5').exists()
+    assert (nt / 'mesh_material_volumes_0.npz').exists()
+    assert (nt / 'mesh_material_volumes_1.npz').exists()
+    act = Path(outdir) / 'activation'
+    assert (act / 'depletion_results.h5').exists()
+    pt = Path(outdir) / 'photon_transport'
+    assert (pt / 'tally_ids.json').exists()
+    assert (pt / 'time_1' / 'statepoint.10.h5').exists()
+
+    # Two meshes, each with 1 element containing both materials →
+    # 2 element-material combinations per mesh, 4 total
+    assert len(r2s.results['mesh_material_volumes']) == 2
+    assert len(r2s.results['mesh_material_volumes'][0]) == 2
+    assert len(r2s.results['mesh_material_volumes'][1]) == 2
+    assert len(r2s.results['fluxes']) == 4
+    assert len(r2s.results['micros']) == 4
+    assert len(r2s.results['activation_materials']) == 4
+    assert len(r2s.results['depletion_results']) == 2
+
+    # Activation material names encode mesh index
+    amats = r2s.results['activation_materials']
+    assert all(m.depletable for m in amats)
+    assert any('Mesh 0' in m.name for m in amats)
+    assert any('Mesh 1' in m.name for m in amats)
+
+    # Check loading results
+    r2s_loaded = R2SManager(model, [mesh1, mesh2])
+    r2s_loaded.load_results(outdir)
+    assert len(r2s_loaded.results['mesh_material_volumes']) == 2
+    assert len(r2s_loaded.results['mesh_material_volumes'][0]) == 2
+    assert len(r2s_loaded.results['mesh_material_volumes'][1]) == 2
+    assert len(r2s_loaded.results['activation_materials']) == 4
+    assert len(r2s_loaded.results['depletion_results']) == 2
+
+
 def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
     model, (c1, c2), _ = simple_model_and_mesh