diff --git a/fitlins/cli/run.py b/fitlins/cli/run.py
index 89e59fad..24163e63 100755
--- a/fitlins/cli/run.py
+++ b/fitlins/cli/run.py
@@ -96,10 +96,13 @@ def get_parser():
                         help="use BOLD files with the provided description label")
 
     g_prep = parser.add_argument_group('Options for preprocessing BOLD series')
-    g_prep.add_argument('-s', '--smoothing', action='store', metavar="TYPE:FWHM",
-                        help="Smooth BOLD series with FWHM mm kernel prior to fitting. "
-                             "Valid types: iso (isotropic); "
-                             "e.g. `--smothing iso:5` will use an isotropic 5mm FWHM kernel")
+    g_prep.add_argument('-s', '--smoothing', action='store', metavar="FWHM[:LEVEL:[TYPE]]",
+                        help="Smooth BOLD series with FWHM mm kernel prior to fitting at LEVEL. "
+                             "Optional analysis LEVEL (default: l1) may be specified numerically "
+                             "(e.g., `l1`) or by name (`run`, `subject`, `session` or `dataset`). "
+                             "Optional smoothing TYPE (default: iso) must be one of: `iso` (isotropic). "
+                             "e.g., `--smoothing 5:dataset:iso` will perform a 5mm FWHM isotropic "
+                             "smoothing on subject-level maps, before evaluating the dataset level.")
 
     g_perfm = parser.add_argument_group('Options to handle performance')
     g_perfm.add_argument('--n-cpus', action='store', default=0, type=int,
diff --git a/fitlins/interfaces/nistats.py b/fitlins/interfaces/nistats.py
index 1ad1e451..70dbfecd 100644
--- a/fitlins/interfaces/nistats.py
+++ b/fitlins/interfaces/nistats.py
@@ -142,6 +142,7 @@ class SecondLevelModelInputSpec(BaseInterfaceInputSpec):
     variance_maps = traits.List(traits.List(File(exists=True)))
     stat_metadata = traits.List(traits.List(traits.Dict), mandatory=True)
     contrast_info = traits.List(traits.Dict, mandatory=True)
+    smoothing_fwhm = traits.Float(desc='Full-width half max (FWHM) in mm for smoothing in mask')
 
 
 class SecondLevelModelOutputSpec(TraitedSpec):
@@ -171,7 +172,11 @@ class SecondLevelModel(NistatsBaseInterface, SimpleInterface):
 
     def _run_interface(self, runtime):
         from nistats import second_level_model as level2
-        model = level2.SecondLevelModel()
+        smoothing_fwhm = self.inputs.smoothing_fwhm
+        if not isdefined(smoothing_fwhm):
+            smoothing_fwhm = None
+        model = level2.SecondLevelModel(smoothing_fwhm=smoothing_fwhm)
+
         effect_maps = []
         variance_maps = []
         stat_maps = []
diff --git a/fitlins/workflows/base.py b/fitlins/workflows/base.py
index c5f4377b..520e9d94 100644
--- a/fitlins/workflows/base.py
+++ b/fitlins/workflows/base.py
@@ -1,4 +1,5 @@
 from pathlib import Path
+import warnings
 from nipype.pipeline import engine as pe
 from nipype.interfaces import utility as niu
 # from nipype.interfaces import fsl
@@ -59,17 +60,36 @@ def init_fitlins_wf(bids_dir, derivatives, out_dir, analysis_level, space,
         name='getter')
 
     if smoothing:
-        smoothing_params = smoothing.split(':', 1)
-        if smoothing_params[0] != 'iso':
-            raise ValueError(f"Unknown smoothing type {smoothing_params[0]}")
-        smoothing_fwhm = float(smoothing_params[1])
+        smoothing_params = smoothing.split(':', 2)
+        # Convert old style and warn; this should turn into an (informative) error around 0.5.0
+        if smoothing_params[0] == 'iso':
+            smoothing_params = (smoothing_params[1], 'l1', smoothing_params[0])
+            warnings.warn(
+                "The format for smoothing arguments has changed. Please use "
+                f"{':'.join(smoothing_params)} instead of {smoothing}.", FutureWarning)
+        # Add defaults to simplify later logic
+        if len(smoothing_params) == 1:
+            smoothing_params.extend(('l1', 'iso'))
+        elif len(smoothing_params) == 2:
+            smoothing_params.append('iso')
+
+        smoothing_fwhm, smoothing_level, smoothing_type = smoothing_params
+        smoothing_fwhm = float(smoothing_fwhm)
+        if smoothing_type not in ('iso'):
+            raise ValueError(f"Unknown smoothing type {smoothing_type}")
+
+        # Check that smmoothing level exists in model
+        if (smoothing_level.startswith("l") and
+                int(smoothing_level.strip("l")) > len(model_dict)):
+            raise ValueError(f"Invalid smoothing level {smoothing_level}")
+        elif (smoothing_level not in
+                [step['Level'] for step in model_dict['Steps']]):
+            raise ValueError(f"Invalid smoothing level {smoothing_level}")
 
     l1_model = pe.MapNode(
         FirstLevelModel(),
         iterfield=['session_info', 'contrast_info', 'bold_file', 'mask_file'],
         name='l1_model')
-    if smoothing:
-        l1_model.inputs.smoothing_fwhm = smoothing_fwhm
 
     # Set up common patterns
     image_pattern = 'reports/[sub-{subject}/][ses-{session}/]figures/[run-{run}/]' \
@@ -161,6 +181,9 @@ def init_fitlins_wf(bids_dir, derivatives, out_dir, analysis_level, space,
 
         level = 'l{:d}'.format(ix + 1)
 
+        if smoothing and smoothing_level in (step, level):
+            model.inputs.smoothing_fwhm = smoothing_fwhm
+
         # TODO: No longer used at higher level, suggesting we can simply return
         # entities from loader as a single list
         select_entities = pe.Node(