From 60f8a3474919fde82663ade00b74d9bc7e205c3b Mon Sep 17 00:00:00 2001 From: Christa van IJzendoorn <52737913+christavanijzendoorn@users.noreply.github.com> Date: Wed, 19 Feb 2025 17:44:35 -0800 Subject: [PATCH 1/3] Dynamic revetment implementation - New runup method included - z0 can be defined based on grid - shear module can handle spatially varying z0 - implementation of sand infilling within cobble - sediment trapping implemented based on roughness density - constants added where needed --- aeolis/bed.py | 187 +++++++++++++++++++++++++++++++++++++++++++- aeolis/constants.py | 16 +++- aeolis/hydro.py | 70 +++++++++++------ aeolis/shear.py | 21 +++-- aeolis/threshold.py | 99 ++++++++++++++++++++++- aeolis/wind.py | 9 ++- 6 files changed, 360 insertions(+), 42 deletions(-) diff --git a/aeolis/bed.py b/aeolis/bed.py index a71214d6..e3c371fd 100644 --- a/aeolis/bed.py +++ b/aeolis/bed.py @@ -114,7 +114,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 @@ -314,8 +323,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 #??? @@ -504,5 +517,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: + print('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 diff --git a/aeolis/constants.py b/aeolis/constants.py index a1c0c33b..cdcaa0ef 100644 --- a/aeolis/constants.py +++ b/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', #NEWCH # [-] Moisture content (volumetric) 'moist_swr', #NEWCH # [-] Moisture content soil water retention relationship (volumetric) @@ -154,6 +158,7 @@ '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 @@ -161,7 +166,8 @@ '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 @@ -195,6 +201,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 'wave_mask' : None, # Filename of ASCII file with mask for wave height 'tide_mask' : None, # Filename of ASCII file with mask for tidal elevation 'runup_mask' : None, # Filename of ASCII file with mask for run-up @@ -246,7 +256,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 @@ -283,6 +293,7 @@ 'Cl_gw' : 0.7, # NEWCH # [m] Groundwater overheight due to runup 'in_gw' : 0, # NEWCH # [m] Initial groundwater level 'GW_stat' : 1, # NEWCH # [m] Landward static groundwater boundary (if static boundary is defined) + 'max_moist' : 10., # NEWCH # [%] Moisture content (volumetric in percent) above which the threshold shear velocity is set to infinity (no transport, default value Delgado-Fernandez, 2010) 'theta_dyn' : 33., # [degrees] Initial Dynamic angle of repose, critical dynamic slope for avalanching 'theta_stat' : 34., # [degrees] Initial Static angle of repose, critical static slope for avalanching 'avg_time' : 86400., # [s] Indication of the time period over which the bed level change is averaged for vegetation growth @@ -321,6 +332,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 diff --git a/aeolis/hydro.py b/aeolis/hydro.py index 620e6296..17096cd8 100644 --- a/aeolis/hydro.py +++ b/aeolis/hydro.py @@ -124,6 +124,32 @@ def interpolate(s, p, t): p['wave_file'][:,0], p['wave_file'][:,2]) + if p['process_runup']: + ny = p['ny'] + + for iy in range(ny + 1): # do this computation seperately on every y for now so alongshore variable wave runup can be added in the future + + hs = s['Hs'][iy][0] + tp = s['Tp'][iy][0] + wl = s['SWL'][iy][0] + + 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 + + if hasattr(s['runup_mask'], "__len__"): + s['eta'][iy][:] = apply_mask(s['eta'][iy][:], s['runup_mask'][iy][:]) + s['R'][iy][:] = apply_mask(s['R'][iy][:], s['runup_mask'][iy][:]) + + s['TWL'][iy][:] = s['SWL'][iy][:] + s['R'][iy][:] + s['DSWL'][iy][:] = s['SWL'][iy][:] + s['eta'][iy][:] # Was s['zs'] before + + # Alters wave height based on maximum wave height over depth ratio, gamma default = 0.5 s['Hs'] = np.minimum(h * p['gamma'], s['Hs']) # apply complex mask @@ -140,32 +166,15 @@ def interpolate(s, p, t): s['Tp'] = apply_mask(s['Tp'], s['wave_mask']) + if p['process_runup']: ny = p['ny'] - - if ('Hs' not in p['external_vars']): - - for iy in range(ny + 1): # do this computation seperately on every y for now so alongshore variable wave runup can be added in the future - - hs = s['Hs'][iy][0] - tp = s['Tp'][iy][0] - wl = s['SWL'][iy][0] - - eta, sigma_s, R = calc_runup_stockdon(hs, tp, p['beach_slope']) - s['R'][iy][:] = R - s['eta'][iy][:] = eta - s['sigma_s'][iy][:] = sigma_s - - if hasattr(s['runup_mask'], "__len__"): - s['eta'][iy][:] = apply_mask(s['eta'][iy][:], s['runup_mask'][iy][:]) - s['R'][iy][:] = apply_mask(s['R'][iy][:], s['runup_mask'][iy][:]) - - s['TWL'][iy][:] = s['SWL'][iy][:] + s['R'][iy][:] - s['DSWL'][iy][:] = s['SWL'][iy][:] + s['eta'][iy][:] # Was s['zs'] before - 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__"): @@ -670,6 +679,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 + diff --git a/aeolis/shear.py b/aeolis/shear.py index dd33ee52..e0763dd3 100644 --- a/aeolis/shear.py +++ b/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,7 @@ 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 diff --git a/aeolis/threshold.py b/aeolis/threshold.py index db8f688d..caecf4e3 100644 --- a/aeolis/threshold.py +++ b/aeolis/threshold.py @@ -94,9 +94,11 @@ def compute(s, p): s = compute_humidity(s, p) if p['th_salt']: s = compute_salt(s, p) - if p['th_sheltering']: + if p['th_sheltering']: # sheltering based on grain size dsitribution s = compute_sheltering(s, p) - + if p['th_sedtrapping']: + s = compute_sedtrapping(s, p) + # apply complex mask s['uth'] = apply_mask(s['uth'], s['threshold_mask']) s['uthf'] = s['uth'].copy() @@ -231,9 +233,11 @@ def compute_moisture(s, p): else: logger.log_and_raise('Unknown moisture formulation [%s]' % p['method_moist'], exc=ValueError) + th_mg = (p['max_moist']/100 * p['rhow'] / (p['rhog'] * (1. - p['porosity']))) # should be .04 according to Pye and Tsoar - # should be .64 according to Delgado-Fernandez (10% vol.) - ix = mg > 0.064 + # should be .064 according to Delgado-Fernandez (10% vol.) + # .097 corresponds to 15% + ix = mg > th_mg s['uth'][ix] = np.inf return s @@ -394,6 +398,93 @@ def compute_sheltering(s, p): return s +def compute_sedtrapping(s, p): + '''Modify wind velocity threshold based on the presence of roughness + elements following Raupach (1993) + + The formulation used in compute_sheltering depends on the grain size + distribution to account for the emergence of roughness elements due + to grain size sorting. However, that formulation is not applicable to + larger scale roughness elements such as cobbles on a dynamic + revetment. + + Here, an altered version of the Raupach (1993) formulation is used + that alters the threshold shear velocity based on a pre-defined + value of lambda. + + :math:`\\lambda` is the roughness density defined as the number of + elements per surface area :math:`\\frac{n}{S}` multiplied by the + frontal area of a roughness element :math:`A_f`, also known as the + frontal area index: + + .. math:: + + \\lambda = \\frac{n b h}{S} = \\frac{n A_f}{S} + + Here, it is a assumed that the values for lambda can be calculated + based on the characteristics of the roughness elements. Lambda can + be assigned to the entire domain by providing a single value. + Alternatively, a grid with lambda values can be provided based on + the location of roughness elements. + + .. math:: + + R_t = \\frac{u_{*,th,s}}{u_{*,th,r}} + = \\sqrt{\\frac{\\tau_s''}{\\tau}} + = \\frac{1}{\\sqrt{\\left( 1 - m \\sigma \\lambda \\right) + \\left( 1 + m \\beta \\lambda \\right)}} + + :math:`m` is a constant smaller or equal to unity that accounts + for the difference between the average stress on the bed surface + :math:`\\tau_s` and the maximum stress on the bed surface + :math:`\\tau_s''`. Here, it is assumed 0.5 for field applications. + + :math:`\\sigma` is the shape coefficient defined as the basal area + divided by the frontal area: :math:`\\frac{A_b}{A_f}`. For + hemispheres :math:`\\sigma = 2`, for spheres :math:`\\sigma = 1`. + :math:`\\sigma = 2` is assumed for a cobble surface. + + :math:`\\beta` is the stress partition coefficient defined as the + ratio between the drag coefficient of the roughness element itself + :math:`C_r` and the drag coefficient of the bare surface without + roughness elements :math:`C_s`. Calibration parameter, assumed 130 + by McKenna and Neumann (2012) for a shell surface, same value used + by Hoonhout and de Vries (2016) to simulate sheltering. + + Parameters + ---------- + s : dict + Spatial grids + p : dict + Model configuration parameters + + Returns + ------- + dict + Spatial grids + + ''' + + nx = p['nx']+1 + ny = p['ny']+1 + + # compute inverse of shear stress ratio based on lambda + Rti = np.sqrt((1. - p['m'] * p['sigma'] * p['lambda']) * (1. + p['m'] * p['beta'] * p['lambda'])) + s['Rti'] = Rti + + # modify shear velocity threshold + if isinstance(p['lambda'], int) or isinstance(p['lambda'], float): + s['uth'] *= Rti + else: + s['uth'] *= Rti.reshape((ny,nx,1)) + + # where more than 10 cm pure cobble is present at the surface, set uth to inf + if p['dcob_file'] is not None and p['doverlap_file'] is not None: + ix = s['dcob']-s['doverlap'] >= 0.1 + s['uth'][ix] = np.inf + + return s + def non_erodible(s,p): '''Modify wind velocity threshold based on the presence of a non-erodible layer. diff --git a/aeolis/wind.py b/aeolis/wind.py index c4f51d9c..36027aa3 100644 --- a/aeolis/wind.py +++ b/aeolis/wind.py @@ -70,7 +70,7 @@ def initialize(s, p): # initialize wind shear model (z0 according to Duran much smaller) # Otherwise no Barchan z0 = calculate_z0(p, s) - + if p['process_shear']: if p['ny'] > 0: s['shear'] = aeolis.shear.WindShear(s['x'], s['y'], s['zb'], @@ -174,12 +174,15 @@ def calculate_z0(p, s): z0 ''' + if p['method_roughness'] == 'constant': z0 = p['k'] # Here, the ks (roughness length) is equal to the z0, this method is implemented to assure backward compatibility. Note, this does not follow the definition of z0 = ks /30 by Nikuradse if p['method_roughness'] == 'constant_nikuradse': z0 = p['k'] / 30 # This equaion follows the definition of the bed roughness as introduced by Nikuradse if p['method_roughness'] == 'mean_grainsize_initial': #(based on Nikuradse and Bagnold, 1941), can only be applied in case with uniform grain size and is most applicable to a flat bed z0 = np.sum(p['grain_size']*p['grain_dist']) / 30. + if p['method_roughness'] == 'z0_grid' and isarray(p['z0_file']): + z0 = p['z0_file'][:,:] if p['method_roughness'] == 'mean_grainsize_adaptive': # makes Nikuradse roughness method variable through time and space depending on grain size variations z0 = calc_mean_grain_size(p, s) / 30. if p['method_roughness'] == 'median_grainsize_adaptive': # based on Sherman and Greenwood, 1982 - only appropriate for naturally occurring grain size distribution @@ -212,14 +215,14 @@ def shear(s,p): # Compute shear velocity field (including separation) if 'shear' in s.keys() and p['process_shear'] and p['ny'] > 0: - + s['shear'](x=s['x'], y=s['y'], z=s['zb'], taux=s['taus'], tauy=s['taun'], u0=s['uw'][0,0], udir=s['udir'][0,0], process_separation = p['process_separation'], c = p['c_b'], mu_b = p['mu_b'], - taus0 = s['taus0'][0,0], taun0 = s['taun0'][0,0], + taus0 = s['taus0'], taun0 = s['taun0'], sep_filter_iterations=p['sep_filter_iterations'], zsep_y_filter=p['zsep_y_filter']) From 8bee6aeacab09db775c0b1e21b11f12dec1e1de6 Mon Sep 17 00:00:00 2001 From: Christa van IJzendoorn <52737913+christavanijzendoorn@users.noreply.github.com> Date: Mon, 21 Jul 2025 15:54:42 -0700 Subject: [PATCH 2/3] Update aeolis/threshold.py spelling corrected Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- aeolis/threshold.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/aeolis/threshold.py b/aeolis/threshold.py index deec7f57..b347ab91 100644 --- a/aeolis/threshold.py +++ b/aeolis/threshold.py @@ -94,7 +94,7 @@ def compute(s, p): s = compute_humidity(s, p) if p['th_salt']: s = compute_salt(s, p) - if p['th_sheltering']: # sheltering based on grain size dsitribution + if p['th_sheltering']: # sheltering based on grain size distribution s = compute_sheltering(s, p) if p['th_sedtrapping']: s = compute_sedtrapping(s, p) From f3dc23de3d787745169c38cfb3556cc8a8602b1e Mon Sep 17 00:00:00 2001 From: Christa van IJzendoorn <52737913+christavanijzendoorn@users.noreply.github.com> Date: Mon, 21 Jul 2025 15:57:25 -0700 Subject: [PATCH 3/3] Update aeolis/bed.py proper warning Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- aeolis/bed.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/aeolis/bed.py b/aeolis/bed.py index 6190b942..093fa5eb 100644 --- a/aeolis/bed.py +++ b/aeolis/bed.py @@ -629,7 +629,7 @@ def update_composite(s, p): 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: - print('Cobble thickness is smaller than overlap thickness') + 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]