Theory

This page describes the equations solved by D-Claw. A user interested in the details of the equation derivation is referred to Iverson and George (2014). Those interested in an explanation of the numerical implementation are referred to George and Iverson (2014).

State variables and auxiliary variables

D-Claw considers the flow of material under gravity (\(\vec{g} = (g_x,g_y,g_z)^\mathrm{T}\)) with thickness \(h\), x- and y- directed velocities \(u\) and \(v\), solid-volume fraction \(m\), and basal pore-fluid pressure \(p_b\) over an arbitrary surface \(b\). Optionally, entrainment of erodible sediment and particle size segregation are implemented. The depth of eroded and entrained material is given by \(\Delta b\) (positive indicating entrainment has occurred). Segregation considers two species fractions \(A\) and \(B\) with \(\chi\) representing the fraction of species \(A\).

The variables which vary in space and time are:

Name	Description	Units
\(h\)	Flow depth	meters
\(hu\)	Flow x-momentum: depth times x-directed velocity	meters squared per second
\(hv\)	Flow y-momentum: depth times y-directed velocity	meters squared per second
\(hm\)	Solid volume: depth times solid-volume fraction	meters
\(p_b\)	Basal pore-fluid pressure	kilograms per meter per time squared
\(h\chi\)	Species A volume: depth times species A fraction	meters
\(\Delta b\)	Depth of erosion (material entrained)	meters

The variables which vary in space only are

Name	Description	Units
\(b\)	Topobathymetric surface elevation	meters
\(\Theta\)	Slope angle in x-direction (used only if bed-normal coordinates are enabled)	degrees
\(h_e\)	Initial thickness of erodible material that can be entrained.	meters

Core equations

D-Claw solves the following set of partial differential equations. See Iverson and George (2014) and George and Iverson (2014) for derivation, explanation, and references.

\[\frac{\partial h}{\partial t} + \frac{\partial (hu)}{\partial x} + \frac{\partial (hv)}{\partial y}= \varphi_1\]

\[\frac{\partial (hu)}{\partial t} + \frac{\partial}{\partial x} \left ( hu^2\right ) + \kappa \frac{\partial}{\partial x}\left (\frac{1}{2} g_z h^2\right ) + \frac{\partial (huv)}{\partial y} + \frac{h(1-\kappa)}{\rho}\frac{\partial p_b}{\partial x} = \varphi_2\]

\[\frac{\partial (hv)}{\partial t} + \frac{\partial (huv)}{\partial x} + \frac{\partial}{\partial y} \left ( hv^2 \right ) + \kappa\frac{\partial}{\partial y}\left (\frac{1}{2} g_z h^2\right ) + \frac{h(1-\kappa)}{\rho}\frac{\partial p_b}{\partial y} = \varphi_3\]

\[\frac{\partial (hm)}{\partial t} + \frac{\partial (hum)}{\partial x} + \frac{\partial (hvm)}{\partial y} = \varphi_4\]

\[\frac{\partial p_b}{\partial t} - \rho g_z\varrho u \frac{\partial h}{\partial x} + \rho g_z\varrho \frac{\partial (hu)}{\partial x} + u\frac{\partial p_b}{\partial x} - \rho g_z\varrho v \frac{\partial h}{\partial y} + \rho g_z\varrho \frac{\partial (hv)}{\partial y} + v\frac{\partial p_b}{\partial y}= \varphi_5\]

where \(\varphi_1\), \(\varphi_2\), \(\varphi_3\), \(\varphi_4\), and \(\varphi_5\) represent the source terms.

The bulk density of the flow, \(\rho\), is calculated based on the fluid and solid densities, \(\rho_f\) and \(\rho_s\), respectively.

\[\rho = \rho_s m + (1-m)\rho_f\]

The lateral pressure coefficient \(\kappa\) is typically set to \(\kappa=1\).

The source terms \(\varphi_1--\varphi_5\) are defined as

\[\varphi_1 = \frac{(\rho-\rho_f)}{\rho} D\]

\[\varphi_2 = hg_x + uD\frac{(\rho-\rho_f)}{\rho}-\frac{(\tau_{s,x} + \tau_{f,x})}{\rho}\]

\[\varphi_3 = hg_y + vD\frac{(\rho-\rho_f)}{\rho}-\frac{(\tau_{s,y} + \tau_{f,y})}{\rho}\]

\[\varphi_4 = -Dm\frac{\rho_f}{\rho}\]

\[\varphi_5 = \zeta D - \frac{3}{\alpha h}||(u,v)^\mathrm{T}||\tan\psi\]

where \(\vec{\tau}_s = (\tau_{s,x},\tau_{s,y})^\mathrm{T}\) and \(\vec{\tau}_f = (\tau_{f,x},\tau_{f,y})^\mathrm{T}\) are the shear tractions exerted by the solid phase and fluid phase, respectively, at the base of the flow, \(\alpha\) is the debris’ elastic compressibility, \(\psi\) is the granular dilatancy angle, and \(D\) is the granular dilation rate.

For brevity, the variables \(\varrho\) in and \(\varrho\) are defined.

\[ \varrho = \frac{(\rho_f + 3\rho)}{4\rho }\]

\[ \zeta = \frac{3}{2\alpha h} + \frac{g_z \rho_f(\rho -\rho_f)}{4\rho}\]

\(\alpha\) is calculated based on \(m\), \(h\), and \(p_b\)

\[\alpha = \frac{a}{m(\rho g_z h - p_b + \sigma_0)}\]

where \(a\) and \(\sigma_0\) are constants (typically set to \(a=0.01-0.3\), and \(\sigma_0=1000\) Pa).

An optional definition for \(\sigma_0\) (see setrun documentation) utilizes, instead of a fixed constant, a diagnostic function:

\[\sigma_0 = \frac{a}{2}\rho_f g_z h \frac{(\rho_s-\rho_f)}{\rho}.\]

Selecting this option for the definition of \(\sigma_0\) may provide enhanced stability of numerical simulations but is still undergoing beta testing.

The dilation rate,\(D\), is postulated to be proportional to the excess fluid pressure (see Iverson and George, 2014)

\[D = -\frac{2k}{h\mu}\left ( p_b - \rho_f g_z h \right )\]

where \(k\) is the hydraulic permeability (a function of \(m\)), and \(\mu\) is the effective shear viscosity of the pore-fluid.

A form for \(k\) is

\[k = k_r\exp{\frac{m-m_r}{0.04}}\]

where \(k_r\) and \(m_r\) are the reference permeability and solid volume fraction.

The fluid and solid components of basal resistance \(\vec{\tau}_f\) and \(\vec{\tau}_f\) are calculated as

\[\vec{\tau}_f = \frac{2\mu\left ( 1-m \right )}{h}(u,v)^\mathrm{T},\]

\[\vec{\tau}_s = \sigma_e \tan(\phi + \psi)\frac{(u,v)^\mathrm{T}}{||(u,v)^\mathrm{T}||}.\]

with \(\sigma_e = \rho g_z h - p_b\), the basal effective normal stress.

\(\psi\) is the granular dilatancy angle:

\[c_1 \tan\psi = m - \frac{m_\mathrm{crit}}{1+\sqrt{N}}\]

\(c_1\) is a dimensionless constant (typically 1), \(m_\mathrm{crit}\) is the critical-state solid fraction, and \(N\) is a dimensionless variable defined as

\[N = \frac{\mu \dot{\gamma}}{\rho_s \dot{\gamma}^2 \delta^2 + \sigma_e}.\]

where \(\dot{\gamma}\) is the shear rate and \(\delta\) is a characteristic length scale associated with grain collisions.

Additional equations

The elements of D-Claw described in this section are experimental and may change. They are optional and are not enabled by default.

Segregation

Fully representing the influence of particle-size segregation on flow dynamics requires (1) an expression of how flow behavior results in the segregation of different particle species and (2) an expression of how flow behavior is affected by resulting particle species ratios.

To treat the first element, the influence of flow behavior on the segregation of particle species, D-Claw implements a model developed by Gray and Kokelaar (2010). This model considers the segregation of two particle species (species \(A\) and \(B\)) in a depth-averaged flow with a linear velocity profile defined as

\[u(z) = (1-\beta)\bar{u} + 2\beta\bar{u}\frac{z}{h}\]

where \(\beta\) is a constant between 0 (plug flow) and 1 (simple shear) and \(\bar{u}\) is the depth-averaged velocity. Note that Gray and Kokelaar (2010) use the symbol \(\alpha=1-\beta\) such that \(\alpha=0\) results in simple shear and \(\alpha=1\) results in plug flow. We use the symbol \(\beta\) to disambiguate from the debris’ elastic compressibility values defined above and because we think that a value of 0 for plug flow and 1 for simple shear is more intutitive.

As the mixture shears, species \(A\) moves to the surface of the flow, and is preferentially advected by the flow. The lateral transport of \(\chi\), the fraction of species \(A\) is described by the following equation, which, for brevity, shows only x-directed transport.

\[\frac{\partial}{\partial t}(h\chi) + \frac{\partial}{\partial x}(h \chi u) - \frac{\partial}{\partial x} \left (\beta h \chi u \left( 1-\chi\right) \right)=0\]

The representation of the feedback between the value of \(\chi\) and flow behavior is highly experimental (c.f., Jones et al., 2023). At present the value for the permeability, \(k\), given above is modified as follows:

\[k = k_{chi} k_r\exp{\frac{m-m_r}{0.04}}\]

where

\[\log_{10} kr_chi = 4(\chi-0.5))\]

This implementation results material that is 50% each of Species \(A\) and \(B\) having a permeability of \(k_r\). Material that is 100% species \(A\) will have a higher permeability and material that is 100% \(B\) will have a lower permeability.

Warning

The implementation of the feedback between segregation and flow behavior is likely to change.

Entrainment

Multiple options for entrainment are present in D-Claw. A user specifies a spatially variable value of the thickness of erodible material that is initially available for entrainment, \(h_e\). An entrainment rate, \(\frac{\partial h_e}{\partial t}\) is calculated if \(h_e-\Delta b \gt 0\) (i.e., erodible material remains at a given location). The available erodible material is assumed to have a constant solid volume fraction, \(m_e\).

The options for \(\frac{\partial h_e}{\partial t}\) are:

Name	Value of entrainment_method
Old-style entrainment	0
New-style entrainment (not yet implemented)	1

Old-style entrainment

With entrainment_method==0 the following is done to calculate the depth of entrainment dh associated with a time step dt.

See entrainment.f90 for additional details and context.

!! me, constant value for solid volume fraction of entrainable material
!! rho_s, solid density
!! rho_f, fluid density
!! coeff, manning coefficient
!! h, flow depth
!! m, solid volume fraction
!! vnorm, (u**2+v**2)**0.5 where u, v are x- and y-directed velocity
!! b_remaining, remaining depth of material that may be eroded.

! calculate entrained material density
rhoe = me*rho_s + (1.d0-me)*rho_f

! calculate top and bottom shear stress.
beta = max(1.d0-m, 0.1d0)
visc = beta*vnorm*2.d0*mu*(1.d0-m)/(tanh(h+1.d-2))
gamma= rho*beta*(vnorm**2)*(beta*gz*coeff**2)/(tanh(h+1.d-2)**(1.0d0/3.0d0))
t1bot = visc + gamma + tau
dh = entrainment_rate(t1bot-t2top)/(rhoe*beta*vnorm)
t2top = min(t1bot,(1.d0-beta*entrainment_rate)*tau)

! calculate dh
dh = entrainment_rate*dt*(t1bot-t2top)/(rhoe*beta*vnorm)
dh = min(dh,b_remaining)

References

George, D.L., and Iverson, R.M., 2014, A depth-averaged debris-flow model that includes the effects of evolving dilatancy—II. Numerical predictions and experimental tests: Proceedings of the Royal Society of London. Series A, v. 470, no. 2170, p. 20130820, https://doi.org/10.1098/rspa.2013.0820.

Gray, J.M.N.T., and Kokelaar, B.P., 2010, Large particle segregation, transport and accumulation in granular free-surface flows: Journal of Fluid Mechanics, v. 652, p. 105–137 https://doi.org/10.1017/S002211201000011X.

Iverson, R.M., and George, D.L., 2014, A depth-averaged debris-flow model that includes the effects of evolving dilatancy—I. Physical basis: Proceedings of the Royal Society of London. Series A, v. 470, no. 2170, p. 20130819, https://doi.org/10.1098/rspa.2013.0819.

Jones, R.P., Rengers, F.K., Barnhart, K.R., George, D.L., Staley, D.M., and Kean, J.W., 2023, Simulating Debris Flow and Levee Formation in the 2D Shallow Flow Model D‐Claw: Channelized and Unconfined Flow: Earth and Space Science, v. 10, p. e2022EA002590, https://doi.org/10.1029/2022EA002590.