Add scr working

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 *
@@ -73,6 +74,12 @@ def initialize(s, p):
     s['zb0'][:,:] = p['bed_file']
     s['zne'][:,:] = p['ne_file']
 
+    #initialize variable shoreline, beach slope, dune toe elevation
+    s['zshore'][:,:] = p['zshoreline']
+    s['xshore'][:,:] = p['xshoreline']
+    s['beachslope'][:,:] = p['beach_slope']
+    s['dte'][:,:] = p['dune_toe_elevation']
+
     #initialize thickness of erodable or dry top layer
     s['zdry'][:,:] = 0.05
 
@@ -125,7 +132,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 +204,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 sediment_supply(s, p):
+    ''' Increase elevation of beach topography.
 
     Parameters
     ----------
@@ -219,19 +252,106 @@ 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_sediment_supply']:
+
+        if p['method_sed_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['method_sed_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)
+            s['zshore'] = s['zshore'] + vrate
+
+            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_sed_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]
+                s['dte'] = beach_z[-1]
 
+                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_sed_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]
+
+                s['dte'][iy,:] = beach_z[-1]
+
+                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]))
+                s['beachslope'][iy,:] = new_slope
+
+                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['xshore'][iy,:] = s['xshore'][0][0] - hrate
+                s['zshore'][iy,:] = s['zshore'][0][0] + vrate
+
+                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 +466,7 @@ def update(s, p):
         s['zb'] += dz
         if p['process_tide']:
             s['zs'] += dz #???
-    
+
     return s
 
 
@@ -496,7 +616,6 @@ def average_change(l, s, p):
     if p['_time'] < p['avg_time']:
         s['dzbveg'] *= 0.
 
-
     return s
 
 @njit
@@ -535,5 +654,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

@@ -142,6 +142,10 @@
         'unST',                            # [NEW] [m/s] Component of the saltation velocity in y-direction for SedTRAILS
         'u0',
         'masstop',                          # [kg/m^2] Sediment mass in bed toplayer, only stored for output
+        'zshore',
+        'xshore',
+        'dte',
+        'beachslope',
     ),
     ('ny','nx','nlayers') : (
         'thlyr',                            # [m] Bed composition layer thickness
@@ -174,7 +178,8 @@
     'process_runup'                 : False,              # Enable the process of wave runup
     '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_bed_reset'         : False,              # Enable the process of bed-reset in the intertidal zone
+    'process_sediment_supply'       : False,              # Enable the process of linear shoreline change rate
     'process_meteo'                 : False,              # Enable the process of meteo
     'process_salt'                  : False,              # Enable the process of salt
     'process_humidity'              : False,              # Enable the process of humidity
@@ -342,6 +347,11 @@
     '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_sed_supply'             :'pos_shoreline_change', # Name of method to comput sediment supply (pos_shoreline_change, wet_bed_reset)
+    'shoreline_change_rate'         : 0,                  # Elevation added to beach during process sed supply (shoreline change rate m/year)
+    'zshoreline'                    :0,
+    'xshoreline'                    :0
+
 }
 
 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. 
 
 
 

aeolis/model.py

@@ -353,8 +353,11 @@ def update(self, dt:float=-1) -> None:
         self.s = aeolis.avalanching.angele_of_repose(self.s, self.p)
         self.s = aeolis.avalanching.avalanche(self.s, self.p)
 
-        # reset original bed in marine zone (wet)
-        self.s = aeolis.bed.wet_bed_reset(self.s, self.p)
+        # # reset original bed in marine zone (wet)
+        # self.s = aeolis.bed.wet_bed_reset(self.s, self.p)
+
+        # either wet bed reset or positive shoreline change 
+        self.s = aeolis.bed.sediment_supply(self.s, self.p)
 
         # calculate average bed level change over time
         self.s = aeolis.bed.average_change(self.l, self.s, self.p)