granular-channel-hydro

Subglacial hydrology model for sedimentary channels
git clone git://src.adamsgaard.dk/granular-channel-hydro
Log | Files | Refs | README | LICENSE

commit 5de4d190454b1c94b1982063f442efcfc30ee649
parent cc1ced89147308322393b4220845eabe5970e39a
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Fri,  3 Mar 2017 09:45:13 -0800

fix comment style

Diffstat:
M1d-channel.py | 55+++++++++++++++++++++++++++++--------------------------
1 file changed, 29 insertions(+), 26 deletions(-)

diff --git a/1d-channel.py b/1d-channel.py @@ -22,7 +22,7 @@ import sys # # Model parameters Ns = 25 # Number of nodes [-] -#Ls = 100e3 # Model length [m] +# 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] @@ -45,11 +45,11 @@ D_s = 0.01 # Mean grain size in sand fraction (<= 2 mm) # Water source term [m/s] 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.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 @@ -75,8 +75,8 @@ 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-2 +# S_min = 1e-1 +# S_min = 1e-2 S_min = 1. @@ -87,8 +87,8 @@ 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 +# slope = 0.1 # Surface slope [%] +# H = 1000. + -slope/100.*s b = numpy.zeros_like(H) @@ -114,7 +114,7 @@ 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 [-] +f_s = numpy.ones_like(S)*sand_fraction # Initial sediment fraction of sand [-] # # Helper functions @@ -145,8 +145,8 @@ def channel_erosion_rate(tau): # 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.) @@ -219,13 +219,15 @@ def channel_sediment_flux_sand(tau, W, f_s, D_s): I = numpy.nonzero((0.1 < f_s) & (f_s <= 0.4)) ref_shear_stress[I] = 0.88 - 2.8*(f_s[I] - 0.1) + # Non-dimensionalize shear stress shields_stress = tau/((rho_s - rho_w)*g*D_s) - #import ipdb; ipdb.set_trace() + # 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 @@ -241,18 +243,19 @@ def channel_sediment_flux_gravel(tau, W, f_g, D_g): I = numpy.nonzero((0.1 < f_g) & (f_g <= 0.4)) ref_shear_stress[I] = 0.04 - 0.1*(f_g[I] - 0.1) + # Non-dimensionalize shear stress 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 + * (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 + * (tau[I]/rho_w)**1.5 \ + * 0.0025*(shields_stress[I]/ref_shear_stress[I])**14.2 return Q_g @@ -351,8 +354,8 @@ def plot_state(step, time, S_, S_max_, title=True): fig.set_size_inches(3.3*1.1, 3.3*1.1) 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_c/1000., P_c/1000., '--r', label='$P_c$') + # ax_Pa.plot(s/1000., N/1000., '--r', label='$N$') 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$') @@ -364,7 +367,7 @@ def plot_state(step, time, S_, S_max_, title=True): plt.title('Day: {:.3}'.format(time/(60.*60.*24.))) ax_Pa.legend(loc=2) ax_m3s.legend(loc=4) - #ax_Pa.set_ylabel('[kPa]') + # ax_Pa.set_ylabel('[kPa]') ax_Pa.set_ylabel('[MPa]') ax_m3s.set_ylabel('[m$^3$/s]') @@ -375,8 +378,8 @@ 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$') @@ -473,12 +476,12 @@ while time <= t_end: # TODO: Update f_s from fluxes # Determine change in channel size for each channel segment - #dSdt = channel_growth_rate(e_dot, d_dot, 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 + # Update channel cross-sectional area and width according to growth + # rate and size limit for each channel segment S, W, S_max = update_channel_size_with_limit(S, S_old, dSdt, dt, N_c) # Find new water fluxes consistent with mass conservation and local