diff --git a/docs/source/v1.10.md.inc b/docs/source/v1.10.md.inc
index f26568f34..eb75418dd 100644
--- a/docs/source/v1.10.md.inc
+++ b/docs/source/v1.10.md.inc
@@ -1,8 +1,8 @@
 ## v1.10.0 (unreleased)
 
-[//]: # (### :new: New features & enhancements)
+### :new: New features & enhancements
 
-[//]: # (- Whatever (#000 by @whoever))
+- It is now possible to use separate MRIs for each session within a subject, as in longitudinal studies. This is achieved by creating separate "subject" folders for each subject-session combination, with the naming convention `sub-XXX_ses-YYY`, in the freesurfer `SUBJECTS_DIR`. (#987 by @drammock)
 
 [//]: # (### :warning: Behavior changes)
 
@@ -25,6 +25,7 @@
 
 - Choose the theme (dark of light) automatically based on the user's operating system setting (#979 by @hoechenberger)
 - Bundle all previously-external JavaScript to better preserve users' privacy (#982 by @hoechenberger)
+
 ### :medical_symbol: Code health
 
 - Switch from using relative to using absolute imports (#969 by @hoechenberger)
diff --git a/mne_bids_pipeline/_config_utils.py b/mne_bids_pipeline/_config_utils.py
index fa3f2accb..b52ab95bc 100644
--- a/mne_bids_pipeline/_config_utils.py
+++ b/mne_bids_pipeline/_config_utils.py
@@ -36,13 +36,17 @@ def get_fs_subjects_dir(config: SimpleNamespace) -> pathlib.Path:
         return pathlib.Path(config.subjects_dir).expanduser().resolve()
 
 
-def get_fs_subject(config: SimpleNamespace, subject: str) -> str:
+def get_fs_subject(
+    config: SimpleNamespace, subject: str, session: str | None = None
+) -> str:
     subjects_dir = get_fs_subjects_dir(config)
 
     if config.use_template_mri is not None:
         return config.use_template_mri
 
-    if (pathlib.Path(subjects_dir) / subject).exists():
+    if session is not None:
+        return f"sub-{subject}_ses-{session}"
+    elif (pathlib.Path(subjects_dir) / subject).exists():
         return subject
     else:
         return f"sub-{subject}"
diff --git a/mne_bids_pipeline/steps/freesurfer/_01_recon_all.py b/mne_bids_pipeline/steps/freesurfer/_01_recon_all.py
index 26f84dda3..634740951 100755
--- a/mne_bids_pipeline/steps/freesurfer/_01_recon_all.py
+++ b/mne_bids_pipeline/steps/freesurfer/_01_recon_all.py
@@ -11,25 +11,33 @@
 
 from mne.utils import run_subprocess
 
-from mne_bids_pipeline._config_utils import get_fs_subjects_dir, get_subjects
+from mne_bids_pipeline._config_utils import (
+    get_fs_subjects_dir,
+    get_sessions,
+    get_subjects,
+)
 from mne_bids_pipeline._logging import gen_log_kwargs, logger
 from mne_bids_pipeline._parallel import get_parallel_backend, parallel_func
 
 fs_bids_app = Path(__file__).parent / "contrib" / "run.py"
 
 
-def run_recon(root_dir, subject, fs_bids_app, subjects_dir) -> None:
+def run_recon(root_dir, subject, fs_bids_app, subjects_dir, session=None) -> None:
     subj_dir = subjects_dir / f"sub-{subject}"
+    sub_ses = f"Subject {subject}"
+    if session is not None:
+        subj_dir = subj_dir.with_name(f"{subj_dir.name}_ses-{session}")
+        sub_ses = f"{sub_ses} session {session}"
 
     if subj_dir.exists():
         msg = (
-            f"Subject {subject} is already present. Please delete the "
+            f"Recon for {sub_ses} is already present. Please delete the "
             f"directory if you want to recompute."
         )
         logger.info(**gen_log_kwargs(message=msg))
         return
     msg = (
-        "Running recon-all on subject {subject}. This will take "
+        "Running recon-all on {sub_ses}. This will take "
         "a LONG time – it's a good idea to let it run over night."
     )
     logger.info(**gen_log_kwargs(message=msg))
@@ -56,10 +64,16 @@ def run_recon(root_dir, subject, fs_bids_app, subjects_dir) -> None:
         f"--license_file={license_file}",
         f"--participant_label={subject}",
     ]
+    if session is not None:
+        cmd += [f"--session_label={session}"]
     logger.debug("Running: " + " ".join(cmd))
     run_subprocess(cmd, env=env, verbose=logger.level)
 
 
+def _has_session_specific_anat(subject, session, subjects_dir):
+    return (subjects_dir / f"sub-{subject}_ses-{session}").exists()
+
+
 def main(*, config) -> None:
     """Run freesurfer recon-all command on BIDS dataset.
 
@@ -77,21 +91,33 @@ def main(*, config) -> None:
 
     You must have freesurfer available on your system.
 
-    Run via the MNE BIDS Pipeline's `run.py`:
+    Run via the MNE BIDS Pipeline's CLI:
 
-    python run.py --steps=freesurfer --config=your_pipeline_config.py
+    mne_bids_pipeline --steps=freesurfer --config=your_pipeline_config.py
 
     """  # noqa
     subjects = get_subjects(config)
+    sessions = get_sessions(config)
     root_dir = config.bids_root
     subjects_dir = Path(get_fs_subjects_dir(config))
     subjects_dir.mkdir(parents=True, exist_ok=True)
 
+    # check for session-specific MRIs within subject, and handle accordingly
+    subj_sess = list()
+    for _subj in subjects:
+        for _sess in sessions:
+            session = (
+                _sess
+                if _has_session_specific_anat(_subj, _sess, subjects_dir)
+                else None
+            )
+            subj_sess.append((_subj, session))
+
     with get_parallel_backend(config.exec_params):
         parallel, run_func = parallel_func(run_recon, exec_params=config.exec_params)
         parallel(
-            run_func(root_dir, subject, fs_bids_app, subjects_dir)
-            for subject in subjects
+            run_func(root_dir, subject, fs_bids_app, subjects_dir, session)
+            for subject, session in subj_sess
         )
 
         # Handle fsaverage
diff --git a/mne_bids_pipeline/steps/freesurfer/_02_coreg_surfaces.py b/mne_bids_pipeline/steps/freesurfer/_02_coreg_surfaces.py
index dbab7d037..c92ba581a 100644
--- a/mne_bids_pipeline/steps/freesurfer/_02_coreg_surfaces.py
+++ b/mne_bids_pipeline/steps/freesurfer/_02_coreg_surfaces.py
@@ -13,6 +13,7 @@
 from mne_bids_pipeline._config_utils import (
     get_fs_subject,
     get_fs_subjects_dir,
+    get_sessions,
     get_subjects,
 )
 from mne_bids_pipeline._logging import gen_log_kwargs, logger
@@ -81,9 +82,14 @@ def make_coreg_surfaces(
     )
 
 
-def get_config(*, config, subject) -> SimpleNamespace:
+def get_config(
+    *,
+    config: SimpleNamespace,
+    subject: str,
+    session: str | None = None,
+) -> SimpleNamespace:
     cfg = SimpleNamespace(
-        fs_subject=get_fs_subject(config, subject),
+        fs_subject=get_fs_subject(config, subject, session=session),
         fs_subjects_dir=get_fs_subjects_dir(config),
     )
     return cfg
@@ -92,6 +98,7 @@ def get_config(*, config, subject) -> SimpleNamespace:
 def main(*, config) -> None:
     # Ensure we're also processing fsaverage if present
     subjects = get_subjects(config)
+    sessions = get_sessions(config)
     if (Path(get_fs_subjects_dir(config)) / "fsaverage").exists():
         subjects.append("fsaverage")
 
@@ -105,10 +112,12 @@ def main(*, config) -> None:
                 cfg=get_config(
                     config=config,
                     subject=subject,
+                    session=session,
                 ),
                 exec_params=config.exec_params,
                 force_run=config.recreate_scalp_surface,
                 subject=subject,
             )
             for subject in subjects
+            for session in sessions
         )
diff --git a/mne_bids_pipeline/steps/source/_01_make_bem_surfaces.py b/mne_bids_pipeline/steps/source/_01_make_bem_surfaces.py
index f9bd60ec9..ff5bf4ff3 100644
--- a/mne_bids_pipeline/steps/source/_01_make_bem_surfaces.py
+++ b/mne_bids_pipeline/steps/source/_01_make_bem_surfaces.py
@@ -123,9 +123,10 @@ def get_config(
     *,
     config: SimpleNamespace,
     subject: str,
+    session: str | None = None,
 ) -> SimpleNamespace:
     cfg = SimpleNamespace(
-        fs_subject=get_fs_subject(config=config, subject=subject),
+        fs_subject=get_fs_subject(config=config, subject=subject, session=session),
         fs_subjects_dir=get_fs_subjects_dir(config=config),
         bem_mri_images=config.bem_mri_images,
         freesurfer_verbose=config.freesurfer_verbose,
@@ -160,12 +161,14 @@ def main(*, config: SimpleNamespace) -> None:
                 cfg=get_config(
                     config=config,
                     subject=subject,
+                    session=session,
                 ),
                 exec_params=config.exec_params,
                 subject=subject,
-                session=get_sessions(config)[0],
+                session=session,
                 force_run=config.recreate_bem,
             )
             for subject in get_subjects(config)
+            for session in get_sessions(config)
         )
     save_logs(config=config, logs=logs)
diff --git a/mne_bids_pipeline/steps/source/_02_make_bem_solution.py b/mne_bids_pipeline/steps/source/_02_make_bem_solution.py
index 5c74e2379..dda6dd2b1 100644
--- a/mne_bids_pipeline/steps/source/_02_make_bem_solution.py
+++ b/mne_bids_pipeline/steps/source/_02_make_bem_solution.py
@@ -12,6 +12,7 @@
     _get_bem_conductivity,
     get_fs_subject,
     get_fs_subjects_dir,
+    get_sessions,
     get_subjects,
 )
 from mne_bids_pipeline._logging import gen_log_kwargs, logger
@@ -23,6 +24,7 @@ def get_input_fnames_make_bem_solution(
     *,
     cfg: SimpleNamespace,
     subject: str,
+    session: str | None = None,
 ) -> dict:
     in_files = dict()
     conductivity, _ = _get_bem_conductivity(cfg)
@@ -37,6 +39,7 @@ def get_output_fnames_make_bem_solution(
     *,
     cfg: SimpleNamespace,
     subject: str,
+    session: str | None = None,
 ) -> dict:
     out_files = dict()
     bem_dir = Path(cfg.fs_subjects_dir) / cfg.fs_subject / "bem"
@@ -56,9 +59,10 @@ def make_bem_solution(
     exec_params: SimpleNamespace,
     subject: str,
     in_files: dict,
+    session: str | None = None,
 ) -> dict:
     msg = "Calculating BEM solution"
-    logger.info(**gen_log_kwargs(message=msg, subject=subject))
+    logger.info(**gen_log_kwargs(message=msg, subject=subject, session=session))
     conductivity, _ = _get_bem_conductivity(cfg)
     bem_model = mne.make_bem_model(
         subject=cfg.fs_subject,
@@ -81,9 +85,10 @@ def get_config(
     *,
     config: SimpleNamespace,
     subject: str,
+    session: str | None = None,
 ) -> SimpleNamespace:
     cfg = SimpleNamespace(
-        fs_subject=get_fs_subject(config=config, subject=subject),
+        fs_subject=get_fs_subject(config=config, subject=subject, session=session),
         fs_subjects_dir=get_fs_subjects_dir(config),
         ch_types=config.ch_types,
         use_template_mri=config.use_template_mri,
@@ -112,11 +117,13 @@ def main(*, config) -> None:
         )
         logs = parallel(
             run_func(
-                cfg=get_config(config=config, subject=subject),
+                cfg=get_config(config=config, subject=subject, session=session),
                 exec_params=config.exec_params,
                 subject=subject,
+                session=session,
                 force_run=config.recreate_bem,
             )
             for subject in get_subjects(config)
+            for session in get_sessions(config)
         )
     save_logs(config=config, logs=logs)
diff --git a/mne_bids_pipeline/steps/source/_03_setup_source_space.py b/mne_bids_pipeline/steps/source/_03_setup_source_space.py
index 54a604f1e..316e8290a 100644
--- a/mne_bids_pipeline/steps/source/_03_setup_source_space.py
+++ b/mne_bids_pipeline/steps/source/_03_setup_source_space.py
@@ -10,6 +10,7 @@
 from mne_bids_pipeline._config_utils import (
     get_fs_subject,
     get_fs_subjects_dir,
+    get_sessions,
     get_subjects,
 )
 from mne_bids_pipeline._logging import gen_log_kwargs, logger
@@ -71,11 +72,12 @@ def get_config(
     *,
     config: SimpleNamespace,
     subject: str,
+    session: str | None = None,
 ) -> SimpleNamespace:
     cfg = SimpleNamespace(
         spacing=config.spacing,
         use_template_mri=config.use_template_mri,
-        fs_subject=get_fs_subject(config=config, subject=subject),
+        fs_subject=get_fs_subject(config=config, subject=subject, session=session),
         fs_subjects_dir=get_fs_subjects_dir(config),
     )
     return cfg
@@ -92,6 +94,7 @@ def main(*, config: SimpleNamespace) -> None:
         subjects = [config.use_template_mri]
     else:
         subjects = get_subjects(config=config)
+    sessions = get_sessions(config=config)
 
     with get_parallel_backend(config.exec_params):
         parallel, run_func = parallel_func(
@@ -102,10 +105,12 @@ def main(*, config: SimpleNamespace) -> None:
                 cfg=get_config(
                     config=config,
                     subject=subject,
+                    session=session,
                 ),
                 exec_params=config.exec_params,
                 subject=subject,
             )
             for subject in subjects
+            for session in sessions
         )
     save_logs(config=config, logs=logs)
diff --git a/mne_bids_pipeline/steps/source/_04_make_forward.py b/mne_bids_pipeline/steps/source/_04_make_forward.py
index 16eda19a5..45167d8a7 100644
--- a/mne_bids_pipeline/steps/source/_04_make_forward.py
+++ b/mne_bids_pipeline/steps/source/_04_make_forward.py
@@ -260,7 +260,7 @@ def get_config(
         source_info_path_update=config.source_info_path_update,
         noise_cov=_sanitize_callable(config.noise_cov),
         ch_types=config.ch_types,
-        fs_subject=get_fs_subject(config=config, subject=subject),
+        fs_subject=get_fs_subject(config=config, subject=subject, session=session),
         fs_subjects_dir=get_fs_subjects_dir(config),
         t1_bids_path=t1_bids_path,
         landmarks_kind=landmarks_kind,
diff --git a/mne_bids_pipeline/steps/source/_05_make_inverse.py b/mne_bids_pipeline/steps/source/_05_make_inverse.py
index f850fe4ef..3a88ab194 100644
--- a/mne_bids_pipeline/steps/source/_05_make_inverse.py
+++ b/mne_bids_pipeline/steps/source/_05_make_inverse.py
@@ -166,6 +166,7 @@ def get_config(
     *,
     config: SimpleNamespace,
     subject: str,
+    session: str | None = None,
 ) -> SimpleNamespace:
     cfg = SimpleNamespace(
         source_info_path_update=config.source_info_path_update,
@@ -178,7 +179,7 @@ def get_config(
         inverse_method=config.inverse_method,
         noise_cov=_sanitize_callable(config.noise_cov),
         report_stc_n_time_points=config.report_stc_n_time_points,
-        fs_subject=get_fs_subject(config=config, subject=subject),
+        fs_subject=get_fs_subject(config=config, subject=subject, session=session),
         fs_subjects_dir=get_fs_subjects_dir(config),
         **_bids_kwargs(config=config),
     )
@@ -199,6 +200,7 @@ def main(*, config: SimpleNamespace) -> None:
                 cfg=get_config(
                     config=config,
                     subject=subject,
+                    session=session,
                 ),
                 exec_params=config.exec_params,
                 subject=subject,
diff --git a/mne_bids_pipeline/steps/source/_99_group_average.py b/mne_bids_pipeline/steps/source/_99_group_average.py
index 2c9acfec7..77b0dae15 100644
--- a/mne_bids_pipeline/steps/source/_99_group_average.py
+++ b/mne_bids_pipeline/steps/source/_99_group_average.py
@@ -232,7 +232,7 @@ def main(*, config: SimpleNamespace) -> None:
                 cfg=cfg,
                 exec_params=exec_params,
                 subject=subject,
-                fs_subject=get_fs_subject(config=cfg, subject=subject),
+                fs_subject=get_fs_subject(config=cfg, subject=subject, session=session),
                 session=session,
             )
             for subject in subjects