Skip to content

Long term dune growth (by Selwyn) #286

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Sierd wants to merge 5 commits into
base: main
Choose a base branch
from
+677,883 −21
Open

Long term dune growth (by Selwyn) #286

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
b617b80
Clean commit of Sewlyns code
Sierd Feb 19, 2026
1f469f7
test wetbedreset
Sierd Feb 19, 2026
208d045
Rename sediment_supply to wet_supply and update related parameters fo…
ssheminway Feb 20, 2026
5db83bb
updated constants.py and vegetation.py
Sierd Feb 23, 2026
eb5d66b
Merge branch 'sed_supply_clean' of https://github.com/openearth/aeoli…
Sierd Feb 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion aeolis/avalanching.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
137 changes: 120 additions & 17 deletions aeolis/bed.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 *
Expand Down Expand Up @@ -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


Expand Down Expand Up @@ -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


Comment on lines +212 to 232
Copy link

Copilot AI Feb 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Large block of commented-out code (lines 204-230) should be removed rather than kept in the codebase. If this code needs to be preserved for reference, it belongs in version control history. Keeping extensive commented code reduces readability and can cause confusion about what is actually being used.

Suggested change
# 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

Copilot uses AI. Check for mistakes.
def wet_supply(s, p):
''' Increase elevation of beach topography.

Parameters
----------
Expand All @@ -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']:
Copy link

Copilot AI Feb 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The naming inconsistency between process flags is confusing. The constants file defines process_wet_supply (line 177), but the configuration file uses process_sediment_supply (line 49 in aeolis.txt). The function checks both p['process_wet_supply'] and p['process_wet_bed_reset'], but the configuration doesn't define process_wet_supply. This inconsistency will likely cause the feature not to work as expected. Use consistent naming throughout.

Copilot uses AI. Check for mistakes.

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)
Copy link

Copilot AI Feb 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using * for multiplication with time conversion in line 262 is correct, but the similar expression in line 261 uses division. The inconsistent notation pattern between (1/365.25/24/3600) and (beach_inc/(365.25*24*3600)) is confusing. Use consistent mathematical notation for time conversions throughout the function.

Suggested change
vrate = (beach_inc*(1/365.25/24/3600))*p['dt'] #(m/timestep)
vrate = (beach_inc/(365.25*24*3600))*p['dt'] # (m/timestep)

Copilot uses AI. Check for mistakes.
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)
Copy link

Copilot AI Feb 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Division by multiple time factors appears unnecessarily complex. The expression (beach_inc/(365.25*24*3600))*p['dt'] could be simplified for clarity. Consider computing the conversion factor once and documenting the units more clearly in comments.

Copilot uses AI. Check for mistakes.
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'])
Copy link

Copilot AI Feb 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Potential issue with accessing xi3[0][0] without first checking if xi3[0] is non-empty. If np.where(beach_z > p['zshoreline']) returns an empty array, accessing xi3[0][0] will raise an IndexError. Add a check to ensure the array is not empty before indexing, or handle the case where no points exceed the shoreline elevation.

Suggested change
xi3 = np.where(beach_z > p['zshoreline'])
xi3 = np.where(beach_z > p['zshoreline'])
# Guard against empty result: if no points exceed zshoreline, skip this profile
if xi3[0].size == 0:
continue

Copilot uses AI. Check for mistakes.
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':
Comment on lines +250 to +302
Copy link

Copilot AI Feb 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing input validation for the method_wet_supply parameter. If an invalid method name is provided in the configuration, the function will silently do nothing without any error or warning. Add validation to check if the method is one of the supported values and raise an informative error if not.

Copilot uses AI. Check for mistakes.
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']
Copy link

Copilot AI Feb 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The code uses both <= (line 312) and < (line 268, 287) when comparing to p['dune_toe_elevation']. This inconsistency could lead to boundary cells being handled differently in different methods. Verify whether the dune toe elevation should be included or excluded, and use consistent comparison operators across all methods.

Suggested change
xi = zb_all <= p['dune_toe_elevation']
xi = zb_all < p['dune_toe_elevation']

Copilot uses AI. Check for mistakes.
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]))
Comment on lines +317 to +321
Copy link

Copilot AI Feb 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same indexing issue as above: accessing xi3[0][0] on line 318 without verifying the array is non-empty could raise an IndexError. This pattern appears in multiple places (lines 292, 318) and should be handled consistently with proper validation.

Suggested change
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]))
xi3 = np.where(beach_z > p['zshoreline'])[0]
if xi3.size == 0:
# No points satisfy beach_z > zshoreline for this transect; skip update
continue
xi3 = xi3[0]
beach_x = x - x[xi3]
xy1 = ((np.min(x[xi3])), np.min(beach_z[xi3]))

Copilot uses AI. Check for mistakes.
xy2 = (np.max(x), np.max(beach_z))

new_slope = (xy2[1]-xy1[1])/(xy2[0]-(xy1[0]))
Comment on lines +321 to +324
Copy link

Copilot AI Feb 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The variable xy1 is assigned but never used (line 321). Similarly, xy2 appears to be defined for calculating slope but the actual slope calculation on line 324 doesn't use these tuples correctly. The slope calculation uses individual array elements when the tuple structure suggests coordinate pairs were intended. Either remove these unused variables or fix the slope calculation to properly use them.

Suggested change
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]))
# compute slope directly between (x[xi3], beach_z[xi3]) and (max(x), max(beach_z))
new_slope = (np.max(beach_z) - np.min(beach_z[xi3])) / (np.max(x) - np.min(x[xi3]))

Copilot uses AI. Check for mistakes.
b = np.min(beach_z[xi3]) + vrate
new_temp_beach = new_slope*(beach_x) + b
Comment on lines +325 to +326
Copy link

Copilot AI Feb 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The variable name b on lines 270, 293, and 325 is not descriptive. This appears to represent a y-intercept or baseline elevation value. Use a more descriptive name like baseline_elevation or y_intercept to improve code readability.

Suggested change
b = np.min(beach_z[xi3]) + vrate
new_temp_beach = new_slope*(beach_x) + b
baseline_elevation = np.min(beach_z[xi3]) + vrate
new_temp_beach = new_slope*(beach_x) + baseline_elevation

Copilot uses AI. Check for mistakes.

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
Expand Down Expand Up @@ -346,7 +450,7 @@ def update(s, p):
s['zb'] += dz
if p['process_tide']:
s['zs'] += dz #???

return s


Expand Down Expand Up @@ -496,13 +600,12 @@ def average_change(l, s, p):
if p['_time'] < p['avg_time']:
s['dzbveg'] *= 0.


return s

@njit
def arrange_layers(m,dm,d,nl,ix_ero,ix_dep):
'''Arranges mass redistrubution between layers.
This function is called in the bed.update fucntion to speed up code using numba
'''Arranges mass redistribution between layers.
This function is called in the bed.update function to speed up code using numba



Expand All @@ -511,7 +614,7 @@ def arrange_layers(m,dm,d,nl,ix_ero,ix_dep):
m : array
mass in layers
dm : array
total mass exchanged between layers derrived from pickup
total mass exchanged between layers derived from pickup
d : array
normalized mass in layers
nl : int
Expand All @@ -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


5 changes: 5 additions & 0 deletions aeolis/constants.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Copy link

Copilot AI Feb 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The parameter name method_wet_supply in constants.py (line 344) doesn't match method_sed_supply used in the configuration file (line 195 of aeolis.txt). This inconsistency will cause the method selection to fail. Ensure the parameter names match between the constants definition and the configuration file.

Suggested change
'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)
'method_sed_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)

Copilot uses AI. Check for mistakes.
'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 }
Copy link

Copilot AI Feb 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The closing brace is missing for the dictionary in constants.py. Line 347 has a closing brace in a comment (# ... }) which suggests it was accidentally commented out. This will cause a syntax error. Uncomment or add the proper closing brace.

Suggested change
'xshoreline' : 0, # Cross-shore position 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

Copilot uses AI. Check for mistakes.
}

REQUIRED_CONFIG = ['nx', 'ny']
Expand Down
2 changes: 2 additions & 0 deletions aeolis/examples/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,9 @@ vanWesten2024/parabolic
<description>


longterm_dune_growth - 1D 27-year simulation from Long Beach Penninsula, Washington, USA (LBP)
Copy link

Copilot AI Feb 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typo in the comment: "Penninsula" should be "Peninsula".

Suggested change
longterm_dune_growth - 1D 27-year simulation from Long Beach Penninsula, Washington, USA (LBP)
longterm_dune_growth - 1D 27-year simulation from Long Beach Peninsula, Washington, USA (LBP)

Copilot uses AI. Check for mistakes.

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.



Expand Down
Loading
Loading