commit 4bb1e340b01ce37312e12942c05715d30529e2b0
parent ed3e36d656e883a7527ac866dad7e6d1ebffcf57
Author: Anders Damsgaard <andersd@riseup.net>
Date: Wed, 19 Apr 2017 20:52:56 -0400
add seperate version with two-phase sediment transport model
Diffstat:
1 file changed, 431 insertions(+), 0 deletions(-)
diff --git a/1d-channel-wilcock-two-phase.py b/1d-channel-wilcock-two-phase.py
@@ -0,0 +1,431 @@
+#!/usr/bin/env python
+
+# # ABOUT THIS FILE
+# The following script uses basic Python and Numpy functionality to solve the
+# coupled systems of equations describing subglacial channel development in
+# soft beds as presented in `Damsgaard et al. "Sediment plasticity controls
+# channelization of subglacial meltwater in soft beds"`, submitted to Journal
+# of Glaciology.
+#
+# High performance is not the goal for this implementation, which is instead
+# intended as a heavily annotated example on the solution procedure without
+# relying on solver libraries, suitable for low-level languages like C, Fortran
+# or CUDA.
+#
+# License: Gnu Public License v3
+# Author: Anders Damsgaard, adamsgaard@ucsd.edu, https://adamsgaard.dk
+
+import numpy
+import matplotlib.pyplot as plt
+import sys
+
+
+# # Model parameters
+Ns = 25 # Number of nodes [-]
+Ls = 1e3 # Model length [m]
+total_days = 60. # Total simulation time [d]
+t_end = 24.*60.*60.*total_days # Total simulation time [s]
+tol_S = 1e-3 # Tolerance criteria for the norm. max. residual for Q
+tol_Q = 1e-3 # Tolerance criteria for the norm. max. residual for Q
+tol_N_c = 1e-3 # Tolerance criteria for the norm. max. residual for N_c
+max_iter = 1e2*Ns # Maximum number of solver iterations before failure
+print_output_convergence = False # Display convergence in nested loops
+print_output_convergence_main = True # Display convergence in main loop
+safety = 0.01 # Safety factor ]0;1] for adaptive timestepping
+plot_interval = 20 # Time steps between plots
+plot_during_iterations = False # Generate plots for intermediate results
+#plot_during_iterations = True # Generate plots for intermediate results
+speedup_factor = 1. # Speed up channel growth to reach steady state faster
+#relax_dSdt = 0.3 # Relaxation parameter for channel growth rate ]0;1]
+relax = 0.05 # Relaxation parameter for effective pressure ]0;1]
+
+# Physical parameters
+rho_w = 1000. # Water density [kg/m^3]
+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 = 5e-3 # Mean grain size in gravel fraction (> 2 mm) [m]
+D_s = 5e-4 # Mean grain size in sand fraction (<= 2 mm) [m]
+#D_g = 1
+#D_g = 0.1
+
+# Boundary conditions
+P_terminus = 0. # Water pressure at terminus [Pa]
+#m_dot = 3.5e-6
+m_dot = numpy.linspace(0., 3.5e-6, Ns-1) # Water source term [m/s]
+Q_upstream = 1e-5 # Water influx upstream (must be larger than 0) [m^3/s]
+
+# Channel hydraulic properties
+manning = 0.1 # Manning roughness coefficient [m^{-1/3} s]
+friction_factor = 0.1 # Darcy-Weisbach friction factor [-]
+
+# Channel growth-limit parameters
+c_1 = -0.118 # [m/kPa]
+c_2 = 4.60 # [m]
+
+# Minimum channel size [m^2], must be bigger than 0
+S_min = 1e-2
+# S_min = 1e-1
+# S_min = 1.
+
+
+# # Initialize model arrays
+# Node positions, terminus at Ls
+s = numpy.linspace(0., Ls, Ns)
+ds = s[1:] - s[:-1]
+
+# Ice thickness [m]
+H = 6.*(numpy.sqrt(Ls - s + 5e3) - numpy.sqrt(5e3)) + 1.0
+# slope = 0.1 # Surface slope [%]
+# H = 1000. + -slope/100.*s
+
+# Bed topography [m]
+b = numpy.zeros_like(H)
+
+N = H*0.1*rho_i*g # Initial effective stress [Pa]
+
+# Initialize arrays for channel segments between nodes
+S = numpy.ones(len(s) - 1)*S_min # Cross-sect. area of channel segments[m^2]
+S_max = numpy.zeros_like(S) # Max. channel size [m^2]
+dSdt = numpy.zeros_like(S) # Transient in channel cross-sect. area [m^2/s]
+W = S/numpy.tan(numpy.deg2rad(theta)) # Assuming no channel floor wedge
+Q = numpy.zeros_like(S) # Water flux in channel segments [m^3/s]
+Q_s = numpy.zeros_like(S) # Sediment flux in channel segments [m^3/s]
+N_c = numpy.zeros_like(S) # Effective pressure in channel segments [Pa]
+P_c = numpy.zeros_like(S) # Water pressure in channel segments [Pa]
+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) # Total sediment flux [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])
+
+
+def avg_midpoint(arr):
+ # Averaged value of neighboring array elements
+ return (arr[:-1] + arr[1:])/2.
+
+
+def channel_hydraulic_roughness(manning, S, W, theta):
+ # Determine hydraulic roughness assuming that the channel has no channel
+ # floor wedge.
+ l = W + W/numpy.cos(numpy.deg2rad(theta)) # wetted perimeter
+ return manning**2.*(l**2./S)**(2./3.)
+
+
+def channel_shear_stress(Q, S):
+ # Determine mean wall shear stress from Darcy-Weisbach friction loss
+ u_bar = Q/S
+ return 1./8.*friction_factor*rho_w*u_bar**2.
+
+
+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)
+
+ # Non-dimensionalize shear stress
+ shields_stress = tau/((rho_s - rho_w)*g*D_s)
+
+ # import ipdb; ipdb.set_trace()
+ Q_c = 11.2*f_s*W/((rho_s - rho_w)/rho_w*g) \
+ * (tau/rho_w)**1.5 \
+ * numpy.maximum(0.0,
+ (1.0 - 0.846*numpy.sqrt(ref_shear_stress/shields_stress))
+ )**4.5
+
+ return Q_c
+
+
+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)
+
+ # 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 \
+ * numpy.maximum(0.0,
+ (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_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(
+ numpy.maximum(
+ 1./4.*(c_1*numpy.maximum(N_c, 0.)/1000. + c_2), 0.)**2. *
+ 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
+ dSdt = S - S_old # adjust dSdt for clipping due to channel size limits
+ return S, W, S_max, dSdt
+
+
+def flux_solver(m_dot, ds):
+ # Iteratively find new water fluxes
+ it = 0
+ max_res = 1e9 # arbitrary large value
+
+ # Iteratively find solution, do not settle for less iterations than the
+ # number of nodes
+ while max_res > tol_Q:
+
+ Q_old = Q.copy()
+ # dQ/ds = m_dot -> Q_out = m*delta(s) + Q_in
+ # Upwind information propagation (upwind)
+ Q[0] = Q_upstream
+ Q[1:] = m_dot[1:]*ds[1:] + Q[:-1]
+ max_res = numpy.max(numpy.abs((Q - Q_old)/(Q + 1e-16)))
+
+ if print_output_convergence:
+ print('it = {}: max_res = {}'.format(it, max_res))
+
+ # import ipdb; ipdb.set_trace()
+ if it >= max_iter:
+ raise Exception('t = {}, step = {}: '.format(time, step) +
+ 'Iterative solution not found for Q')
+ it += 1
+
+ return Q
+
+
+def pressure_solver(psi, f, Q, S):
+ # Iteratively find new water pressures
+ # dN_c/ds = f*rho_w*g*Q^2/S^{8/3} - psi (Kingslake and Ng 2013)
+
+ it = 0
+ max_res = 1e9 # arbitrary large value
+ while max_res > tol_N_c:
+
+ N_c_old = N_c.copy()
+
+ # Dirichlet BC (fixed pressure) at terminus
+ N_c[-1] = rho_i*g*H_c[-1] - P_terminus
+
+ N_c[:-1] = N_c[1:] \
+ + psi[:-1]*ds[:-1] \
+ - f[:-1]*rho_w*g*Q[:-1]*numpy.abs(Q[:-1]) \
+ /(S[:-1]**(8./3.))*ds[:-1]
+
+ max_res = numpy.max(numpy.abs((N_c - N_c_old)/(N_c + 1e-16)))
+
+ if print_output_convergence:
+ print('it = {}: max_res = {}'.format(it, max_res))
+
+ if it >= max_iter:
+ raise Exception('t = {}, step = {}:'.format(time, step) +
+ 'Iterative solution not found for N_c')
+ it += 1
+
+ return N_c
+ #return N_c_old*(1 - relax_N_c) + N_c*relax_N_c
+
+
+def plot_state(step, time, S_, S_max_, title=True):
+ # Plot parameters along profile
+ fig = plt.gcf()
+ fig.set_size_inches(3.3*1.1, 3.3*1.1*1.5)
+
+ ax_Pa = plt.subplot(3, 1, 1) # axis with Pascals as y-axis unit
+ 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, '.-b', label='$Q$')
+
+ if title:
+ plt.title('Day: {:.3}'.format(time/(60.*60.*24.)))
+ ax_Pa.legend(loc=2)
+ ax_m3s.legend(loc=4)
+ ax_Pa.set_ylabel('[MPa]')
+ ax_m3s.set_ylabel('[m$^3$/s]')
+
+ ax_m3s_sed = plt.subplot(3, 1, 2, sharex=ax_Pa)
+ ax_m3s_sed.plot(s_c/1000., Q_t, '-', label='$Q_{total}$')
+ ax_m3s_sed.plot(s_c/1000., Q_s, ':', label='$Q_{sand}$')
+ ax_m3s_sed.plot(s_c/1000., Q_g, '--', label='$Q_{gravel}$')
+ ax_m3s_sed.set_ylabel('[m$^3$/s]')
+ ax_m3s_sed.legend(loc=2)
+
+ ax_m2 = plt.subplot(3, 1, 3, sharex=ax_Pa)
+ ax_m2.plot(s_c/1000., S_, '-k', label='$S$')
+ ax_m2.plot(s_c/1000., S_max_, '--', color='#666666', label='$S_{max}$')
+
+ ax_m2s = ax_m2.twinx()
+ ax_m2s.plot(s_c/1000., dSdt, ':', label='$dS/dt$')
+
+ ax_m2.legend(loc=2)
+ ax_m2s.legend(loc=3)
+ ax_m2.set_xlabel('$s$ [km]')
+ ax_m2.set_ylabel('[m$^2$]')
+ ax_m2s.set_ylabel('[m$^2$/s]')
+
+ ax_Pa.set_xlim([s.min()/1000., s.max()/1000.])
+
+ plt.setp(ax_Pa.get_xticklabels(), visible=False)
+ plt.tight_layout()
+ if step == -1:
+ plt.savefig('chan-0.init.pdf')
+ else:
+ plt.savefig('chan-' + str(step) + '.pdf')
+ plt.clf()
+ plt.close()
+
+
+def find_new_timestep(ds, Q, Q_t, S):
+ # Determine the timestep using the Courant-Friedrichs-Lewy condition
+ dt = safety*numpy.minimum(60.*60.*24.,
+ numpy.min(numpy.abs(ds/(Q*S),\
+ ds/(Q_t*S)+1e-16)))
+
+ if dt < 1.0:
+ raise Exception('Error: Time step less than 1 second at step '
+ + '{}, time '.format(step)
+ + '{:.3} s/{:.3} d'.format(time, time/(60.*60.*24.)))
+
+ return dt
+
+
+def print_status_to_stdout(step, time, dt):
+ sys.stdout.write('\rstep = {}, '.format(step) +
+ 't = {:.2} s or {:.4} d, dt = {:.2} s '
+ .format(time, time/(60.*60.*24.), dt))
+ sys.stdout.flush()
+
+s_c = avg_midpoint(s) # Channel section midpoint coordinates [m]
+H_c = avg_midpoint(H)
+
+# Water-pressure gradient from geometry [Pa/m]
+psi = -rho_i*g*gradient(H, s) - (rho_w - rho_i)*g*gradient(b, s)
+
+# Prepare figure object for plotting during the simulation
+fig = plt.figure('channel')
+plot_state(-1, 0.0, S, S_max)
+
+
+# # Time loop
+time = 0.
+step = 0
+while time <= t_end:
+
+ # Determine time step length from water flux
+ dt = find_new_timestep(ds, Q, Q_t, S)
+
+ # Display current simulation status
+ print_status_to_stdout(step, time, dt)
+
+ it = 0
+
+ # Initialize the maximum normalized residual for S to an arbitrary large
+ # value
+ max_res = 1e9
+
+ S_old = S.copy()
+ # Iteratively find solution with the Jacobi relaxation method
+ while max_res > tol_S:
+
+ S_prev_it = S.copy()
+
+ # Find new water fluxes consistent with mass conservation and local
+ # meltwater production (m_dot)
+ Q = flux_solver(m_dot, ds)
+
+ # Find average shear stress from water flux for each channel segment
+ tau = channel_shear_stress(Q, S)
+
+ # Determine sediment fluxes for each size fraction
+ 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
+
+ # Determine change in channel size for each channel segment.
+ # Use backward differences and assume dS/dt=0 in first segment.
+ dSdt[1:] = channel_growth_rate_sedflux(Q_t, porosity, s_c)
+ #dSdt *= speedup_factor * relax
+
+ # Update channel cross-sectional area and width according to growth
+ # rate and size limit for each channel segment
+ #S_prev = S.copy()
+ S, W, S_max, dSdt = \
+ update_channel_size_with_limit(S, S_old, dSdt, dt, N_c)
+ #S = S_prev*(1.0 - relax) + S*relax
+
+
+ # Find hydraulic roughness
+ f = channel_hydraulic_roughness(manning, S, W, theta)
+
+ # Find new water pressures consistent with the flow law
+ N_c = pressure_solver(psi, f, Q, S)
+
+ # Find new effective pressure in channel segments
+ P_c = rho_i*g*H_c - N_c
+
+ if plot_during_iterations:
+ plot_state(step + it/1e4, time, S, S_max)
+
+ # Find new maximum normalized residual value
+ max_res = numpy.max(numpy.abs((S - S_prev_it)/(S + 1e-16)))
+ if print_output_convergence_main:
+ print('it = {}: max_res = {}'.format(it, max_res))
+
+ #import ipdb; ipdb.set_trace()
+ if it >= max_iter:
+ raise Exception('t = {}, step = {}: '.format(time, step) +
+ 'Iterative solution not found')
+ it += 1
+
+ # Generate an output figure for every n time steps
+ if step % plot_interval == 0:
+ plot_state(step, time, S, S_max)
+
+ # Update time
+ time += dt
+ step += 1
