diff --git a/nibabies/cli/parser.py b/nibabies/cli/parser.py
index 557ab6c0..60c63910 100644
--- a/nibabies/cli/parser.py
+++ b/nibabies/cli/parser.py
@@ -456,6 +456,19 @@ def _slice_time_ref(value, parser):
         default=None,
         help='Initialize the random seed for the workflow',
     )
+    g_conf.add_argument(
+        '--me-t2s-fit-method',
+        action='store',
+        default='curvefit',
+        choices=['curvefit', 'loglin'],
+        help=(
+            'The method by which to estimate T2* and S0 for multi-echo data. '
+            "'curvefit' uses nonlinear regression. "
+            "It is more memory intensive, but also may be more accurate, than 'loglin'. "
+            "'loglin' uses log-linear regression. "
+            'It is faster and less memory intensive, but may be less accurate.'
+        ),
+    )
 
     # Confounds options
     g_confounds = parser.add_argument_group('Specific options for estimating confounds')
diff --git a/nibabies/config.py b/nibabies/config.py
index 9babae3e..d76917f0 100644
--- a/nibabies/config.py
+++ b/nibabies/config.py
@@ -630,6 +630,8 @@ class workflow(_Config):
     use_syn_sdc = None
     """Run *fieldmap-less* susceptibility-derived distortions estimation
     in the absence of any alternatives."""
+    me_t2s_fit_method = 'curvefit'
+    """The method by which to estimate T2*/S0 for multi-echo data"""
 
 
 class loggers:
diff --git a/nibabies/workflows/bold/fit.py b/nibabies/workflows/bold/fit.py
index c7f8f9e4..4ee2a1be 100644
--- a/nibabies/workflows/bold/fit.py
+++ b/nibabies/workflows/bold/fit.py
@@ -882,6 +882,7 @@ def init_bold_native_wf(
             echo_times=echo_times,
             mem_gb=mem_gb['filesize'],
             omp_nthreads=config.nipype.omp_nthreads,
+            me_t2s_fit_method=config.workflow.me_t2s_fit_method,
             name='bold_t2smap_wf',
         )
 
diff --git a/nibabies/workflows/bold/t2s.py b/nibabies/workflows/bold/t2s.py
index 951b5b12..0bec409b 100644
--- a/nibabies/workflows/bold/t2s.py
+++ b/nibabies/workflows/bold/t2s.py
@@ -46,6 +46,7 @@ def init_bold_t2s_wf(
     echo_times: ty.Sequence[float],
     mem_gb: float,
     omp_nthreads: int,
+    me_t2s_fit_method: ty.Literal['curvefit', 'loglin'] = 'curvefit',
     name: str = 'bold_t2s_wf',
 ):
     r"""
@@ -89,14 +90,16 @@ def init_bold_t2s_wf(
     from niworkflows.interfaces.morphology import BinaryDilation
 
     workflow = Workflow(name=name)
-    if config.workflow.me_t2s_fit_method == 'curvefit':
+    if me_t2s_fit_method == 'curvefit':
         fit_str = (
             'nonlinear regression. '
             'The T2<sup>★</sup>/S<sub>0</sub> estimates from a log-linear regression fit '
             'were used for initial values'
         )
-    else:
+    elif me_t2s_fit_method == 'loglin':
         fit_str = 'log-linear regression'
+    else:
+        fit_str = f'unknown method: {me_t2s_fit_method}'
 
     workflow.__desc__ = f"""\
 A T2<sup>★</sup> map was estimated from the preprocessed EPI echoes, by voxel-wise fitting
@@ -116,7 +119,7 @@ def init_bold_t2s_wf(
     dilate_mask = pe.Node(BinaryDilation(radius=2), name='dilate_mask')
 
     t2smap_node = pe.Node(
-        T2SMap(echo_times=list(echo_times), fittype=config.workflow.me_t2s_fit_method),
+        T2SMap(echo_times=list(echo_times), fittype=me_t2s_fit_method),
         name='t2smap_node',
         mem_gb=2.5 * mem_gb * len(echo_times),
     )