diff --git a/aeolis/bed.py b/aeolis/bed.py
index 36e97741..c5bbe216 100644
--- a/aeolis/bed.py
+++ b/aeolis/bed.py
@@ -108,6 +108,14 @@ def initialize(s, p):
     else:
         s['mass'][:,:,:,:] = p['bedcomp_file'].reshape(s['mass'].shape)                
 
+
+    # Store mass for mixing (only for reset_initial, default is layer_average)
+    if p['method_mixing'] == 'reset_initial': 
+        if p['bedcomp_mixing_file'] is not None:
+            s['mass_mixing'][:,:,:,:] = p['bedcomp_mixing_file'].reshape(s['mass_mixing'].shape)
+        else:
+            s['mass_mixing'][:] = s['mass'][:]         
+
     # initialize masks
     for k, v in p.items():
         if k.endswith('_mask'):
@@ -185,8 +193,17 @@ def mixtoplayer(s, p):
                                        # .repeat(nx, axis=1) \
                                        # .repeat(nl, axis=2)
 
-            mass1 = np.nanmean(mass, axis=2, keepdims=True).repeat(nl, axis=2)
-            # mass2 = gd * s['thlyr'][:,:,:,np.newaxis].repeat(nf, axis=-1)
+            # Two different methods for mixing are implemented. 
+            # The first method is the default method, which is layer_average (default), which is used to average the grain size distribution over the layers that are above the depth of disturbance. 
+            # The second method is reset_initial, which is used to reset the initial grain size distribution in the top layer of the bed.
+            if p['method_mixing'] == 'reset_initial':
+                mass1 = s['mass_mixing'][:]
+            elif p['method_mixing'] == 'layer_average':
+                mass1 = np.nanmean(mass, axis=2, keepdims=True).repeat(nl, axis=2)
+            else:
+                logger.warning(f'Unknown method for mixing: {p["method_mixing"]}')
+                mass1 = np.nanmean(mass, axis=2, keepdims=True).repeat(nl, axis=2)
+
             mass = mass1 * f + mass * (1. - f)
 
         
diff --git a/aeolis/constants.py b/aeolis/constants.py
index a07fd23a..9a527ae1 100644
--- a/aeolis/constants.py
+++ b/aeolis/constants.py
@@ -146,6 +146,7 @@
     ),
     ('ny','nx','nlayers','nfractions') : (
         'mass',                             # [kg/m^2] Sediment mass in bed
+        'mass_mixing',                      # [kg/m^2] Sediment mass in bed used for mixing
     ),
 }
 
@@ -192,6 +193,7 @@
     'wave_file'                     : None,               # Filename of ASCII file with time series of wave heights
     'meteo_file'                    : None,               # Filename of ASCII file with time series of meteorlogical conditions
     'bedcomp_file'                  : None,               # Filename of ASCII file with initial bed composition
+    'bedcomp_mixing_file'           : None,               # Filename of ASCII file with initial bed composition
     'threshold_file'                : None,               # Filename of ASCII file with shear velocity threshold
     'fence_file'                    : None,               # Filename of ASCII file with sand fence location/height (above the bed)
     'ne_file'                       : None,               # Filename of ASCII file with non-erodible layer
@@ -307,6 +309,7 @@
     'boundary_gw'                   : 'no_flow',          # Landward groundwater boundary, dGw/dx = 0 (or 'static')
     'method_moist_threshold'        : 'belly_johnson',    # Name of method to compute wind velocity threshold based on soil moisture content
     'method_moist_process'          : 'infiltration',     # Name of method to compute soil moisture content(infiltration or surface_moisture)
+    'method_mixing'                 : 'layer_average',    # Name of method to mixing the sediment bed composition in the mixed layers
     'offshore_flux'                 : 0.,                 # [-] Factor to determine offshore boundary flux as a function of Q0 (= 1 for saturated flux , = 0 for noflux)
     'constant_offshore_flux'        : 0.,                 # [kg/m/s] Constant input flux at offshore boundary
     'onshore_flux'                  : 0.,                 # [-] Factor to determine onshore boundary flux as a function of Q0 (= 1 for saturated flux , = 0 for noflux)