From 8ed4e50153f017463cce3859d1220f89d3fa075c Mon Sep 17 00:00:00 2001 From: bartvanwesten Date: Wed, 24 Jul 2024 23:15:50 +0200 Subject: [PATCH] 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 | 21 +++++++++++++++++++-- aeolis/constants.py | 3 +++ 2 files changed, 22 insertions(+), 2 deletions(-) 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)