Long term dune growth (by Selwyn)

aeolis/avalanching.py

@@ -240,4 +240,4 @@ def avalanche_loop(zb, zne, ds, dn, nx, ny, E, max_iter_ava, tan_dyn):
         # update bed level without non-erodible layer       
         zb += E * (q_in - q_out)
 
-    return zb, grad_h
+    return zb, grad_h

aeolis/bed.py

@@ -32,6 +32,7 @@
 import aeolis.gridparams
 from matplotlib import pyplot as plt
 from numba import njit
+import math
 
 # package modules
 from aeolis.utils import *
@@ -125,7 +126,7 @@ def initialize(s, p):
     # initialize threshold
     if p['threshold_file'] is not None:
         s['uth'] = p['threshold_file'][:,:,np.newaxis].repeat(nf, axis=-1)
-        
+
     return s
 
 
@@ -197,14 +198,40 @@ def mixtoplayer(s, p):
 
 
             s['mass'][ix] = mass[ix]
-
     return s
 
 
-def wet_bed_reset(s, p):
-    ''' Reset wet bed to initial bed level if the total water level is above the bed level.
+# def wet_bed_reset(s, p): **** Now part of def sediment_supply ****
+#     ''' Reset wet bed to initial bed level if the total water level is above the bed level.
+
+
+
+#     Parameters
+#     ----------
+#     s : dict
+#         Spatial grids
+#     p : dict
+#         Model configuration parameters
+
+#     Returns
+#     -------
+#     dict
+#         Spatial grids
+
+#     '''
+
+#     if p['process_wet_bed_reset']:
+
+#         Tbedreset = p['dt_opt'] / p['Tbedreset']
+
+#         ix = s['TWL'] > (s['zb'])
+#         s['zb'][ix] += (s['zb0'][ix] - s['zb'][ix]) * Tbedreset
+
+#     return s
 
 
+def wet_supply(s, p):
+    ''' Increase elevation of beach topography.
 
     Parameters
     ----------
@@ -219,19 +246,96 @@ def wet_bed_reset(s, p):
         Spatial grids
 
     '''
-
-    if p['process_wet_bed_reset']:
-
-        Tbedreset = p['dt_opt'] / p['Tbedreset']
 
-        ix = s['TWL'] > (s['zb'])
-        s['zb'][ix] += (s['zb0'][ix] - s['zb'][ix]) * Tbedreset
+    if p['process_wet_supply'] or p['process_wet_bed_reset']:
+
+        if p['method_wet_supply'] == 'wet_bed_reset':            
+            Tbedreset = p['dt_opt'] / p['Tbedreset']
 
-    return s
+            ix = s['TWL'] > (s['zb'])
+            s['zb'][ix] += (s['zb0'][ix] - s['zb'][ix]) * Tbedreset
+
+    if p['process_wet_supply']:
+
+        if p['method_wet_supply'] == 'vertical_beach_growth':
+            beach_inc = p['shoreline_change_rate']*math.cos((math.pi/2)-math.atan(p['beach_slope']))
+            vrate = (beach_inc*(1/365.25/24/3600))*p['dt'] #(m/timestep)
+            ny, nx = s['zb'].shape  
+
+            for iy in range(ny):
+                x_all = s['x'][iy,:]
+                zb_all = s['zb'][iy,:]               
+                xi =  (zb_all < p['dune_toe_elevation'])
+                beach_z = zb_all[xi]
+                b = beach_z[0] + vrate
+                x = x_all[xi] 
+                slope = p['beach_slope']
+                new_beach = slope*x + b 
+                s['zb'][iy,xi] = new_beach    
+
+        if p['method_wet_supply'] == 'constant_SCR_constant_tanB':
+
+            beach_inc = p['shoreline_change_rate']*math.cos((math.pi/2)-math.atan(p['beach_slope']))
+            vrate = (beach_inc/(365.25*24*3600))*p['dt'] #(m/timestep)
+            ny, nx = s['zb'].shape  
+
+            for iy in range(ny):
 
+                x_all = s['x'][iy,:]
+                zb_all = s['zb'][iy,:]
+
+                xi = zb_all < p['dune_toe_elevation']
+                beach_z = zb_all[xi]
+                x = x_all[xi]
+
+                xi3 = np.where(beach_z > p['zshoreline'])
+                xi3 = xi3[0][0]
+                b = beach_z[xi3] + vrate
+
+                new_temp_beach = p['beach_slope']*(x-x[xi3]) + b
+
+                xi2 = new_temp_beach <= np.min(zb_all)
+                new_temp_beach[xi2] = np.min(zb_all)
+
+                s['zb'][iy,xi]= new_temp_beach
+
+        if p['method_wet_supply'] == 'constant_SCR_variable_tanB':
+            beach_inc = p['shoreline_change_rate']*math.cos((math.pi/2)-math.atan(p['beach_slope']))
+            vrate = (beach_inc/(365.25*24*3600))*p['dt'] #(m/timestep)
+            hrate = (p['shoreline_change_rate']/(365.25*24*3600))*p['dt'] #(m/timestep)
+            ny, nx = s['zb'].shape        
+
+            for iy in range(ny):
+                x_all = s['x'][iy,:]
+                zb_all = s['zb'][iy,:]
+
+                xi = zb_all <= p['dune_toe_elevation']
+                beach_z = zb_all[xi]
+
+                x = x_all[xi]
+
+                xi3 = np.where(beach_z > p['zshoreline'])
+                xi3 = xi3[0][0]
+                beach_x = x-x[xi3]
+
+                xy1 = ((np.min(x[xi3])),np.min(beach_z[xi3]))
+                xy2 = (np.max(x), np.max(beach_z))
+
+                new_slope = (xy2[1]-xy1[1])/(xy2[0]-(xy1[0]))
+                b = np.min(beach_z[xi3]) + vrate
+                new_temp_beach = new_slope*(beach_x) + b
+
+                xi2 = new_temp_beach <= np.min(zb_all)
+                new_temp_beach[xi2] = np.min(zb_all)
+
+                xi4 = new_temp_beach > p['dune_toe_elevation']
+                new_temp_beach[xi4] = p['dune_toe_elevation']
+                s['zb'][iy,xi]= new_temp_beach
+    return s
 
 
 def update(s, p):
+
     '''Update bathymetry and bed composition
 
     Update bed composition by moving sediment fractions between bed
@@ -346,7 +450,7 @@ def update(s, p):
         s['zb'] += dz
         if p['process_tide']:
             s['zs'] += dz #???
-    
+
     return s
 
 
@@ -496,7 +600,6 @@ def average_change(l, s, p):
     if p['_time'] < p['avg_time']:
         s['dzbveg'] *= 0.
 
-
     return s
 
 @njit
@@ -535,5 +638,5 @@ def arrange_layers(m,dm,d,nl,ix_ero,ix_dep):
     m[ix_dep,-1,:] -= dm[ix_dep,:] * d[ix_dep,-1,:]
 
     return m
-    
+
 

aeolis/constants.py

@@ -174,6 +174,7 @@
     'process_moist'                 : False,              # Enable the process of moist
     'process_mixtoplayer'           : False,              # Enable the process of mixing 
     'process_wet_bed_reset'         : False,              # Enable the process of bed-reset in the intertidal zone
+    'process_wet_supply'            : False,              # Enable the process of beach sediment supply
     'process_meteo'                 : False,              # Enable the process of meteo
     'process_salt'                  : False,              # Enable the process of salt
     'process_humidity'              : False,              # Enable the process of humidity
@@ -340,6 +341,10 @@
     'rhoveg_max'                    : 0.5,                #maximum vegetation density, only used in duran and moore 14 formulation
     't_veg'                         : 3,                  #time scale of vegetation growth (days), only used in duran and moore 14 formulation
     'v_gam'                         : 1,                  # only used in duran and moore 14 formulation
+    'method_wet_supply'             :'wet_bed_reset',       # Name of method to comput sediment supply (wet_bed_reset, vertical_beach_growth, constant_SCR_constant_tanB, constant_SCR_variable_tanB)
+    'shoreline_change_rate'         : 0,                  # Elevation added to beach during process sed supply (shoreline change rate m/year)
+    'zshoreline'                    : 0,                  # Elevation where beach profile ends and shoreline change rate is applied           
+    'xshoreline'                    : 0,                  # Cross-shore position where beach profile ends and shoreline change rate is applied  }
 }
 
 REQUIRED_CONFIG = ['nx', 'ny']

aeolis/examples/README.md

@@ -69,7 +69,9 @@ vanWesten2024/parabolic
 <description>
 
 
+longterm_dune_growth - 1D 27-year simulation from Long Beach Penninsula, Washington, USA (LBP)
 
+The beach-dune system on LBP has been rapidly prograding over the past three decades. This example uses the newly process_sediment_supply function in AeoLiS to hindcast almost 3 decades of dune evolution.