diff --git a/fmriprep/data/boilerplate.bib b/fmriprep/data/boilerplate.bib
index 4cc802f25..66b9c540d 100644
--- a/fmriprep/data/boilerplate.bib
+++ b/fmriprep/data/boilerplate.bib
@@ -365,3 +365,17 @@ @article{patriat_improved_2017
     keywords = {Motion, Correction, Methods, Rs-fMRI},
     pages = {74--82},
 }
+
+@article{onavg,
+	author = {Feilong, Ma and Jiahui, Guo and Gobbini, Maria Ida and Haxby, James V.},
+	title = {A cortical surface template for human neuroscience},
+	url = {https://www.nature.com/articles/s41592-024-02346-y},
+	journal = {Nature Methods},
+	issn = {1548-7105},
+	number = {9},
+	volume = {21},
+	year = {2024},
+	month = sep,
+	pages = {1736--1742},
+	doi = {10.1038/s41592-024-02346-y},
+}
diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index bcd282d50..0cf2263b6 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/fmriprep/workflows/bold/base.py
@@ -499,6 +499,115 @@ def init_bold_wf(
             ]),
         ])  # fmt:skip
 
+    surf_std = spaces.get_standard(dim=(2,))
+    if surf_std:  # Probably ensure reconall and msmsulc are run
+        workflow.__postdesc__ += """\
+Non-gridded (surface) resamplings were performed using the Connectome
+Workbench.
+"""
+        config.loggers.workflow.debug('Creating BOLD surface workbench resampling workflow.')
+        from smriprep.workflows.surfaces import init_resample_surfaces_wb_wf
+
+        from .resampling import (
+            init_goodvoxels_bold_mask_wf,
+            init_wb_surf_surf_wf,
+            init_wb_vol_surf_wf,
+        )
+
+        wb_vol_surf_wf = init_wb_vol_surf_wf(
+            name='wb_vol_surf_wf',
+            omp_nthreads=omp_nthreads,
+            mem_gb=mem_gb['resampled'],
+            dilate=True,
+        )
+        workflow.connect([
+            (inputnode, wb_vol_surf_wf,[
+                ('white', 'inputnode.white'),
+                ('pial', 'inputnode.pial'),
+                ('midthickness', 'inputnode.midthickness'),
+            ]),
+            (bold_anat_wf, wb_vol_surf_wf, [
+                ('outputnode.bold_file', 'inputnode.bold_file'),
+            ]),
+        ])  # fmt:skip
+
+        if config.workflow.project_goodvoxels:
+            goodvoxels_bold_mask_wf = init_goodvoxels_bold_mask_wf(mem_gb['resampled'])
+
+            workflow.connect([
+                (inputnode, goodvoxels_bold_mask_wf, [('anat_ribbon', 'inputnode.anat_ribbon')]),
+                (bold_anat_wf, goodvoxels_bold_mask_wf, [
+                    ('outputnode.bold_file', 'inputnode.bold_file'),
+                ]),
+                (goodvoxels_bold_mask_wf, wb_vol_surf_wf, [
+                    ('outputnode.goodvoxels_mask', 'inputnode.volume_roi'),
+                ]),
+            ])  # fmt:skip
+            workflow.__desc__ += """\
+A "goodvoxels" mask was applied during volume-to-surface sampling, excluding
+voxels whose time-series have a locally high coefficient of variation.
+"""
+
+        for ref_ in surf_std:
+            space, den = ref_.space, ref_.spec['den']
+
+            resample_surfaces_wb_wf = init_resample_surfaces_wb_wf(
+                name=f'resample_surfaces_wb_wf_{space}_{den}',
+                surfaces=['midthickness'],
+                space=space,
+                density=den,
+            )
+
+            wb_surf_surf_wf = init_wb_surf_surf_wf(
+                space=space,
+                density=den,
+                name=f'wb_surf_surf_wf_{space}_{den}',
+                omp_nthreads=omp_nthreads,
+                mem_gb=mem_gb['resampled'],
+            )
+
+            ds_bold_surf_wb = pe.Node(
+                DerivativesDataSink(
+                    base_directory=fmriprep_dir,
+                    hemi=['L', 'R'],
+                    dismiss_entities=dismiss_echo(),
+                    space=space,
+                    density=den,
+                    suffix='bold',
+                    # compress=False, # not sure if needed for gii.
+                    TaskName=all_metadata[0].get('TaskName'),
+                    extension='.func.gii',
+                    **prepare_timing_parameters(all_metadata[0]),
+                ),
+                iterfield=('in_file', 'hemi'),
+                name=f'ds_bold_surf_wb_{space}_{den}',
+                run_without_submitting=True,
+            )
+            ds_bold_surf_wb.inputs.source_file = bold_file
+
+            workflow.connect([
+                (inputnode, resample_surfaces_wb_wf, [
+                    ('midthickness', 'inputnode.midthickness'),
+                    ('sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR'),
+                ]),
+                (wb_vol_surf_wf, wb_surf_surf_wf, [
+                    ('outputnode.bold_fsnative', 'inputnode.bold_fsnative'),
+                ]),
+                (inputnode, wb_surf_surf_wf, [
+                    ('midthickness', 'inputnode.midthickness'),
+                    # # TODO: check inputnode.midthickness_resampled
+                    # ('midthickness_resampled', 'inputnode.midthickness_resampled'),
+                    ('sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR'),
+                ]),
+                (resample_surfaces_wb_wf, wb_surf_surf_wf, [
+                    ('outputnode.midthickness_resampled', 'inputnode.midthickness_resampled'),
+                ]),
+                (wb_surf_surf_wf, ds_bold_surf_wb, [
+                    ('outputnode.bold_resampled', 'in_file'),
+                    # TODO: json metadata?
+                ]),
+            ])  # fmt:skip
+
     if config.workflow.run_reconall and freesurfer_spaces:
         workflow.__postdesc__ += """\
 Non-gridded (surface) resamplings were performed using `mri_vol2surf`
diff --git a/fmriprep/workflows/bold/resampling.py b/fmriprep/workflows/bold/resampling.py
index d725da6f0..0d6bce20e 100644
--- a/fmriprep/workflows/bold/resampling.py
+++ b/fmriprep/workflows/bold/resampling.py
@@ -25,6 +25,8 @@
 ++++++++++++++++++++
 
 .. autofunction:: init_bold_surf_wf
+.. autofunction:: init_wb_vol_surf_wf
+.. autofunction:: init_wb_surf_surf_wf
 .. autofunction:: init_bold_fsLR_resampling_wf
 .. autofunction:: init_bold_grayords_wf
 .. autofunction:: init_goodvoxels_bold_mask_wf
@@ -514,6 +516,325 @@ def _calc_lower_thr(in_stats):
     return workflow
 
 
+def init_wb_vol_surf_wf(
+    omp_nthreads: int,
+    mem_gb: float,
+    name: str = 'wb_vol_surf_wf',
+    dilate: bool = True,
+):
+    """Resample volume to native surface and dilate it using the Workbench.
+
+    This workflow performs the first two steps of surface resampling:
+    1. Resample volume to native surface using "ribbon-constrained" method
+    2. Dilate the resampled surface to fix small holes using nearest neighbors
+
+    The output of this workflow can be reused to resample to multiple template
+    spaces and resolutions.
+
+    Workflow Graph
+        .. workflow::
+            :graph2use: colored
+            :simple_form: yes
+
+            from fmriprep.workflows.bold.resampling import init_wb_vol_surf_wf
+            wf = init_wb_vol_surf_wf(omp_nthreads=1, mem_gb=1)
+
+
+    Parameters
+    ----------
+    omp_nthreads : :class:`int`
+        Maximum number of threads an individual process may use.
+    mem_gb : :class:`float`
+        Size of BOLD file in GB.
+    name : :class:`str`
+        Name of workflow (default: ``wb_vol_surf_wf``).
+
+    Inputs
+    ------
+    bold_file : :class:`str`
+        Path to BOLD file resampled into T1 space
+    white : :class:`list` of :class:`str`
+        Path to left and right hemisphere white matter GIFTI surfaces.
+    pial : :class:`list` of :class:`str`
+        Path to left and right hemisphere pial GIFTI surfaces.
+    midthickness : :class:`list` of :class:`str`
+        Path to left and right hemisphere midthickness GIFTI surfaces.
+    volume_roi : :class:`str` or Undefined
+        Pre-calculated goodvoxels mask. Not required.
+
+    Outputs
+    -------
+    bold_fsnative : :class:`list` of :class:`str`
+        Path to BOLD series resampled as functional GIFTI files in native
+        surface space.
+    """
+    from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+    from niworkflows.interfaces.utility import KeySelect
+
+    from fmriprep.interfaces.workbench import VolumeToSurfaceMapping
+
+    workflow = Workflow(name=name)
+    workflow.__desc__ = """\
+The BOLD time-series were resampled onto the native surface of the subject
+using the "ribbon-constrained" method
+"""
+    workflow.__desc__ += ' and then dilated by 10 mm.' if dilate else '.'
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                'bold_file',
+                'white',
+                'pial',
+                'midthickness',
+                'volume_roi',
+            ]
+        ),
+        name='inputnode',
+    )
+
+    hemisource = pe.Node(
+        niu.IdentityInterface(fields=['hemi']),
+        name='hemisource_vol_surf',
+        iterables=[('hemi', ['L', 'R'])],
+    )
+
+    joinnode = pe.JoinNode(
+        niu.IdentityInterface(fields=['bold_fsnative']),
+        name='joinnode_vol_surf',
+        joinsource='hemisource_vol_surf',
+    )
+
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=['bold_fsnative']),
+        name='outputnode',
+    )
+
+    select_surfaces = pe.Node(
+        KeySelect(
+            fields=[
+                'white',
+                'pial',
+                'midthickness',
+            ],
+            keys=['L', 'R'],
+        ),
+        name='select_surfaces',
+        run_without_submitting=True,
+    )
+
+    volume_to_surface = pe.Node(
+        VolumeToSurfaceMapping(method='ribbon-constrained'),
+        name='volume_to_surface',
+        mem_gb=mem_gb * 3,
+        n_procs=omp_nthreads,
+    )
+    if dilate:
+        metric_dilate = pe.Node(
+            MetricDilate(distance=10, nearest=True),
+            name='metric_dilate',
+            mem_gb=1,
+            n_procs=omp_nthreads,
+        )
+
+    workflow.connect([
+        (inputnode, select_surfaces, [
+            ('white', 'white'),
+            ('pial', 'pial'),
+            ('midthickness', 'midthickness'),
+        ]),
+        (hemisource, select_surfaces, [('hemi', 'key')]),
+        (inputnode, volume_to_surface, [
+            ('bold_file', 'volume_file'),
+            ('volume_roi', 'volume_roi'),
+        ]),
+        (select_surfaces, volume_to_surface, [
+            ('midthickness', 'surface_file'),
+            ('white', 'inner_surface'),
+            ('pial', 'outer_surface'),
+        ]),
+    ])  # fmt:skip
+    if dilate:
+        workflow.connect([
+            (select_surfaces, metric_dilate, [
+                ('midthickness', 'surf_file'),
+            ]),
+            (volume_to_surface, metric_dilate, [
+                ('out_file', 'in_file'),
+            ]),
+            (metric_dilate, joinnode, [
+                ('out_file', 'bold_fsnative'),
+            ]),
+            (joinnode, outputnode, [
+                ('bold_fsnative', 'bold_fsnative'),
+            ]),
+        ])  # fmt:skip
+    else:
+        workflow.connect([
+            (volume_to_surface, joinnode, [
+                ('out_file', 'bold_fsnative'),
+            ]),
+            (joinnode, outputnode, [
+                ('bold_fsnative', 'bold_fsnative'),
+            ]),
+        ])  # fmt:skip
+
+    return workflow
+
+
+def init_wb_surf_surf_wf(
+    space: str,
+    density: ty.Literal['10k', '32k', '41k'],
+    omp_nthreads: int,
+    mem_gb: float,
+    name: str = 'wb_surf_surf_wf',
+):
+    """Resample BOLD time series from native surface to template surface.
+
+    This workflow performs the third step of surface resampling:
+    3. Resample the native surface to the template surface using the
+       Connectome Workbench
+
+    Workflow Graph
+        .. workflow::
+            :graph2use: colored
+            :simple_form: yes
+
+            from fmriprep.workflows.bold.resampling import init_wb_surf_surf_wf
+            wf = init_wb_surf_surf_wf(
+                space='fsLR',
+                density='32k',
+                omp_nthreads=1,
+                mem_gb=1,
+            )
+
+    Parameters
+    ----------
+    space : :class:`str`
+        Surface template space, such as ``"onavg"`` or ``"fsLR"``.
+    density : :class:`str`
+        Either ``"10k"``, ``"32k"``, or ``"41k"``, representing the number of
+        vertices per hemisphere.
+    omp_nthreads : :class:`int`
+        Maximum number of threads an individual process may use.
+    mem_gb : :class:`float`
+        Size of BOLD file in GB.
+    name : :class:`str`
+        Name of workflow (default: ``wb_surf_surf_wf``).
+
+    Inputs
+    ------
+    bold_fsnative : :class:`list` of :class:`str`
+        Path to BOLD series resampled as functional GIFTI files in native
+        surface space.
+    midthickness : :class:`list` of :class:`str`
+        Path to left and right hemisphere midthickness GIFTI surfaces.
+    midthickness_resampled : :class:`list` of :class:`str`
+        Path to left and right hemisphere midthickness GIFTI surfaces resampled
+        into the output space.
+    sphere_reg_fsLR : :class:`list` of :class:`str`
+        Path to left and right hemisphere sphere.reg GIFTI surfaces, mapping
+        from subject to fsLR.
+
+    Outputs
+    -------
+    bold_resampled : :class:`list` of :class:`str`
+        Path to BOLD series resampled as functional GIFTI files in the output
+        template space.
+    """
+    import templateflow.api as tf
+    from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+    from niworkflows.interfaces.utility import KeySelect
+
+    workflow = Workflow(name=name)
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=[
+                'bold_fsnative',
+                'midthickness',
+                'midthickness_resampled',
+                'sphere_reg_fsLR',
+            ]
+        ),
+        name='inputnode',
+    )
+
+    hemisource = pe.Node(
+        niu.IdentityInterface(fields=['hemi']),
+        name=name + '_hemisource_surf_surf',
+        iterables=[('hemi', ['L', 'R'])],
+    )
+
+    joinnode = pe.JoinNode(
+        niu.IdentityInterface(fields=['bold_resampled']),
+        name=name + '_joinnode_surf_surf',
+        joinsource=name + '_hemisource_surf_surf',
+    )
+
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=['bold_resampled']),
+        name='outputnode',
+    )
+
+    select_surfaces = pe.Node(
+        KeySelect(
+            fields=[
+                'bold_fsnative',
+                'midthickness',
+                'midthickness_resampled',
+                'sphere_reg_fsLR',
+                'template_sphere',
+            ],
+            keys=['L', 'R'],
+        ),
+        name=name + '_select_surfaces',
+        run_without_submitting=True,
+    )
+    select_surfaces.inputs.template_sphere = [
+        str(sphere)
+        for sphere in tf.get(
+            template=space,
+            space=('fsLR' if space != 'fsLR' else None),
+            density=density,
+            suffix='sphere',
+            extension='.surf.gii',
+        )
+    ]
+
+    resample_to_template = pe.Node(
+        MetricResample(method='ADAP_BARY_AREA', area_surfs=True),
+        name=name + '_resample_to_template',
+        mem_gb=1,
+        n_procs=omp_nthreads,
+    )
+
+    workflow.connect([
+        (inputnode, select_surfaces, [
+            ('bold_fsnative', 'bold_fsnative'),
+            ('midthickness', 'midthickness'),
+            ('midthickness_resampled', 'midthickness_resampled'),
+            ('sphere_reg_fsLR', 'sphere_reg_fsLR'),
+        ]),
+        (hemisource, select_surfaces, [('hemi', 'key')]),
+        (select_surfaces, resample_to_template, [
+            ('bold_fsnative', 'in_file'),
+            ('sphere_reg_fsLR', 'current_sphere'),
+            ('template_sphere', 'new_sphere'),
+            ('midthickness', 'current_area'),
+            ('midthickness_resampled', 'new_area'),
+        ]),
+        (resample_to_template, joinnode, [
+            ('out_file', 'bold_resampled'),
+        ]),
+        (joinnode, outputnode, [
+            ('bold_resampled', 'bold_resampled'),
+        ]),
+    ])  # fmt:skip
+
+    return workflow
+
+
 def init_bold_fsLR_resampling_wf(
     grayord_density: ty.Literal['91k', '170k'],
     omp_nthreads: int,