diff --git a/docs/workflows.rst b/docs/workflows.rst
index 35e2e616c..04757cbb6 100644
--- a/docs/workflows.rst
+++ b/docs/workflows.rst
@@ -598,6 +598,12 @@ and optimally weighted combination of all supplied single echo time series.
 This optimally combined time series is then carried forward for all subsequent
 preprocessing steps.
 
+The method by which T2* and S0 are estimated is determined by the ``--me-t2s-fit-method`` parameter.
+The default method is "curvefit", which uses nonlinear regression to estimate T2* and S0.
+The other option is "loglin", which uses log-linear regression.
+The "loglin" option is faster and less memory intensive,
+but it may be less accurate than "curvefit".
+
 References
 ----------
 
diff --git a/fmriprep/cli/parser.py b/fmriprep/cli/parser.py
index 433836bbc..215dc77c4 100644
--- a/fmriprep/cli/parser.py
+++ b/fmriprep/cli/parser.py
@@ -359,6 +359,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."
+        ),
+    )
 
     g_outputs = parser.add_argument_group("Options for modulating outputs")
     g_outputs.add_argument(
diff --git a/fmriprep/config.py b/fmriprep/config.py
index d384250c6..fac24856c 100644
--- a/fmriprep/config.py
+++ b/fmriprep/config.py
@@ -571,6 +571,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/fmriprep/workflows/bold/t2s.py b/fmriprep/workflows/bold/t2s.py
index 2d3a26053..f0b087716 100644
--- a/fmriprep/workflows/bold/t2s.py
+++ b/fmriprep/workflows/bold/t2s.py
@@ -89,12 +89,19 @@ def init_bold_t2s_wf(
     from niworkflows.interfaces.morphology import BinaryDilation
 
     workflow = Workflow(name=name)
-    workflow.__desc__ = """\
+    if config.workflow.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:
+        fit_str = "log-linear regression"
+
+    workflow.__desc__ = f"""\
 A T2<sup>★</sup> map was estimated from the preprocessed EPI echoes, by voxel-wise fitting
 the maximal number of echoes with reliable signal in that voxel to a monoexponential signal
-decay model with nonlinear regression.
-The T2<sup>★</sup>/S<sub>0</sub> estimates from a log-linear regression fit were used for
-initial values.
+decay model with {fit_str}.
 The calculated T2<sup>★</sup> map was then used to optimally combine preprocessed BOLD across
 echoes following the method described in [@posse_t2s].
 The optimally combined time series was carried forward as the *preprocessed BOLD*.
@@ -109,7 +116,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)),
+        T2SMap(echo_times=list(echo_times), fittype=config.workflow.me_t2s_fit_method),
         name='t2smap_node',
         mem_gb=2.5 * mem_gb * len(echo_times),
     )