Dynamic revetment implementation

aeolis/bed.py

@@ -125,7 +125,16 @@ def initialize(s, p):
     # initialize threshold
     if p['threshold_file'] is not None:
         s['uth'] = p['threshold_file'][:,:,np.newaxis].repeat(nf, axis=-1)
-
+
+    # initialize sand and cobble layers for composite beaches
+    s['zsand'][:,:] = p['zsand_file']
+    s['dcob'][:,:] = p['dcob_file']
+    s['dsandcover'][:,:] = p['dsandcover_file']
+    s['doverlap'][:,:] = p['doverlap_file']
+    if p['process_bedupdate_comp'] is True and (p['zsand_file'] is None or p['dcob_file'] is None or p['dsandcover_file'] is None or p['doverlap_file'] is None):
+        logger.log_and_raise('Process bedupdate for composite beaches is turned on but no initial sand and cobble locations are provided. Please'
+                             'provide a zsand_file, dcob_file, dsandcover_file and doverlap_file', exc=ValueError) 
+
     return s
 
 
@@ -342,8 +351,12 @@ def update(s, p):
         # s['dzb'] = dm[:, 0].reshape((ny + 1, nx + 1))
         s['dzb'] = dz.copy()
 
+        if p['process_bedupdate_comp']: 
+            s = update_composite(s, p)
+        else:
+            s['zb'] += dz
+
         # redistribute sediment from inactive zone to marine interaction zone
-        s['zb'] += dz
         if p['process_tide']:
             s['zs'] += dz #???
 
@@ -535,5 +548,173 @@ def arrange_layers(m,dm,d,nl,ix_ero,ix_dep):
     m[ix_dep,-1,:] -= dm[ix_dep,:] * d[ix_dep,-1,:]
 
     return m
-
+
+def update_composite(s, p):
+    '''Update bed and sand level and cobble infilling for composite beaches
+
+    Update XX by YYY
+    layers. 
+    Initial sand level, cobble thickness, overlap and sandcover are defined 
+    in the model configuration file by ``zsand``, ``dcob``, ``doverlap`` and
+    ``dsandcover``. The bathymetry is updated if the sand cover increases.
+
+    Parameters
+    ----------
+    s : dict
+        Spatial grids
+    p : dict
+        Model configuration parameters
+
+    Returns
+    -------
+    dict
+        Spatial grids
+
+    '''
+
+    dz = s['dzb']
+    por = p['porosity'] # Assumes that cobble porosity is the same as sand porosity
+
+    # FIND LOCATIONS
+    # with deposition
+    ix_depo = dz >0
+    # with erosion
+    ix_ero = dz < 0
+    # with 0 bed level change
+    ix_none = dz == 0.
+
+    # where no cobble is present, i.e. pure sand
+    ix_nocob = s['dcob'] == 0.
+    # where open cobbles are present
+    ix_cob = s['dcob'] > 0.
+    # where sand cover is present
+    ix_sc = s['dsandcover'] > 0.
+    # where no sand cover is present
+    ix_nosc = s['dsandcover'] == 0.
+    # where elevation change is bigger than sand cover, these are locations that can accommodate erosion
+    ix_sc_accom = dz + s['dsandcover'] >= 0.
+    # where elevation change is smaller than sand cover, these are locations that cannot accommodate erosion
+    ix_sc_noaccom = dz + s['dsandcover'] < 0.
+    # where open cobbles are present at the top
+    ix_opencob = s['dcob']-s['doverlap'] > 0.
+    # where elevation change is smaller than pore space in open cobble, these are locations that can accommodate deposition
+    ix_oc_accom = (s['dcob']-s['doverlap']) - dz/(1-por) >= 0.
+    # where elevation change is bigger than pore space in open cobble, these are locations that can't accommodate deposition
+    ix_oc_noaccom = (s['dcob']-s['doverlap']) - dz/(1-por) < 0.
+    # where cobbles are fully filled with sand, so doverlap == dcob
+    ix_fullcob = s['dcob']==s['doverlap']
+    # where filled cobbles are present at the top, dcob > 0 & doverlap == dcob & dsandcover == 0 
+    ix_fillcob = (ix_cob & ix_fullcob & ix_nosc)
+
+    # where no bed level change occurs and cobble is present
+    ix_step0 = ix_none & ix_cob
+    # where deposition takes place and sand cover is present
+    ix_step1 = ix_depo & ix_sc
+    # where deposition takes place and filled cobble is present  
+    ix_step2 = ix_depo & ix_cob & ix_fillcob
+    # where deposition takes place and open cobble is present with enough accommodation space
+    ix_step3 = ix_depo & ix_cob & ix_opencob & ix_oc_accom
+    # where deposition takes place and open cobble is present without enough accommodation space
+    ix_step4 = ix_depo & ix_cob & ix_opencob & ix_oc_noaccom
+    # where erosion takes place and open cobble is present
+    ix_step5 = ix_ero & ix_cob & ix_opencob
+    # where erosion takes place and filled cobble is present
+    ix_step6 = ix_ero & ix_cob & ix_fillcob
+    # where erosion takes place, sand cover is present and there is enough accommodation space
+    ix_step7 = ix_ero & ix_sc & ix_sc_accom
+    # where erosion takes place, sand cover is present and there is not enough accommodation space
+    ix_step8 = ix_ero & ix_sc & ix_sc_noaccom
+
+    # Check if all locations are assigned to a single step
+    ix_steps = np.vstack((ix_nocob[1,:], ix_step0[1,:], ix_step1[1,:], ix_step2[1,:], ix_step3[1,:], ix_step4[1,:], ix_step5[1,:], ix_step6[1,:], ix_step7[1,:], ix_step8[1,:]))
+    check = np.sum(ix_steps, axis=0)
+    if np.sum(s['dcob'] - s['doverlap'] < 0) > 0:
+        logger.warning('Cobble thickness is smaller than overlap thickness')
+    if np.sum(check!=1) > 0:
+        args = np.argwhere(check!=1)
+        nocob_check = ix_nocob[1,args]
+        step0_check = ix_step0[1,args]
+        step1_check = ix_step1[1,args]
+        step2_check = ix_step2[1,args]
+        step3_check = ix_step3[1,args]
+        step4_check = ix_step4[1,args]
+        step5_check = ix_step5[1,args]
+        step6_check = ix_step6[1,args]
+        step7_check = ix_step7[1,args]
+        step8_check = ix_step8[1,args]
+        print(check)
+
+    # if no cobble is present, 
+    # change sand level. This is independent of erosion or deposition
+    s['zsand'][ix_nocob] += dz[ix_nocob]
+
+    # if no bed level change occurs and cobble is present,
+    # keep all variables the same.
+    s['zsand'][ix_step0] = s['zsand'][ix_step0]
+
+    # if deposition takes place and sand cover is present 
+    # change sand cover and sand level. 
+    s['dsandcover'][ix_step1] += dz[ix_step1]
+    s['zsand'][ix_step1] += dz[ix_step1]
+
+    # if deposition takes place and filled cobble is present 
+    # change sand cover and sand level. 
+    s['dsandcover'][ix_step2] += dz[ix_step2]
+    s['zsand'][ix_step2] += dz[ix_step2]
+
+    # if deposition takes place and open cobble is present with enough accommodation space, 
+    # change sand level and overlap. 
+    s['zsand'][ix_step3] += dz[ix_step3]/(1-por)
+    s['doverlap'][ix_step3] += dz[ix_step3]/(1-por)
+
+    if np.sum(s['dcob'] - s['doverlap'] < 0) > 0:
+        print('Cobble thickness is smaller than overlap thickness')
+
+    # if deposition takes place and open cobble is present without enough accommodation space, 
+    # change sand cover, sand level and overlap. 
+    s['dsandcover'][ix_step4] += dz[ix_step4] - (s['dcob'][ix_step4] - s['doverlap'][ix_step4])*(1-por)
+    s['zsand'][ix_step4] += (s['dcob'][ix_step4] - s['doverlap'][ix_step4]) + s['dsandcover'][ix_step4]
+    s['doverlap'][ix_step4] = s['dcob'][ix_step4]
+
+    if np.sum(s['dcob'] - s['doverlap'] < 0) > 0:
+        np.argwhere(s['dcob'] - s['doverlap'] < 0)
+        print('Cobble thickness is smaller than overlap thickness')
+
+    # if erosion takes place and open cobble is present, 
+    # change sand level and overlap. 
+    s['zsand'][ix_step5] += dz[ix_step5]/(1-por)
+    s['doverlap'][ix_step5] += dz[ix_step5]/(1-por)
+
+    if np.sum(s['dcob'] - s['doverlap'] < 0) > 0:
+        print('Cobble thickness is smaller than overlap thickness')
+
+    # if erosion takes place and filled cobble is present, 
+    # change sand level and overlap. 
+    s['zsand'][ix_step6] += dz[ix_step6]/(1-por)
+    s['doverlap'][ix_step6] += dz[ix_step6]/(1-por)
+
+    # if erosion takes place, sand cover is present and there is enough accommodation space, 
+    # change sand cover and sand level. 
+    s['dsandcover'][ix_step7] += dz[ix_step7]
+    s['zsand'][ix_step7] += dz[ix_step7]
+
+    # print(s['dcob'][0,380], s['doverlap'][0,380], s['dsandcover'][0,380], s['zsand'][0,380], dz[0,380])
+
+    # if erosion takes place, sand cover is present and there is not enough accommodation space, 
+    # change overlap, sand level and sand cover. 
+    s['doverlap'][ix_step8] += (dz[ix_step8] + s['dsandcover'][ix_step8])/(1-por)
+    s['zsand'][ix_step8] += (dz[ix_step8] + s['dsandcover'][ix_step8])/(1-por) - s['dsandcover'][ix_step8]
+    s['dsandcover'][ix_step8] = 0
+
+    # print(s['dcob'][0,380], s['doverlap'][0,380], s['dsandcover'][0,380], s['zsand'][0,380], dz[0,380])
+
+    if np.sum(s['dcob'] - s['doverlap'] < 0) > 0:
+        arg = np.argwhere(s['dcob'] - s['doverlap'] < 0)
+        print('Cobble thickness is smaller than overlap thickness')
+    # for now, assume not enough erosion takes place to halt erosion through changing the roughness. Requires small enough time step.
+
+    # For all locations calculate the new bed level
+    s['zb'] = s['zsand'] + (s['dcob']-s['doverlap'])
+
+    return s
 

aeolis/constants.py

@@ -78,6 +78,10 @@
         'dzb',                              # [m/dt] Bed level change per time step (computed after avalanching!)
         'dzbyear',                          # [m/yr] Bed level change translated to m/y
         'dzbavg',                           # [m/year] Bed level change averaged over collected time steps
+        'zsand',                            # [m] Level of sand above reference
+        'dcob',                             # [m] Thickness of cobble layer
+        'dsandcover',                      # [m] Thickness of sand layer on top of cobble
+        'doverlap',                         # [m] Overlap between sand an cobble layer
         'S',                                # [-] Level of saturation
         'moist',                            # [-] Moisture content (volumetric)
         'moist_swr',                        # [-] Moisture content soil water retention relationship (volumetric)
@@ -158,14 +162,16 @@
     'process_wind'                  : True,               # Enable the process of wind
     'process_transport'             : True,               # Enable the process of transport
     'process_bedupdate'             : True,               # Enable the process of bed updating
+    'process_bedupdate_comp'        : False,               # Enable the process of bed updating for composite beaches
     'process_threshold'             : True,               # Enable the process of threshold
     'th_grainsize'                  : True,               # Enable wind velocity threshold based on grainsize
     'th_bedslope'                   : False,              # Enable wind velocity threshold based on bedslope
     'th_moisture'                   : False,              # Enable wind velocity threshold based on moisture
     'th_drylayer'                   : False,              # Enable threshold based on drying of layer
     'th_humidity'                   : False,              # Enable wind velocity threshold based on humidity
     'th_salt'                       : False,              # Enable wind velocity threshold based on salt
-    'th_sheltering'                 : False,              # Enable wind velocity threshold based on sheltering by roughness elements
+    'th_sheltering'                 : False,               # Enable wind velocity threshold based on sheltering by roughness elements using the grain size distrbution
+    'th_sedtrapping'                : False,               # Enable wind velocity threshold based on sediment trapping by roughness elements based on roughness density (lambda)
     'th_nelayer'                    : False,              # Enable wind velocity threshold based on a non-erodible layer
     'process_avalanche'             : False,              # Enable the process of avalanching
     'process_shear'                 : False,              # Enable the process of wind shear
@@ -201,6 +207,10 @@
     '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
     'veg_file'                      : None,               # Filename of ASCII file with initial vegetation density
+    'zsand_file'                    : None,               # Filename of ASCII file with sand level
+    'dcob_file'                     : None,               # Filename of ASCII file with cobble layer thickness
+    'dsandcover_file'              : None,               # Filename of ASCII file with thickness of sand cover on top of cobble layer
+    'doverlap_file'                  : None,               # Filename of ASCII file with overlap between sand level and cobble layer
     'supply_file'                   : None,               # Filename of ASCII file with a manual definition of sediment supply (mainly used in academic cases)
     'wave_mask'                     : None,               # Filename of ASCII file with mask for wave height
     'tide_mask'                     : None,               # Filename of ASCII file with mask for tidal elevation
@@ -254,7 +264,7 @@
     'Ck'                            : 2.78,               # [-] Constant in kawamura formulation for equilibrium sediment concentration
     'Cl'                            : 6.7,                # [-] Constant in lettau formulation for equilibrium sediment concentration
     'Cdk'                           : 5.,                 # [-] Constant in DK formulation for equilibrium sediment concentration
-    # 'm'                             : 0.5,                # [-] Factor to account for difference between average and maximum shear stress
+    'm'                             : 0.5,                # [-] Factor to account for difference between average and maximum shear stress
 #    'alpha'                         : 0.4,               # [-] Relation of vertical component of ejection velocity and horizontal velocity difference between impact and ejection
     'kappa'                         : 0.41,               # [-] Von Kármán constant
     'sigma'                         : 4.2,                # [-] Ratio between basal area and frontal area of roughness elements
@@ -270,7 +280,6 @@
     'facDOD'                        : .1,                 # [-] Ratio between depth of disturbance and local wave height
     'csalt'                         : 35e-3,              # [-] Maximum salt concentration in bed surface layer
     'cpair'                         : 1.0035e-3,          # [MJ/kg/oC] Specific heat capacity air
-
     'fc'                            : 0.11,               # [-] Moisture content at field capacity (volumetric)
     'w1_5'                          : 0.02,               # [-] Moisture content at wilting point (gravimetric)
     'resw_moist'                    : 0.01,               # [-] Residual soil moisture content (volumetric) 
@@ -332,6 +341,7 @@
     'alfa'                          : 0,                  # [deg] Real-world grid cell orientation wrt the North (clockwise)
     'dune_toe_elevation'            : 3,                  # Choose dune toe elevation, only used in the PH12 dune erosion solver
     'beach_slope'                   : 0.1,                # Define the beach slope, only used in the PH12 dune erosion solver
+    'method_runup'                  : 'stockdon',         # Define the equation used in the wave runup calculations
     'veg_min_elevation'             : 3,                  # Choose the minimum elevation where vegetation can grow
     'vegshear_type'                 : 'raupach',          # Choose the Raupach grid based solver (1D or 2D) or the Okin approach (1D only)
     'okin_c1_veg'                   : 0.48,               #x/h spatial reduction factor in Okin model for use with vegetation

aeolis/hydro.py

@@ -127,7 +127,11 @@ def interpolate(s, p, t):
                     tp = s['Tp'][iy][0]
                     wl = s['SWL'][iy][0]
 
-                    eta, sigma_s, R = calc_runup_stockdon(hs, tp, p['beach_slope'])
+                    if p['method_runup'] == 'stockdon':
+                        eta, sigma_s, R = calc_runup_stockdon(hs, tp, p['beach_slope'])
+                    elif p['method_runup'] == 'ruggiero':
+                        eta, sigma_s, R = calc_runup_ruggiero(hs)
+
                     s['R'][iy][:] = R
                     s['eta'][iy][:] = eta
                     s['sigma_s'][iy][:] = sigma_s
@@ -156,11 +160,15 @@ def interpolate(s, p, t):
         s['Tp'] = apply_mask(s['Tp'], s['wave_mask'])
 
 
+
     if p['process_runup']:
         ny = p['ny']
         if ('Hs' in p['external_vars']):
 
-            eta, sigma_s, R = calc_runup_stockdon(s['Hs'], s['Tp'], p['beach_slope'])
+            if p['method_runup'] == 'stockdon':
+                eta, sigma_s, R = calc_runup_stockdon(s['Hs'], s['Tp'], p['beach_slope'])
+            if p['method_runup'] == 'ruggiero':
+                eta, sigma_s, R = calc_runup_ruggiero(s['Hs'])
             s['R'][:] = R
 
             if hasattr(s['runup_mask'], "__len__"):
@@ -820,6 +828,23 @@ def calc_runup_stockdon(Ho, Tp, beta):
     return eta, sigma_s, R
 
 
+def calc_runup_ruggiero(Ho):
+    """
+    Calculate runup according to /Ruggiero et al 2004.
+    """
+    if Ho > 0:
+        R = 0.33 * Ho + 0.33 #formula for dissipative conditions
+        eta = 0
+        sigma_s = 0
+        # print(Ho, R)
+
+    else:
+        R = 0
+        eta = 0
+        sigma_s = 0
+
+    return eta, sigma_s, R
+
 
 
 

aeolis/shear.py

@@ -197,8 +197,13 @@ def __call__(self, x, y, z, taux, tauy, u0, udir,
         gc['z'] = self.interpolate(gi['x'], gi['y'], gi['z'], gc['x'], gc['y'], 0)
 
         # Project the taus0 and taun0 on the computational grid
-        gc['taux'] = np.full(np.shape(gc['x']), taus0)
-        gc['tauy'] = np.full(np.shape(gc['x']), taun0)
+        if np.all(taus0 == taus0[0,0]):
+            gc['taux'] = np.full(np.shape(gc['x']), taus0[0,0])
+            gc['tauy'] = np.full(np.shape(gc['x']), taun0[0,0])
+        else:
+            gc['taux'] = self.interpolate(gi['x'], gi['y'], taus0, gc['x'], gc['y'], 0)
+            gc['tauy'] = self.interpolate(gi['x'], gi['y'], taun0, gc['x'], gc['y'], 0)
+
 
         if plot:
             self.plot(ax=axs[0,1], cmap='Reds', stride=10, computational_grid=True)
@@ -261,8 +266,8 @@ def __call__(self, x, y, z, taux, tauy, u0, udir,
         # =====================================================================
 
         # Interpolate wind shear results to real grid
-        gi['taux'] = self.interpolate(gc['x'], gc['y'], gc['taux'], gi['x'], gi['y'], taus0)
-        gi['tauy'] = self.interpolate(gc['x'], gc['y'], gc['tauy'], gi['x'], gi['y'], taun0) 
+        gi['taux'] = self.interpolate(gc['x'], gc['y'], gc['taux'], gi['x'], gi['y'], taus0[0,0])
+        gi['tauy'] = self.interpolate(gc['x'], gc['y'], gc['tauy'], gi['x'], gi['y'], taus0[0,0])
 
         if process_separation:
             gi['hsep'] = self.interpolate(gc['x'], gc['y'], gc['hsep'], gi['x'], gi['y'], 0. )
@@ -556,7 +561,7 @@ def compute_shear(self, u0, nfilter=(1., 2.)):
 
         # Inner layer height
         l  = self.l    
- 
+
         # interpolate roughness length z0 to computational grid
         if np.size(z0)>1:
             z0new = self.interpolate(gi['x'], gi['y'], z0, gc['x'], gc['y'], 0)
@@ -836,13 +841,13 @@ def rotate(x, y, alpha, origin=(0,0)):
     def interpolate(self, x, y, z, xi, yi, z0):
         '''Interpolate one grid to an other'''
 
-        
+
         # First compute angle with horizontal
         dx = x[0,1] - x[0,0]
         dy = y[0,1] - y[0,0]
 
         angle = np.rad2deg(np.arctan(dy/dx))
-        
+
         if dx <= 0 and dy<=0:
             angle += 180.
 
@@ -877,7 +882,6 @@ def interpolate(self, x, y, z, xi, yi, z0):
             inter = scipy.interpolate.RegularGridInterpolator((y_pad[:,0].copy(order='C'), x_pad[0,:].copy(order='C')), z_pad, bounds_error = False, fill_value = z0)
             zi = inter(xyi).reshape(xi.shape)
 
-
         return zi