diff --git a/smriprep/interfaces/surf.py b/smriprep/interfaces/surf.py
index 8b378f8925..bd91db0966 100644
--- a/smriprep/interfaces/surf.py
+++ b/smriprep/interfaces/surf.py
@@ -33,6 +33,8 @@
     SimpleInterface,
     File,
     isdefined,
+    InputMultiObject,
+    traits,
 )
 
 
@@ -114,6 +116,45 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class AggregateSurfacesInputSpec(TraitedSpec):
+    surfaces = InputMultiObject(File(exists=True), desc="Input surfaces")
+    morphometrics = InputMultiObject(File(exists=True), desc="Input morphometrics")
+
+
+class AggregateSurfacesOutputSpec(TraitedSpec):
+    pial = traits.List(File(), maxlen=2, desc="Pial surfaces")
+    white = traits.List(File(), maxlen=2, desc="White surfaces")
+    inflated = traits.List(File(), maxlen=2, desc="Inflated surfaces")
+    midthickness = traits.List(File(), maxlen=2, desc="Midthickness (or graymid) surfaces")
+    thickness = traits.List(File(), maxlen=2, desc="Cortical thickness maps")
+    sulc = traits.List(File(), maxlen=2, desc="Sulcal depth maps")
+    curv = traits.List(File(), maxlen=2, desc="Curvature maps")
+
+
+class AggregateSurfaces(SimpleInterface):
+    """Aggregate and group surfaces & morphometrics into left/right pairs."""
+    input_spec = AggregateSurfacesInputSpec
+    output_spec = AggregateSurfacesOutputSpec
+
+    def _run_interface(self, runtime):
+        from collections import defaultdict
+        import os
+        import re
+
+        container = defaultdict(list)
+        inputs = (self.inputs.surfaces or []) + (self.inputs.morphometrics or [])
+        findre = re.compile(
+            r'(?:^|[^d])(?P<name>white|pial|inflated|midthickness|thickness|sulc|curv)'
+        )
+        for surface in sorted(inputs, key=os.path.basename):
+            match = findre.search(os.path.basename(surface))
+            if match:
+                container[match.group('name')].append(surface)
+        for name, files in container.items():
+            self._results[name] = files
+        return runtime
+
+
 def normalize_surfs(in_file: str, transform_file: str, newpath: Optional[str] = None) -> str:
     """
     Update GIFTI metadata and apply rigid coordinate correction.
diff --git a/smriprep/workflows/surfaces.py b/smriprep/workflows/surfaces.py
index b074c3d741..71b4dbecad 100644
--- a/smriprep/workflows/surfaces.py
+++ b/smriprep/workflows/surfaces.py
@@ -209,6 +209,13 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"):
                 "out_aseg",
                 "out_aparc",
                 "morphometrics",
+                "midthickness",
+                "pial",
+                "white",
+                "inflated",
+                "thickness",
+                "sulc",
+                "curv",
             ]
         ),
         name="outputnode",
@@ -294,7 +301,14 @@ def init_surface_recon_wf(*, omp_nthreads, hires, name="surface_recon_wf"):
         (autorecon_resume_wf, outputnode, [('outputnode.subjects_dir', 'subjects_dir'),
                                            ('outputnode.subject_id', 'subject_id')]),
         (gifti_surface_wf, outputnode, [('outputnode.surfaces', 'surfaces'),
-                                        ('outputnode.morphometrics', 'morphometrics')]),
+                                        ('outputnode.morphometrics', 'morphometrics'),
+                                        ('outputnode.midthickness', 'midthickness'),
+                                        ('outputnode.pial', 'pial'),
+                                        ('outputnode.white', 'white'),
+                                        ('outputnode.inflated', 'inflated'),
+                                        ('outputnode.thickness', 'thickness'),
+                                        ('outputnode.sulc', 'sulc'),
+                                        ('outputnode.curv', 'curv')]),
         (t1w2fsnative_xfm, outputnode, [('out_lta', 't1w2fsnative_xfm')]),
         (fsnative2t1w_xfm, outputnode, [('out_reg_file', 'fsnative2t1w_xfm')]),
         (refine, outputnode, [('out_file', 'out_brainmask')]),
@@ -574,15 +588,29 @@ def init_gifti_surface_wf(*, name="gifti_surface_wf"):
 
     Outputs
     -------
+    midthickness
+        Left and right midthickness (or graymid) surface GIFTIs
+    pial
+        Left and right pial surface GIFTIs
+    white
+        Left and right white surface GIFTIs
+    inflated
+        Left and right inflated surface GIFTIs
     surfaces
         GIFTI surfaces for gray/white matter boundary, pial surface,
         midthickness (or graymid) surface, and inflated surfaces
+    thickness
+        Left and right cortical thickness GIFTIs
+    sulc
+        Left and right sulcal depth map GIFTIs
+    curv
+        Left and right curvature map GIFTIs
     morphometrics
         GIFTIs of cortical thickness, curvature, and sulcal depth
 
     """
     from ..interfaces.freesurfer import MRIsConvertData
-    from ..interfaces.surf import NormalizeSurf
+    from ..interfaces.surf import NormalizeSurf, AggregateSurfaces
 
     workflow = Workflow(name=name)
 
@@ -590,7 +618,21 @@ def init_gifti_surface_wf(*, name="gifti_surface_wf"):
         niu.IdentityInterface(["subjects_dir", "subject_id", "fsnative2t1w_xfm"]),
         name="inputnode",
     )
-    outputnode = pe.Node(niu.IdentityInterface(["surfaces", "morphometrics"]), name="outputnode")
+    outputnode = pe.Node(
+        niu.IdentityInterface([
+            "pial",
+            "white",
+            "inflated",
+            "midthickness",
+            "thickness",
+            "sulc",
+            "curv",
+            # Preserve grouping
+            "surfaces",
+            "morphometrics",
+        ]),
+        name="outputnode"
+    )
 
     get_surfaces = pe.Node(nio.FreeSurferSource(), name="get_surfaces")
 
@@ -624,6 +666,8 @@ def init_gifti_surface_wf(*, name="gifti_surface_wf"):
         name="morphs2gii",
     )
 
+    agg_surfaces = pe.Node(AggregateSurfaces(), name="agg_surfaces")
+
     # fmt:off
     workflow.connect([
         (inputnode, get_surfaces, [('subjects_dir', 'subjects_dir'),
@@ -648,6 +692,16 @@ def init_gifti_surface_wf(*, name="gifti_surface_wf"):
                                         ('curv', 'in3')]),
         (surfmorph_list, morphs2gii, [('out', 'scalarcurv_file')]),
         (morphs2gii, outputnode, [('converted', 'morphometrics')]),
+        # Output individual surfaces as well
+        (fix_surfs, agg_surfaces, [('out_file', 'surfaces')]),
+        (morphs2gii, agg_surfaces, [('converted', 'morphometrics')]),
+        (agg_surfaces, outputnode, [('pial', 'pial'),
+                                    ('white', 'white'),
+                                    ('inflated', 'inflated'),
+                                    ('midthickness', 'midthickness'),
+                                    ('thickness', 'thickness'),
+                                    ('sulc', 'sulc'),
+                                    ('curv', 'curv')]),
     ])
     # fmt:on
     return workflow