Add mixing option

Before, after hydrodynamic mixing, the bed composition of the mixed layers was reset to the average of those layers.

On shorter term scale, this approach seems to work, but after longer simulation periods, the approach appears to result in an underestimation of fine particles emerging in the top layer. If armouring occurs for a long period, the average bed composition of all layers starts to become more coarse. After several years, during mixing no small particles were emerging to the top layer anymore.

To solve this, an optional mixing method could be implemented that does not reset the mixed bed to the current average of the layers, but to the original bed composition.

method_mixing = layer_average (default)

OR

method_mixing = reset_initial (new)

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)
 
 

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)