commit cc1ced89147308322393b4220845eabe5970e39a
parent 97a5cad4554cc24b354750d75d98d7e661bab5ca
Author: Anders Damsgaard <andersd@riseup.net>
Date: Wed, 1 Mar 2017 20:43:10 -0800
begin implementing Wilcock's two-fraction sediment transport model
Diffstat:
| M | 1d-channel.py | | | 167 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------- |
1 file changed, 125 insertions(+), 42 deletions(-)
diff --git a/1d-channel.py b/1d-channel.py
@@ -22,13 +22,16 @@ import sys
# # Model parameters
Ns = 25 # Number of nodes [-]
-Ls = 100e3 # Model length [m]
-t_end = 24.*60.*60.*30. # Total simulation time [s]
+#Ls = 100e3 # Model length [m]
+Ls = 1e3 # Model length [m]
+total_days = 60. # Total simulation time [d]
+t_end = 24.*60.*60.*total_days # Total simulation time [s]
tol_Q = 1e-3 # Tolerance criteria for the normalized max. residual for Q
tol_P_c = 1e-3 # Tolerance criteria for the normalized max residual for P_c
max_iter = 1e2*Ns # Maximum number of solver iterations before failure
output_convergence = False # Display convergence statistics during run
safety = 0.1 # Safety factor ]0;1] for adaptive timestepping
+plot_interval = 20 # Time steps between plots
# Physical parameters
rho_w = 1000. # Water density [kg/m^3]
@@ -36,18 +39,23 @@ rho_i = 910. # Ice density [kg/m^3]
rho_s = 2600. # Sediment density [kg/m^3]
g = 9.8 # Gravitational acceleration [m/s^2]
theta = 30. # Angle of internal friction in sediment [deg]
+sand_fraction = 0.5 # Initial volumetric fraction of sand relative to gravel
+D_g = 1. # Mean grain size in gravel fraction (> 2 mm)
+D_s = 0.01 # Mean grain size in sand fraction (<= 2 mm)
# Water source term [m/s]
-# m_dot = 7.93e-11
-m_dot = 4.5e-7
-# m_dot = 4.5e-6
-# m_dot = 5.79e-5
+m_dot = 7.93e-11
+#m_dot = 1.0e-7
+#m_dot = 2.0e-6
+#m_dot = 4.5e-7
+#m_dot = 5.79e-5
+#m_dot = 5.0e-6
# m_dot = 1.8/(1000.*365.*24.*60.*60.) # Whillan's melt rate from Joughin 2004
# Walder and Fowler 1994 sediment transport parameters
K_e = 6.0 # Erosion constant [-], disabled when 0.0
# K_d = 6.0 # Deposition constant [-], disabled when 0.0
-K_d = K_e # Deposition constant [-], disabled when 0.0
+K_d = 0.01*K_e # Deposition constant [-], disabled when 0.0
alpha = 1e5 # Geometric correction factor (Carter et al 2017)
# D50 = 1e-3 # Median grain size [m]
# tau_c = 0.5*g*(rho_s - rho_i)*D50 # Critical shear stress for transport
@@ -67,7 +75,9 @@ c_1 = -0.118 # [m/kPa]
c_2 = 4.60 # [m]
# Minimum channel size [m^2], must be bigger than 0
-S_min = 1e-1
+#S_min = 1e-1
+#S_min = 1e-2
+S_min = 1.
# # Initialize model arrays
@@ -76,9 +86,9 @@ s = numpy.linspace(0., Ls, Ns)
ds = s[1:] - s[:-1]
# Ice thickness and bed topography
-# H = 6.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0 # glacier
-slope = 0.1 # Surface slope [%]
-H = 1000. + -slope/100.*s
+H = 6.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0 # glacier
+#slope = 0.1 # Surface slope [%]
+#H = 1000. + -slope/100.*s
b = numpy.zeros_like(H)
@@ -101,13 +111,17 @@ c_bar = numpy.zeros_like(S) # Vertically integrated sediment concentration [-]
tau = numpy.zeros_like(S) # Avg. shear stress from current [Pa]
porosity = numpy.ones_like(S)*0.3 # Sediment porosity [-]
res = numpy.zeros_like(S) # Solution residual during solver iterations
+Q_t = numpy.zeros_like(S) # Sediment flux where D <= 2 mm [m3/s]
+Q_s = numpy.zeros_like(S) # Sediment flux where D <= 2 mm [m3/s]
+Q_g = numpy.zeros_like(S) # Sediment flux where D > 2 mm [m3/s]
+f_s = numpy.ones_like(S)*sand_fraction # Initial sediment fraction of sand [-]
# # Helper functions
def gradient(arr, arr_x):
# Central difference gradient of an array ``arr`` with node positions at
# ``arr_x``.
- return (arr[:-1] - arr[1:])/(arr_x[:-1] - arr_x[1:])
+ return (arr[1:] - arr[:-1])/(arr_x[1:] - arr_x[:-1])
def avg_midpoint(arr):
@@ -129,27 +143,21 @@ def channel_shear_stress(Q, S):
def channel_erosion_rate(tau):
# Parker 1979, Walder and Fowler 1994
# return K_e*v_s*(tau - tau_c).clip(min=0.)/(g*(rho_s - rho_w)*d15)**(3./2)
+
# Carter et al 2017
- return K_e*v_s/alpha*(tau - tau_c).clip(min=0.) / \
- (g*(rho_s - rho_w)*d15)**(3./2.)
+ #return K_e*v_s/alpha*(tau - tau_c).clip(min=0.) / \
+ #(g*(rho_s - rho_w)*d15)**(3./2.)
+
+ # Ng 2000
+ return 0.092*(tau/(2.*(rho_s - rho_w)*g*d15))**(3./2.)
def channel_deposition_rate_kernel(tau, c_bar, ix):
# Parker 1979, Walder and Fowler 1994
- # result = K_d*v_s*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
+ # return K_d*v_s*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
# Carter et al. 2017
- result = K_d*v_s/alpha*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
-
- '''
- print('tau[{}] = {}'.format(ix, tau[ix]))
- print('c_bar[{}] = {}'.format(ix, c_bar[ix]))
- print('e_dot[{}] = {}'.format(ix, e_dot[ix]))
- print('d_dot[{}] = {}'.format(ix, result))
- print('')
- '''
-
- return result
+ return K_d*v_s/alpha*c_bar[ix]*(g*(rho_s - rho_w)*d15/tau[ix])**0.5
def channel_deposition_rate_kernel_ng(c_bar, ix):
@@ -196,22 +204,81 @@ def channel_deposition_rate(tau, c_bar, d_dot, Ns):
return d_dot, c_bar
-def channel_growth_rate(e_dot, d_dot, porosity, W):
+def channel_sediment_flux_sand(tau, W, f_s, D_s):
+ # Parker 1979, Wilcock 1997, 2001, Egholm 2013
+ # tau: Shear stress by water flow
+ # W: Channel width
+ # f_s: Sand volume fraction
+ # D_s: Mean sand fraction grain size
+
+ # Piecewise linear functions for nondimensional critical shear stresses
+ # dependent on sand fraction from Gasparini et al 1999 of Wilcock 1997
+ # data.
+ ref_shear_stress = numpy.ones_like(f_s)*0.04
+ ref_shear_stress[numpy.nonzero(f_s <= 0.1)] = 0.88
+ I = numpy.nonzero((0.1 < f_s) & (f_s <= 0.4))
+ ref_shear_stress[I] = 0.88 - 2.8*(f_s[I] - 0.1)
+
+ shields_stress = tau/((rho_s - rho_w)*g*D_s)
+
+ #import ipdb; ipdb.set_trace()
+ return 11.2*f_s*W/((rho_s - rho_w)/rho_w*g) \
+ * (tau/rho_w)**1.5 \
+ * (1.0 - 0.846*numpy.sqrt(ref_shear_stress/shields_stress))**4.5
+
+def channel_sediment_flux_gravel(tau, W, f_g, D_g):
+ # Parker 1979, Wilcock 1997, 2001, Egholm 2013
+ # tau: Shear stress by water flow
+ # W: Channel width
+ # f_g: Gravel volume fraction
+ # D_g: Mean gravel fraction grain size
+
+ # Piecewise linear functions for nondimensional critical shear stresses
+ # dependent on sand fraction from Gasparini et al 1999 of Wilcock 1997
+ # data.
+ ref_shear_stress = numpy.ones_like(f_g)*0.01
+ ref_shear_stress[numpy.nonzero(f_g <= 0.1)] = 0.04
+ I = numpy.nonzero((0.1 < f_g) & (f_g <= 0.4))
+ ref_shear_stress[I] = 0.04 - 0.1*(f_g[I] - 0.1)
+
+ shields_stress = tau/((rho_s - rho_w)*g*D_g)
+
+ # From Wilcock 2001, eq. 3
+ Q_g = 11.2*f_g*W/((rho_s - rho_w)/rho_w*g) \
+ * (tau/rho_w)**1.5 \
+ * (1.0 - 0.846*ref_shear_stress/shields_stress)**4.5
+
+ # From Wilcock 2001, eq. 4
+ I = numpy.nonzero(ref_shear_stress/shields_stress < 1.)
+ Q_g[I] = f_g[I]*W[I]/((rho_s - rho_w)/rho_w*g) \
+ * (tau[I]/rho_w)**1.5 \
+ * 0.0025*(shields_stress[I]/ref_shear_stress[I])**14.2
+
+ return Q_g
+
+
+def channel_growth_rate(e_dot, d_dot, W):
# Damsgaard et al, in prep
return (e_dot - d_dot)*W
-def update_channel_size_with_limit(S, S_old, dSdt, dt, N):
+def channel_growth_rate_sedflux(Q_t, porosity, s_c):
+ # Damsgaard et al, in prep
+ return 1./porosity[1:] * gradient(Q_t, s_c)
+
+
+def update_channel_size_with_limit(S, S_old, dSdt, dt, N_c):
# Damsgaard et al, in prep
- S_max = numpy.maximum((c_1*numpy.maximum(N, 0.)/1000. + c_2)**2./4. *
- numpy.tan(numpy.deg2rad(theta)), S_min)
+ S_max = numpy.maximum(
+ numpy.maximum((c_1*numpy.maximum(N_c, 0.)/1000. + c_2), 0.)**2./4. *
+ numpy.tan(numpy.deg2rad(theta)), S_min)
S = numpy.maximum(numpy.minimum(S_old + dSdt*dt, S_max), S_min)
W = S/numpy.tan(numpy.deg2rad(theta)) # Assume no channel floor wedge
return S, W, S_max
def flux_solver(m_dot, ds):
- # Iteratively find new fluxes
+ # Iteratively find new water fluxes
it = 0
max_res = 1e9 # arbitrary large value
@@ -286,16 +353,17 @@ def plot_state(step, time, S_, S_max_, title=True):
ax_Pa = plt.subplot(2, 1, 1) # axis with Pascals as y-axis unit
#ax_Pa.plot(s_c/1000., P_c/1000., '--r', label='$P_c$')
#ax_Pa.plot(s/1000., N/1000., '--r', label='$N$')
- ax_Pa.plot(s/1000., H*rho_i*g/1e6, '--r', label='$P_i$')
- ax_Pa.plot(s_c/1000., P_c/1e6, ':b', label='$P_c$')
+ ax_Pa.plot(s_c/1000., N_c/1e6, '-k', label='$N$')
+ ax_Pa.plot(s_c/1000., H_c*rho_i*g/1e6, '--r', label='$P_i$')
+ ax_Pa.plot(s_c/1000., P_c/1e6, ':y', label='$P_c$')
ax_m3s = ax_Pa.twinx() # axis with m3/s as y-axis unit
- ax_m3s.plot(s_c/1000., Q, '-k', label='$Q$')
+ ax_m3s.plot(s_c/1000., Q, '.-b', label='$Q$')
if title:
plt.title('Day: {:.3}'.format(time/(60.*60.*24.)))
ax_Pa.legend(loc=2)
- ax_m3s.legend(loc=3)
+ ax_m3s.legend(loc=4)
#ax_Pa.set_ylabel('[kPa]')
ax_Pa.set_ylabel('[MPa]')
ax_m3s.set_ylabel('[m$^3$/s]')
@@ -307,8 +375,12 @@ def plot_state(step, time, S_, S_max_, title=True):
# ax_m.semilogy(s_c/1000., S_max, '--k', label='$S_{max}$')
ax_ms = ax_m2.twinx()
- ax_ms.plot(s_c/1000., e_dot, '--r', label='$\dot{e}$')
- ax_ms.plot(s_c/1000., d_dot, ':b', label='$\dot{d}$')
+ #ax_ms.plot(s_c/1000., e_dot, '--r', label='$\dot{e}$')
+ #ax_ms.plot(s_c/1000., d_dot, ':b', label='$\dot{d}$')
+ ax_ms.plot(s_c/1000., Q_t, label='$Q_t$')
+ ax_ms.plot(s_c/1000., Q_g, label='$Q_g$')
+ ax_ms.plot(s_c/1000., Q_s, label='$Q_s$')
+ # TODO: check units on sediment fluxes: m2/s or m3/s ?
ax_m2.legend(loc=2)
ax_ms.legend(loc=3)
@@ -351,7 +423,7 @@ hydro_pot_grad = gradient(hydro_pot, s)
# Find field values at the middle of channel segments
N_c = avg_midpoint(N)
-H_c = avg_midpoint(N)
+H_c = avg_midpoint(H)
# Find fluxes in channel segments [m^3/s]
Q = channel_water_flux(S, hydro_pot_grad)
@@ -389,11 +461,21 @@ while time <= t_end:
tau = channel_shear_stress(Q, S)
# Find sediment erosion and deposition rates for each channel segment
- e_dot = channel_erosion_rate(tau)
- d_dot, c_bar = channel_deposition_rate(tau, c_bar, d_dot, Ns)
+ # e_dot = channel_erosion_rate(tau)
+ # d_dot, c_bar = channel_deposition_rate(tau, c_bar, d_dot, Ns)
+
+ f_g = 1./f_s # gravel volume fraction is reciprocal to sand
+ Q_s = channel_sediment_flux_sand(tau, W, f_s, D_s)
+ Q_g = channel_sediment_flux_gravel(tau, W, f_g, D_g)
+ Q_t = Q_s + Q_g
+ break
+
+ # TODO: Update f_s from fluxes
# Determine change in channel size for each channel segment
- dSdt = channel_growth_rate(e_dot, d_dot, porosity, W)
+ #dSdt = channel_growth_rate(e_dot, d_dot, W)
+ # Use backward differences and assume dS/dt=0 in first segment
+ dSdt[1:] = channel_growth_rate_sedflux(Q_t, porosity, s_c)
# Update channel cross-sectional area and width according to growth rate
# and size limit for each channel segment
@@ -424,7 +506,8 @@ while time <= t_end:
'Iterative solution not found for Q')
it += 1
- if step % 10 == 0:
+ # Generate an output figure for every n time steps
+ if step % plot_interval == 0:
plot_state(step, time, S, S_max)
# import ipdb; ipdb.set_trace()
