diff --git a/docs/faq.rst b/docs/faq.rst
index 5c54a360c..8e62be29f 100644
--- a/docs/faq.rst
+++ b/docs/faq.rst
@@ -290,8 +290,8 @@ Some examples follow:
 
 * Surgery: use only pre-operation sessions for the anatomical data. This will typically be done
   by omitting post-operation sessions from the inputs to *fMRIPrep*.
-* Developing and elderly populations: there is currently no standard way of processing these.
-  However, `as suggested by U. Tooley at NeuroStars.org
+* Developing and elderly populations: ``--subject-anatomical-reference sessionwise`` will process each session independently.
+  Additionally, `as suggested by U. Tooley at NeuroStars.org
   <https://neurostars.org/t/fmriprep-how-to-reuse-longitudinal-and-pre-run-freesurfer/4585/15>`__,
   it is theoretically possible to leverage the *anatomical fast-track* along with the
   ``--bids-filters`` option to process sessions fully independently, or grouped by some study-design
diff --git a/docs/workflows.rst b/docs/workflows.rst
index 135a5ab76..482763dfd 100644
--- a/docs/workflows.rst
+++ b/docs/workflows.rst
@@ -174,15 +174,22 @@ Your ``.bidsignore`` file should include the following line::
 
 Longitudinal processing
 ~~~~~~~~~~~~~~~~~~~~~~~
-In the case of multiple T1w images (across sessions and/or runs), T1w images are
-merged into a single template image using FreeSurfer's `mri_robust_template`_.
-This template may be *unbiased*, or equidistant from all source images, or
-aligned to the first image (determined lexicographically by session label).
+In the case of multiple T1w images (across sessions and/or runs), *fMRIPrep* provides
+a few choices on how to generate the reference anatomical space.
+
+If ``--subject-anatomical-reference first-lex`` is used, all T1w images are
+merged into a single template image using FreeSurfer's `mri_robust_template`_,
+aligned to the first image (determined lexicographically by session label). This is
+the default behavior.
+
+If ``--subject-anatomical-reference unbiased`` is used, all T1w images are merged into
+an *unbiased* template, equidistant from all source images.
 For two images, the additional cost of estimating an unbiased template is trivial,
 but aligning three or more images is too expensive to justify being the default behavior.
-For consistency, in the case of multiple images, *fMRIPrep* constructs
-templates aligned to the first image, unless passed the ``--longitudinal``
-flag, which forces the estimation of an unbiased template.
+
+If ``--subject-anatomical-reference sessionwise`` is used, a reference template will be
+generated for each session independently. If multiple T1w images are found within a session,
+the images will be aligned to the first image, sorted lexicographically, from that session.
 
 .. note::
 
diff --git a/fmriprep/cli/parser.py b/fmriprep/cli/parser.py
index a4cf741a7..9f437ada3 100644
--- a/fmriprep/cli/parser.py
+++ b/fmriprep/cli/parser.py
@@ -23,10 +23,14 @@
 """Parser."""
 
 import sys
+import typing as ty
 from pathlib import Path
 
 from .. import config
 
+if ty.TYPE_CHECKING:
+    from bids import BIDSLayout
+
 
 def _build_parser(**kwargs):
     """Build parser object.
@@ -51,6 +55,7 @@ def _build_parser(**kwargs):
         'force_bbr': ('--force bbr', '26.0.0'),
         'force_no_bbr': ('--force no-bbr', '26.0.0'),
         'force_syn': ('--force syn-sdc', '26.0.0'),
+        'longitudinal': ('--subject-anatomical-reference unbiased', '26.1.0'),
     }
 
     class DeprecatedAction(Action):
@@ -220,9 +225,14 @@ def _fallback_trt(value, parser):
         help='A space delimited list of participant identifiers or a single '
         'identifier (the sub- prefix can be removed)',
     )
-    # Re-enable when option is actually implemented
-    # g_bids.add_argument('-s', '--session-id', action='store', default='single_session',
-    #                     help='Select a specific session to be processed')
+
+    g_bids.add_argument(
+        '--session-label',
+        nargs='+',
+        type=lambda label: label.removeprefix('ses-'),
+        help='A space delimited list of session identifiers or a single '
+        'identifier (the ses- prefix can be removed)',
+    )
     # Re-enable when option is actually implemented
     # g_bids.add_argument('-r', '--run-id', action='store', default='single_run',
     #                     help='Select a specific run to be processed')
@@ -235,6 +245,17 @@ def _fallback_trt(value, parser):
         type=int,
         help='Select a specific echo to be processed in a multiecho series',
     )
+    g_bids.add_argument(
+        '--subject-anatomical-reference',
+        choices=['first-lex', 'unbiased', 'sessionwise'],
+        default='first-lex',
+        help='Method to produce the reference anatomical space:\n'
+        '\t"first-lex" will use the first image in lexicographical order\n'
+        '\t"unbiased" will construct an unbiased template from all images '
+        '(previously "--longitudinal")\n'
+        '\t"sessionwise" will independently process each session. If multiple runs are '
+        'found, the behavior will be similar to "first-lex" for each session.',
+    )
     g_bids.add_argument(
         '--bids-filter-file',
         dest='bids_filters',
@@ -384,8 +405,8 @@ def _fallback_trt(value, parser):
     )
     g_conf.add_argument(
         '--longitudinal',
-        action='store_true',
-        help='Treat dataset as longitudinal - may increase runtime',
+        action=DeprecatedAction,
+        help='Deprecated - use `--subject-anatomical-reference unbiased` instead',
     )
     g_conf.add_argument(
         '--bold2anat-init',
@@ -769,6 +790,7 @@ def parse_args(args=None, namespace=None):
     """Parse args and run further checks on the command line."""
     import logging
 
+    from niworkflows.utils.bids import collect_participants
     from niworkflows.utils.spaces import Reference, SpatialReferences
 
     parser = _build_parser()
@@ -779,6 +801,14 @@ def parse_args(args=None, namespace=None):
         config.load(opts.config_file, skip=skip, init=False)
         config.loggers.cli.info(f'Loaded previous configuration file {opts.config_file}')
 
+    if opts.longitudinal:
+        opts.subject_anatomical_reference = 'unbiased'
+        msg = (
+            'The `--longitudinal` flag is deprecated - use '
+            '`--subject-anatomical-reference unbiased` instead.'
+        )
+        config.loggers.cli.warning(msg)
+
     config.execution.log_level = int(max(25 - 5 * opts.verbose_count, logging.DEBUG))
     config.from_dict(vars(opts), init=['nipype'])
 
@@ -927,19 +957,67 @@ def parse_args(args=None, namespace=None):
     work_dir.mkdir(exist_ok=True, parents=True)
 
     # Force initialization of the BIDSLayout
+    config.loggers.cli.debug('Initializing BIDS Layout')
     config.execution.init()
-    all_subjects = config.execution.layout.get_subjects()
+
+    # Please note this is the input folder's dataset_description.json
+    dset_desc_path = config.execution.bids_dir / 'dataset_description.json'
+    if dset_desc_path.exists():
+        from hashlib import sha256
+
+        desc_content = dset_desc_path.read_bytes()
+        config.execution.bids_description_hash = sha256(desc_content).hexdigest()
+
+    # First check that bids_dir looks like a BIDS folder
+    subject_list = collect_participants(
+        config.execution.layout, participant_label=config.execution.participant_label
+    )
     if config.execution.participant_label is None:
-        config.execution.participant_label = all_subjects
+        config.execution.participant_label = subject_list
 
-    participant_label = set(config.execution.participant_label)
-    missing_subjects = participant_label - set(all_subjects)
-    if missing_subjects:
-        parser.error(
-            'One or more participant labels were not found in the BIDS directory: {}.'.format(
-                ', '.join(missing_subjects)
+    session_list = config.execution.session_label or []
+    subject_session_list = create_processing_groups(
+        config.execution.layout,
+        subject_list,
+        session_list,
+        config.workflow.subject_anatomical_reference,
+    )
+    config.execution.processing_groups = subject_session_list
+    config.execution.participant_label = sorted(subject_list)
+    config.workflow.skull_strip_template = config.workflow.skull_strip_template[0]
+
+
+def create_processing_groups(
+    layout: 'BIDSLayout',
+    subject_list: list,
+    session_list: list | str | None,
+    subject_anatomical_reference: str,
+) -> list[tuple[str]]:
+    """Generate a list of subject-session pairs to be processed."""
+    from bids.layout import Query
+
+    subject_session_list = []
+
+    for subject in subject_list:
+        sessions = (
+            layout.get_sessions(
+                scope='raw',
+                subject=subject,
+                session=session_list or Query.OPTIONAL,
             )
+            or None
         )
 
-    config.execution.participant_label = sorted(participant_label)
-    config.workflow.skull_strip_template = config.workflow.skull_strip_template[0]
+        if subject_anatomical_reference == 'sessionwise':
+            if sessions is None:
+                config.loggers.cli.warning(
+                    '`--subject-anatomical-reference sessionwise` was requested, but no sessions '
+                    f'found for subject {subject}... treating as single-session.'
+                )
+                subject_session_list.append((subject, None))
+            else:
+                subject_session_list.extend((subject, session) for session in sessions)
+        else:
+            subject_session_list.append((subject, sessions))
+
+    return subject_session_list
diff --git a/fmriprep/cli/workflow.py b/fmriprep/cli/workflow.py
index 29f2bfdaf..2f3d7b518 100644
--- a/fmriprep/cli/workflow.py
+++ b/fmriprep/cli/workflow.py
@@ -35,15 +35,13 @@
 def build_workflow(config_file, retval):
     """Create the Nipype Workflow that supports the whole execution graph."""
 
-    from niworkflows.utils.bids import collect_participants
     from niworkflows.utils.misc import check_valid_fs_license
 
+    from fmriprep import config, data
     from fmriprep.reports.core import generate_reports
     from fmriprep.utils.bids import check_pipeline_version
-
-    from .. import config, data
-    from ..utils.misc import check_deps
-    from ..workflows.base import init_fmriprep_wf
+    from fmriprep.utils.misc import check_deps, fmt_subjects_sessions
+    from fmriprep.workflows.base import init_fmriprep_wf
 
     config.load(config_file)
     build_log = config.loggers.workflow
@@ -69,27 +67,20 @@ def build_workflow(config_file, retval):
     if msg is not None:
         build_log.warning(msg)
 
-    # Please note this is the input folder's dataset_description.json
-    dset_desc_path = config.execution.bids_dir / 'dataset_description.json'
-    if dset_desc_path.exists():
-        from hashlib import sha256
-
-        desc_content = dset_desc_path.read_bytes()
-        config.execution.bids_description_hash = sha256(desc_content).hexdigest()
-
-    # First check that bids_dir looks like a BIDS folder
-    subject_list = collect_participants(
-        config.execution.layout, participant_label=config.execution.participant_label
-    )
-
     # Called with reports only
     if config.execution.reports_only:
-        build_log.log(25, 'Running --reports-only on participants %s', ', '.join(subject_list))
-        session_list = (
-            config.execution.bids_filters.get('bold', {}).get('session')
-            if config.execution.bids_filters
-            else None
+        build_log.log(
+            25,
+            'Running --reports-only on %s',
+            fmt_subjects_sessions(config.execution.processing_groups),
         )
+        session_list = config.execution.session_label
+        if not session_list:
+            session_list = (
+                config.execution.bids_filters.get('bold', {}).get('session')
+                if config.execution.bids_filters
+                else None
+            )
 
         failed_reports = generate_reports(
             config.execution.participant_label,
@@ -110,7 +101,7 @@ def build_workflow(config_file, retval):
     init_msg = [
         "Building fMRIPrep's workflow:",
         f'BIDS dataset path: {config.execution.bids_dir}.',
-        f'Participant list: {subject_list}.',
+        f'Participants and sessions: {fmt_subjects_sessions(config.execution.processing_groups)}.',
         f'Run identifier: {config.execution.run_uuid}.',
         f'Output spaces: {config.execution.output_spaces}.',
     ]
diff --git a/fmriprep/config.py b/fmriprep/config.py
index c9f6efd49..11ca372a3 100644
--- a/fmriprep/config.py
+++ b/fmriprep/config.py
@@ -232,6 +232,12 @@ def load(cls, settings, init=True, ignore=None):
                 else:
                     setattr(cls, k, Path(v).absolute())
             elif hasattr(cls, k):
+                match k:
+                    # Handle special deserializations
+                    case 'processing_groups':
+                        v = _deserialize_pg(v)
+                    case _:
+                        pass
                 setattr(cls, k, v)
 
         if init:
@@ -434,8 +440,12 @@ class execution(_Config):
     """Only build the reports, based on the reportlets found in a cached working directory."""
     run_uuid = f'{strftime("%Y%m%d-%H%M%S")}_{uuid4()}'
     """Unique identifier of this particular run."""
+    processing_groups = None
+    """List of tuples (participant, session(s)) that will be preprocessed."""
     participant_label = None
     """List of participant identifiers that are to be preprocessed."""
+    session_label = None
+    """List of session identifiers that are to be preprocessed."""
     task_id = None
     """Select a particular task from all available in the dataset."""
     templateflow_home = _templateflow_home
@@ -588,8 +598,6 @@ class workflow(_Config):
     """Force particular steps for *fMRIPrep*."""
     level = None
     """Level of preprocessing to complete. One of ['minimal', 'resampling', 'full']."""
-    longitudinal = False
-    """Run FreeSurfer ``recon-all`` with the ``-logitudinal`` flag."""
     run_msmsulc = True
     """Run Multimodal Surface Matching surface registration."""
     medial_surface_nan = None
@@ -619,6 +627,11 @@ class workflow(_Config):
     spaces = None
     """Keeps the :py:class:`~niworkflows.utils.spaces.SpatialReferences`
     instance keeping standard and nonstandard spaces."""
+    subject_anatomical_reference = 'first-lex'
+    """Method to produce the reference anatomical space. Available options are:
+    `first-lex` will use the first image in lexicographical order, `unbiased` will
+    construct an unbiased template from all available images (previously --longitudinal),
+    and `sessionwise` will independently process each session."""
     use_aroma = None
     """Run ICA-:abbr:`AROMA (automatic removal of motion artifacts)`."""
     use_bbr = None
@@ -776,6 +789,7 @@ def get(flat=False):
         'nipype': nipype.get(),
         'seeds': seeds.get(),
     }
+
     if not flat:
         return settings
 
@@ -790,7 +804,12 @@ def dumps():
     """Format config into toml."""
     from toml import dumps
 
-    return dumps(get())
+    settings = get()
+    # Serialize to play nice with TOML
+    if pg := settings['execution'].get('processing_groups'):
+        settings['execution']['processing_groups'] = _serialize_pg(pg)
+
+    return dumps(settings)
 
 
 def to_filename(filename):
@@ -827,3 +846,53 @@ def init_spaces(checkpoint=True):
 
     # Make the SpatialReferences object available
     workflow.spaces = spaces
+
+
+def _serialize_pg(value: list[tuple[str, list[str] | None]]) -> list[str]:
+    """
+    Serialize a list of participant-session tuples to be TOML-compatible.
+
+    Examples
+    --------
+    >>> _serialize_pg([('01', ['pre']), ('01', ['post'])])
+    ['sub-01_ses-pre', 'sub-01_ses-post']
+    >>> _serialize_pg([('01', ['pre', 'post']), ('02', ['post'])])
+    ['sub-01_ses-pre,post', 'sub-02_ses-post']
+    >>> _serialize_pg([('01', None), ('02', ['pre'])])
+    ['sub-01', 'sub-02_ses-pre']
+    """
+    serial = []
+    for val in value:
+        if val[1] is None:
+            serial.append(f'sub-{val[0]}')
+        else:
+            if not isinstance(val[1], list):
+                val[1] = [val[1]]
+            serial.append(f'sub-{val[0]}_ses-{",".join(val[1])}')
+    return serial
+
+
+def _deserialize_pg(value: list[str]) -> list[tuple[str, list[str] | None]]:
+    """
+    Deserialize a list of participant-session tuples to be TOML-compatible.
+
+    Examples
+    --------
+    >>> _deserialize_pg(['sub-01_ses-pre', 'sub-01_ses-post'])
+    [('01', ['pre']), ('01', ['post'])]
+    >>> _deserialize_pg(['sub-01_ses-pre,post', 'sub-02_ses-post'])
+    [('01', ['pre', 'post']), ('02', ['post'])]
+    >>> _deserialize_pg(['sub-01', 'sub-02_ses-pre'])
+    [('01', None), ('02', ['pre'])]
+    """
+    deserial = []
+    for val in value:
+        vals = val.split('_', 1)
+        vals[0] = vals[0].removeprefix('sub-')
+        if len(vals) == 1:
+            vals.append(None)
+        else:
+            vals[1] = vals[1].removeprefix('ses-')
+            vals[1] = vals[1].split(',')
+        deserial.append(tuple(vals))
+    return deserial
diff --git a/fmriprep/data/tests/config.toml b/fmriprep/data/tests/config.toml
index 3af10cc50..89974c7da 100644
--- a/fmriprep/data/tests/config.toml
+++ b/fmriprep/data/tests/config.toml
@@ -24,21 +24,19 @@ output_spaces = "MNI152NLin2009cAsym:res-2 MNI152NLin2009cAsym:res-native fsaver
 reports_only = false
 run_uuid = "20200306-105302_d365772b-fd60-4741-a722-372c2f558b50"
 participant_label = [ "01",]
+processing_groups = [ "sub-01",]
 templateflow_home = "~/.cache/templateflow"
 work_dir = "work/"
 write_graph = false
 
 [workflow]
 anat_only = false
-aroma_err_on_warn = false
-aroma_melodic_dim = -200
 bold2anat_dof = 6
 fmap_bspline = false
 force = []
 force_syn = false
 hires = true
 ignore = []
-longitudinal = false
 medial_surface_nan = false
 project_goodvoxels = false
 regressors_all_comps = false
@@ -47,8 +45,8 @@ regressors_fd_th = 0.5
 run_reconall = true
 skull_strip_fixed_seed = false
 skull_strip_template = "OASIS30ANTs"
+subject_anatomical_reference = "first-lex"
 t2s_coreg = false
-use_aroma = false
 
 [nipype]
 crashfile_format = "txt"
diff --git a/fmriprep/interfaces/bids.py b/fmriprep/interfaces/bids.py
index 87ec5d630..eff5a3a06 100644
--- a/fmriprep/interfaces/bids.py
+++ b/fmriprep/interfaces/bids.py
@@ -5,6 +5,7 @@
 from bids.utils import listify
 from nipype.interfaces.base import (
     DynamicTraitedSpec,
+    File,
     SimpleInterface,
     TraitedSpec,
     isdefined,
@@ -60,3 +61,130 @@ def _run_interface(self, runtime):
         self._results['out'] = out
 
         return runtime
+
+
+class _BIDSSourceFileInputSpec(TraitedSpec):
+    bids_info = traits.Dict(
+        mandatory=True,
+        desc='BIDS information dictionary',
+    )
+    precomputed = traits.Dict({}, usedefault=True, desc='Precomputed BIDS information')
+    sessionwise = traits.Bool(False, usedefault=True, desc='Keep session information')
+    anat_type = traits.Enum('t1w', 't2w', usedefault=True, desc='Anatomical reference type')
+
+
+class _BIDSSourceFileOutputSpec(TraitedSpec):
+    source_file = File(desc='Source file')
+
+
+class BIDSSourceFile(SimpleInterface):
+    input_spec = _BIDSSourceFileInputSpec
+    output_spec = _BIDSSourceFileOutputSpec
+
+    def _run_interface(self, runtime):
+        src = self.inputs.bids_info[self.inputs.anat_type]
+
+        if not src and self.inputs.precomputed.get(f'{self.inputs.anat_type}_preproc'):
+            src = self.inputs.bids_info['bold']
+            self._results['source_file'] = _create_multi_source_file(src)
+            return runtime
+
+        self._results['source_file'] = _create_multi_source_file(
+            src,
+            sessionwise=self.inputs.sessionwise,
+        )
+        return runtime
+
+
+class _CreateFreeSurferIDInputSpec(TraitedSpec):
+    subject_id = traits.Str(mandatory=True, desc='BIDS Subject ID')
+    session_id = traits.Str(desc='BIDS session ID')
+
+
+class _CreateFreeSurferIDOutputSpec(TraitedSpec):
+    subject_id = traits.Str(desc='FreeSurfer subject ID')
+
+
+class CreateFreeSurferID(SimpleInterface):
+    input_spec = _CreateFreeSurferIDInputSpec
+    output_spec = _CreateFreeSurferIDOutputSpec
+
+    def _run_interface(self, runtime):
+        self._results['subject_id'] = _create_fs_id(
+            self.inputs.subject_id,
+            self.inputs.session_id or None,
+        )
+        return runtime
+
+
+def _create_multi_source_file(in_files, sessionwise=False):
+    """
+    Create a generic source name from multiple input files.
+
+    If sessionwise is True, session information from the first file is retained in the name.
+
+    Examples
+    --------
+    >>> _create_multi_source_file([
+    ...     '/path/to/sub-045_ses-test_T1w.nii.gz',
+    ...     '/path/to/sub-045_ses-retest_T1w.nii.gz'])
+    '/path/to/sub-045_T1w.nii.gz'
+    >>> _create_multi_source_file([
+    ...     '/path/to/sub-045_ses-1_run-1_T1w.nii.gz',
+    ...     '/path/to/sub-045_ses-1_run-2_T1w.nii.gz'],
+    ...     sessionwise=True)
+    '/path/to/sub-045_ses-1_T1w.nii.gz'
+    """
+    import re
+    from pathlib import Path
+
+    from nipype.utils.filemanip import filename_to_list
+
+    if not isinstance(in_files, tuple | list):
+        return in_files
+    elif len(in_files) == 1:
+        return in_files[0]
+
+    p = Path(filename_to_list(in_files)[0])
+    try:
+        subj = re.search(r'(?<=^sub-)[a-zA-Z0-9]*', p.name).group()
+        suffix = re.search(r'(?<=_)\w+(?=\.)', p.name).group()
+    except AttributeError as e:
+        raise AttributeError('Could not extract BIDS information') from e
+
+    prefix = f'sub-{subj}'
+
+    if sessionwise:
+        ses = re.search(r'(?<=_ses-)[a-zA-Z0-9]*', p.name)
+        if ses:
+            prefix += f'_ses-{ses.group()}'
+    return str(p.parent / f'{prefix}_{suffix}.nii.gz')
+
+
+def _create_fs_id(subject_id, session_id=None):
+    """
+    Create FreeSurfer subject ID.
+
+    Examples
+    --------
+    >>> _create_fs_id('01')
+    'sub-01'
+    >>> _create_fs_id('sub-01')
+    'sub-01'
+    >>> _create_fs_id('01', 'pre')
+    'sub-01_ses-pre'
+    """
+
+    if not subject_id.startswith('sub-'):
+        subject_id = f'sub-{subject_id}'
+
+    if session_id:
+        ses_str = session_id
+        if isinstance(session_id, list):
+            from smriprep.utils.misc import stringify_sessions
+
+            ses_str = stringify_sessions(session_id)
+        if not ses_str.startswith('ses-'):
+            ses_str = f'ses-{ses_str}'
+        subject_id += f'_{ses_str}'
+    return subject_id
diff --git a/fmriprep/reports/core.py b/fmriprep/reports/core.py
index a7b04196f..4966ed891 100644
--- a/fmriprep/reports/core.py
+++ b/fmriprep/reports/core.py
@@ -67,7 +67,13 @@ def run_reports(
 
 
 def generate_reports(
-    subject_list, output_dir, run_uuid, session_list=None, bootstrap_file=None, work_dir=None
+    subject_list: list[str] | str,
+    output_dir: Path | str,
+    run_uuid: str,
+    session_list: list[str] | str | None = None,
+    bootstrap_file: Path | str | None = None,
+    work_dir: Path | str | None = None,
+    sessionwise: bool = False,
 ):
     """Generate reports for a list of subjects."""
     reportlets_dir = None
@@ -100,21 +106,22 @@ def generate_reports(
             bootstrap_file = data.load('reports-spec-anat.yml')
             html_report = f'sub-{subject_label}_anat.html'
 
-        report_error = run_reports(
-            output_dir,
-            subject_label,
-            run_uuid,
-            bootstrap_file=bootstrap_file,
-            out_filename=html_report,
-            reportlets_dir=reportlets_dir,
-            errorname=f'report-{run_uuid}-{subject_label}.err',
-            subject=subject_label,
-        )
-        # If the report generation failed, append the subject label for which it failed
-        if report_error is not None:
-            errors.append(report_error)
-
-        if n_ses > config.execution.aggr_ses_reports:
+        if not sessionwise:
+            report_error = run_reports(
+                output_dir,
+                subject_label,
+                run_uuid,
+                bootstrap_file=bootstrap_file,
+                out_filename=html_report,
+                reportlets_dir=reportlets_dir,
+                errorname=f'report-{run_uuid}-{subject_label}.err',
+                subject=subject_label,
+            )
+            # If the report generation failed, append the subject label for which it failed
+            if report_error is not None:
+                errors.append(report_error)
+
+        if (n_ses > config.execution.aggr_ses_reports) or sessionwise:
             # Beyond a certain number of sessions per subject,
             # we separate the functional reports per session
             if session_list is None:
@@ -124,11 +131,15 @@ def generate_reports(
                     subject=subject_label, **filters
                 )
 
-            session_list = [ses.removeprefix('ses-') for ses in session_list]
-
             for session_label in session_list:
-                bootstrap_file = data.load('reports-spec-func.yml')
-                html_report = f'sub-{subject_label}_ses-{session_label}_func.html'
+                session_label = session_label.removeprefix('ses-')
+                if sessionwise:
+                    # Include the anatomical as well
+                    bootstrap_file = data.load('reports-spec.yml')
+                    html_report = f'sub-{subject_label}_ses-{session_label}.html'
+                else:
+                    bootstrap_file = data.load('reports-spec-func.yml')
+                    html_report = f'sub-{subject_label}_ses-{session_label}_func.html'
 
                 report_error = run_reports(
                     output_dir,
diff --git a/fmriprep/utils/misc.py b/fmriprep/utils/misc.py
index 27b4a60fa..f3c803644 100644
--- a/fmriprep/utils/misc.py
+++ b/fmriprep/utils/misc.py
@@ -66,3 +66,28 @@ def estimate_bold_mem_usage(bold_fname: str) -> tuple[int, dict]:
     }
 
     return bold_tlen, mem_gb
+
+
+def fmt_subjects_sessions(subses: list[tuple[str]], concat_limit: int = 1):
+    """
+    Format a list of subjects and sessions to be printed.
+
+    Example
+    -------
+    >>> fmt_subjects_sessions([('01', 'A'), ('02', ['A', 'B']), ('03', None), ('04', ['A'])])
+    'sub-01 ses-A, sub-02 (2 sessions), sub-03, sub-04 ses-A'
+    """
+    output = []
+    for subject, session in subses:
+        if isinstance(session, list):
+            if len(session) > concat_limit:
+                output.append(f'sub-{subject} ({len(session)} sessions)')
+                continue
+            session = session[0]
+
+        if session is None:
+            output.append(f'sub-{subject}')
+        else:
+            output.append(f'sub-{subject} ses-{session}')
+
+    return ', '.join(output)
diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py
index a409a9ba5..1058c863b 100644
--- a/fmriprep/workflows/base.py
+++ b/fmriprep/workflows/base.py
@@ -52,7 +52,7 @@ def init_fmriprep_wf():
     Build *fMRIPrep*'s pipeline.
 
     This workflow organizes the execution of FMRIPREP, with a sub-workflow for
-    each subject.
+    each processing group.
 
     If FreeSurfer's ``recon-all`` is to be run, a corresponding folder is created
     and populated with any needed template subjects under the derivatives folder.
@@ -91,12 +91,23 @@ def init_fmriprep_wf():
         if config.execution.fs_subjects_dir is not None:
             fsdir.inputs.subjects_dir = str(config.execution.fs_subjects_dir.absolute())
 
-    for subject_id in config.execution.participant_label:
-        single_subject_wf = init_single_subject_wf(subject_id)
+    for subject_id, session_ids in config.execution.processing_groups:
+        log_dir = config.execution.fmriprep_dir / f'sub-{subject_id}'
+        sessions = listify(session_ids)
+        ses_str = ''
 
-        single_subject_wf.config['execution']['crashdump_dir'] = str(
-            config.execution.fmriprep_dir / f'sub-{subject_id}' / 'log' / config.execution.run_uuid
-        )
+        if isinstance(sessions, list):
+            from smriprep.utils.misc import stringify_sessions
+
+            ses_str = stringify_sessions(sessions)
+            log_dir /= f'ses-{ses_str}'
+
+        log_dir = log_dir / 'log' / config.execution.run_uuid
+
+        wf_name = '_'.join(['sub', subject_id, *(('ses', ses_str) if ses_str else ()), 'wf'])
+        single_subject_wf = init_single_subject_wf(subject_id, sessions, name=wf_name)
+
+        single_subject_wf.config['execution']['crashdump_dir'] = str(log_dir)
         for node in single_subject_wf._get_all_nodes():
             node.config = deepcopy(single_subject_wf.config)
         if freesurfer:
@@ -105,16 +116,17 @@ def init_fmriprep_wf():
             fmriprep_wf.add_nodes([single_subject_wf])
 
         # Dump a copy of the config file into the log directory
-        log_dir = (
-            config.execution.fmriprep_dir / f'sub-{subject_id}' / 'log' / config.execution.run_uuid
-        )
         log_dir.mkdir(exist_ok=True, parents=True)
         config.to_filename(log_dir / 'fmriprep.toml')
 
     return fmriprep_wf
 
 
-def init_single_subject_wf(subject_id: str):
+def init_single_subject_wf(
+    subject_id: str,
+    session_id: str | list[str] | None = None,
+    name: str | None = None,
+):
     """
     Organize the preprocessing pipeline for a single subject.
 
@@ -139,6 +151,12 @@ def init_single_subject_wf(subject_id: str):
     ----------
     subject_id : :obj:`str`
         Subject label for this single-subject workflow.
+    session_id
+        Session label(s) for this workflow.
+    name
+        Name of the workflow.
+        If not provided, will be set to ``sub_{subject_id}_ses_{session_id}_wf``.
+
 
     Inputs
     ------
@@ -151,7 +169,6 @@ def init_single_subject_wf(subject_id: str):
     from niworkflows.interfaces.nilearn import NILEARN_VERSION
     from niworkflows.interfaces.utility import KeySelect
     from niworkflows.utils.bids import collect_data
-    from niworkflows.utils.misc import fix_multi_T1w_source_name
     from niworkflows.utils.spaces import Reference
     from smriprep.workflows.anatomical import init_anat_fit_wf
     from smriprep.workflows.outputs import (
@@ -167,9 +184,13 @@ def init_single_subject_wf(subject_id: str):
         init_resample_surfaces_wf,
     )
 
+    from fmriprep.interfaces.bids import BIDSSourceFile, CreateFreeSurferID
     from fmriprep.workflows.bold.base import init_bold_wf
 
-    workflow = Workflow(name=f'sub_{subject_id}_wf')
+    if name is None:
+        name = f'sub_{subject_id}_wf'
+
+    workflow = Workflow(name=name)
     workflow.__desc__ = f"""
 Results included in this manuscript come from preprocessing
 performed using *fMRIPrep* {config.environment.version}
@@ -204,6 +225,7 @@ def init_single_subject_wf(subject_id: str):
     subject_data = collect_data(
         config.execution.layout,
         subject_id,
+        session_id=session_id,
         task=config.execution.task_id,
         echo=config.execution.echo_idx,
         bids_filters=config.execution.bids_filters,
@@ -255,6 +277,7 @@ def init_single_subject_wf(subject_id: str):
                     derivatives_dir=deriv_dir,
                     subject_id=subject_id,
                     std_spaces=std_spaces,
+                    session_id=session_id,
                 )
             )
 
@@ -265,15 +288,25 @@ def init_single_subject_wf(subject_id: str):
             subject_data=subject_data,
             anat_only=config.workflow.anat_only,
             subject_id=subject_id,
-            anat_derivatives=anatomical_cache if anatomical_cache else None,
+            anat_derivatives=anatomical_cache or None,
         ),
         name='bidssrc',
     )
 
+    src_file = pe.Node(
+        BIDSSourceFile(
+            precomputed=anatomical_cache,
+            sessionwise=config.workflow.subject_anatomical_reference == 'sessionwise',
+        ),
+        name='source_anatomical',
+    )
+
     bids_info = pe.Node(
         BIDSInfo(bids_dir=config.execution.bids_dir, bids_validate=False), name='bids_info'
     )
 
+    create_fs_id = pe.Node(CreateFreeSurferID(), name='create_fs_id')
+
     summary = pe.Node(
         SubjectSummary(
             std_spaces=spaces.get_spaces(nonstandard=False),
@@ -322,7 +355,7 @@ def init_single_subject_wf(subject_id: str):
         freesurfer=config.workflow.run_reconall,
         hires=config.workflow.hires,
         fs_no_resume=config.workflow.fs_no_resume,
-        longitudinal=config.workflow.longitudinal,
+        longitudinal=config.workflow.subject_anatomical_reference == 'unbiased',
         msm_sulc=msm_sulc,
         t1w=subject_data['t1w'],
         t2w=subject_data['t2w'],
@@ -342,17 +375,15 @@ def init_single_subject_wf(subject_id: str):
             'No T1w image found; using precomputed T1w image: %s', anatomical_cache['t1w_preproc']
         )
         workflow.connect([
-            (bidssrc, bids_info, [(('bold', fix_multi_T1w_source_name), 'in_file')]),
             (anat_fit_wf, summary, [('outputnode.t1w_preproc', 't1w')]),
             (anat_fit_wf, ds_report_summary, [('outputnode.t1w_preproc', 'source_file')]),
             (anat_fit_wf, ds_report_about, [('outputnode.t1w_preproc', 'source_file')]),
         ])  # fmt:skip
     else:
         workflow.connect([
-            (bidssrc, bids_info, [(('t1w', fix_multi_T1w_source_name), 'in_file')]),
-            (bidssrc, summary, [('t1w', 't1w')]),
-            (bidssrc, ds_report_summary, [(('t1w', fix_multi_T1w_source_name), 'source_file')]),
-            (bidssrc, ds_report_about, [(('t1w', fix_multi_T1w_source_name), 'source_file')]),
+            (src_file, summary, [('source_file', 't1w')]),
+            (src_file, ds_report_summary, [('source_file', 'source_file')]),
+            (src_file, ds_report_about, [('source_file', 'source_file')]),
         ])  # fmt:skip
 
     workflow.connect([
@@ -363,7 +394,13 @@ def init_single_subject_wf(subject_id: str):
             ('roi', 'inputnode.roi'),
             ('flair', 'inputnode.flair'),
         ]),
-        (bids_info, anat_fit_wf, [(('subject', _prefix), 'inputnode.subject_id')]),
+        (bidssrc, src_file, [('out_dict', 'bids_info')]),
+        (src_file, bids_info, [('source_file', 'in_file')]),
+        (bids_info, create_fs_id, [
+            ('subject', 'subject_id'),
+            ('session', 'session_id'),
+        ]),
+        (create_fs_id, anat_fit_wf, [('subject_id', 'inputnode.subject_id')]),
         # Reporting connections
         (inputnode, summary, [('subjects_dir', 'subjects_dir')]),
         (bidssrc, summary, [('t2w', 't2w'), ('bold', 'bold')]),
@@ -927,10 +964,6 @@ def map_fieldmap_estimation(
     return fmap_estimators, estimator_map
 
 
-def _prefix(subid):
-    return subid if subid.startswith('sub-') else f'sub-{subid}'
-
-
 def clean_datasinks(workflow: pe.Workflow) -> pe.Workflow:
     # Overwrite ``out_path_base`` of smriprep's DataSinks
     for node in workflow.list_node_names():
diff --git a/fmriprep/workflows/tests/test_base.py b/fmriprep/workflows/tests/test_base.py
index df96afdd1..6f61a19d1 100644
--- a/fmriprep/workflows/tests/test_base.py
+++ b/fmriprep/workflows/tests/test_base.py
@@ -117,6 +117,7 @@ def _make_params(
     freesurfer: bool = True,
     ignore: list[str] = None,
     force: list[str] = None,
+    subject_anatomical_reference: str = 'first-lex',
     bids_filters: dict = None,
 ):
     if ignore is None:
@@ -138,6 +139,7 @@ def _make_params(
         freesurfer,
         ignore,
         force,
+        subject_anatomical_reference,
         bids_filters,
     )
 
@@ -158,6 +160,7 @@ def _make_params(
         'freesurfer',
         'ignore',
         'force',
+        'subject_anatomical_reference',
         'bids_filters',
     ),
     [
@@ -189,6 +192,8 @@ def _make_params(
         # _make_params(freesurfer=False, bold2anat_init="header", force=['no-bbr']),
         # Regression test for gh-3154:
         _make_params(bids_filters={'sbref': {'suffix': 'sbref'}}),
+        _make_params(subject_anatomical_reference='unbiased'),
+        _make_params(subject_anatomical_reference='sessionwise'),
     ],
 )
 def test_init_fmriprep_wf(
@@ -208,6 +213,7 @@ def test_init_fmriprep_wf(
     freesurfer: bool,
     ignore: list[str],
     force: list[str],
+    subject_anatomical_reference: str,
     bids_filters: dict,
 ):
     with mock_config(bids_dir=bids_root):