diff --git a/opendrift/models/leeway.py b/opendrift/models/leeway.py
index d66f18681..dcfa2c19d 100644
--- a/opendrift/models/leeway.py
+++ b/opendrift/models/leeway.py
@@ -154,6 +154,15 @@ class Leeway(OpenDriftSimulation):
         'y_sea_water_velocity': {
             'fallback': None
         },
+        'sea_surface_wave_stokes_drift_x_velocity': {
+            'fallback': 0,
+            'important': False
+        },
+        'sea_surface_wave_stokes_drift_y_velocity': {
+            'fallback': 0,
+            'important': False
+        },
+
         'land_binary_mask': {
             'fallback': None
         },
@@ -268,6 +277,10 @@ def __init__(self, d=None, *args, **kwargs):
                 'units': 'm/s',
                 'level': CONFIG_LEVEL_BASIC
             },
+                'drift:stokes_drift': {'type': 'bool', 'default': False,
+                    'description': 'Advection elements with surface Stokes drift (wave orbital motion). Note that this is originally considered to be implicit in Leeway coefficients.',
+                    'level': CONFIG_LEVEL_ADVANCED},
+
         })
 
         self._set_config_default('general:time_step_minutes', 10)
@@ -472,6 +485,12 @@ def update(self):
         logger.debug('Jibing %i out of %i elements.' %
                      (np.sum(jib), self.num_elements_active()))
 
+        # Move elements with Stokes drift, if activated.
+        # Note: Stokesdrift is originally considered to be implicit in Leeway coefficients,
+        # however, this study by Sutherland (2024) indicates it should be added explicitly:
+        # https://link.springer.com/article/10.1007/s10236-024-01600-3
+        self.stokes_drift()
+
     def export_ascii(self, filename):
         '''Export output to ASCII format of original version'''
 
diff --git a/opendrift/models/physics_methods.py b/opendrift/models/physics_methods.py
index b33076ba1..ad76febce 100644
--- a/opendrift/models/physics_methods.py
+++ b/opendrift/models/physics_methods.py
@@ -942,7 +942,7 @@ def stokes_drift(self, factor=1):
             logger.debug('Stokes drift is available, but not Tp: using Tp=8 for Stokes profile')
             wave_period = 8
 
-        stokes_profile = self.get_config('drift:stokes_drift_profile')
+        stokes_profile = self.get_config('drift:stokes_drift_profile', default='monochromatic')
         if stokes_profile == 'monochromatic':
             stokes_u, stokes_v, s = stokes_drift_profile_monochromatic(
                 self.environment.sea_surface_wave_stokes_drift_x_velocity,