Validation of a two-layer depth-averaged model by comparison with an experimental dilute stratified pyroclastic density current

Numerical results of a two-layer depth-averaged model of pyroclastic density currents (PDCs) were compared with an experimental PDC generated at the international eruption simulator facility (the Pyroclastic flow Eruption Large-scale Experiment (PELE)) to establish a minimal dynamical model of PDCs with stratification of particle concentrations. In the present two-layer model, the stratification in PDCs is modeled as a voluminous suspended-load layer with low particle volume fractions (≲10-3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lesssim {10}^{-3})$$\end{document} and a thin basal bed-load layer with higher particle volume fractions (∼10-2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim {10}^{-2}$$\end{document}) on the basis of the source condition in the experiment. Numerical results for the suspended load quantitatively reproduce the time evolutions of the front position and flow thickness in the experimental PDC. The numerical results of the bed-load and deposit thicknesses depend on an assumed value of settling speed at the bottom of the bed load (WsH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${W}_{\mathrm{sH}}$$\end{document}). We show that the thicknesses of bed load and deposit in the simulations agree well with the experimental data, when WsH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${W}_{\mathrm{sH}}$$\end{document} is set to about 1.25×10-2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.25\times {10}^{-2}$$\end{document} m/s. This value of the settling speed is two orders of magnitude smaller than that predicted by a hindered-settling model. The small value of WsH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${W}_{\mathrm{sH}}$$\end{document} is considered to result from decreasing in the effective deposition speed due to the erosion process accompanied by saltating/rolling of particles at the bottom of the bed load.


Introduction
Pyroclastic density currents (PDCs) are a frequent and hazardous process during volcanic eruptions. They occur when a hot mixture of volcanic particles and gas is ejected from the vent, but fails to become buoyant and instead propagates outwards as a ground-hugging gravity current (see the reviews by Branney and Kokelaar (2002), Sulpizio et al. (2014), and Lube et al. (2020)). The flow dynamics of PDCs depends on various factors: eruption conditions such as magma discharge rate (e.g., Bursik and Woods 1996;Dufek and Bergantz 2007;Shimizu et al. 2019;Roche et al. 2021); physical processes of PDCs such as ambient air entrainment, particle sedimentation, and basal friction (e.g., Roche et al. 2008;Andrews and Manga 2012;Lube et al. 2019); and topography (e.g., Esposti Ongaro et al. 2008;Andrews and Manga 2011;Kelfoun 2017). Because of the interplay between these different factors, the fluid dynamical features of PDCs are highly variable and form a wide range of deposit characteristics (e.g., Fisher and Schmincke 1984;Cas and Wright 1987;Branney and Kokelaar 2002).
One major reason behind the wide range of the dynamics of PDCs arises due to strong vertical stratification of the currents with respect to particle concentration (Valentine 1987;Branney and Kokelaar 2002;Burgisser and Bergantz 2002;Breard et al. 2016). PDCs are composed mainly of an upper voluminous dilute turbulent suspension region with low particle volume fractions ( ≲ 10 −3 ; referred to as a "suspended load") and a thinner lower region with higher particle volume fractions ( ∼ 10 −2 -0.5 ). As the particle volume fractions change, the controlling factors of the flow dynamics also change. The flow dynamics of the upper region is controlled mainly by the settling of particles in the current, entrainment of ambient air into the current, thermal expansion of the entrained air, and resistance of ambient air at the flow front (e.g., Sparks et al. 1993;Andrews and Manga 2012;Benage et al. 2016). On the other hand, the flow dynamics of the lower region is controlled mainly by gas-particle interaction, particle-particle collision, frictional interaction between the current and the ground, and deposition/erosion processes at the base (e.g., Roche et al. 2008;Girolami et al. 2010;Lube et al. 2019;Brosch and Lube 2020). The behavior of the whole stratified current is determined by the dynamics of both the upper and lower regions and the interactions between them (i.e., transfers of mass, momentum, and energy from one to the other). The effects of these physical processes on the flow dynamics depend on the source conditions and topography. In particular, the lower region has various characteristics depending on the source conditions (for instance, the particle concentration at the source; e.g., Lube et al. 2015;Breard et al. 2018;Valentine 2020); the lower region behaves as a dense gas-pore pressure-modified (i.e., fluidized) granular flow with very high particle volume fractions ( ∼ 0.4 ; referred to as a "dense underflow"; e.g., Breard et al. 2016;Roche et al. 2016;Lube et al. 2019) or as a flow of saltating/rolling particles with relatively low particle volume fractions ( ∼ 10 −2 ; referred to as a "bed load"; e.g., Valentine 1987;Dufek and Bergantz 2007;Brosch and Lube 2020).
Numerical two-layer depth-averaged models have been developed as a minimal dynamical model to describe global features of stratified PDCs (e.g., Doyle et al. 2008;Kelfoun 2017;Shimizu et al. 2019). In the two-layer models, the continuous stratification of particle concentration and density in PDCs is modeled as upper and lower depth-averaged layers on the basis of the idea that the upper-and lower-region flows in stratified PDCs are controlled by different physical processes. In upper layers, the effects of particle settling, air entrainment, thermal expansion, and frontal air resistance on flow dynamics are mainly taken into account on the basis of experiments of particle-water dilute turbulent suspension flows (e.g., Parker et al. 1987;Bonnecaze et al. 1993;Sparks et al. 1993). In lower layers, the effects of basal friction and deposition on flow dynamics are mainly considered on the basis of experiments of particle-air dense fluidized granular flows (e.g., Girolami et al. 2008Girolami et al. , 2010Roche et al. 2008). The two layers are coupled through mass and momentum exchanges such as inter-layer particle transfer. Although the concept of a two-layer model is useful for systematically assessing the effects of the various physical processes on the flow dynamics and resulting deposits of stratified PDCs, the quantitative agreement of its numerical results with experimental observations needs to be tested.
A community-driven effort is currently underway to compare numerical PDC models with experimental data for the purposes of validation (assessing how well a numerical model represents the physical problem) and benchmarking (comparison of different numerical models with one another) (Valentine 2019;Esposti Ongaro et al. 2020). As a part of the effort, a large-scale experiment was conducted at the international eruption simulator facility (the Pyroclastic flow Eruption Large-scale Experiment (PELE); Lube et al. 2015). This experiment involved the controlled gravitational collapse of a heated suspension of natural volcanic particles and air into an instrumented inclined run-out section. The resulting continuously stratified density current simulated a fully dilute, fully turbulent PDC (i.e., a pyroclastic surge) comprising a thick upper suspended-load region and a thin lower bed-load region. The detailed conditions of the experiment and the characteristics of the spatially and temporally evolving flow structure and deposit are described in Brosch and Lube (2020). The benchmark conditions are described in Supplementary Information S1 in Electronic Supplementary Material (ESM) 1. This paper compares a numerical two-layer PDC model (Shimizu et al. 2019) with the experimental data from PELE for the benchmarking and validation in order to establish a minimal dynamical model of stratified PDCs. We assess how well the two-layer model reproduces the experimental stratified PDC to clarify its applicability and limitation. We also discuss the sedimentation process in the experimental bed load.

Method
We conducted a series of numerical simulations of a twolayer PDC model under the conditions defined in the benchmark. A two-layer PDC flows into run-out sections comprising proximal 6 • inclined and distal horizontal channels at x = 0-9.68 m and x > 9.68 m, respectively, where x represents the distance in a direction parallel to the basal surface ( Fig. 1). We set the source conditions using the experimental data at x = 0 (Profile 1) and compare the numerical results with the experimental data in the distal areas ( x > 0 ) particularly at x = 2.65 and 7.78 m (Profiles 2 and 3, respectively). The basic equations of the two-layer PDC model and the source conditions and input parameters in the simulations are shown below.

Basic equations
The present two-layer model is based on that of Shimizu et al. (2019); the two layers are coupled through mass and momentum exchanges as suspended particles in the upper layer settle into the lower layer, and a deposit progressively aggrades upward from the bottom of the lower (or upper) layer (see Fig. 1b). Shimizu et al. (2019) assumed an axisymmetric PDC spreading radially from the source along the horizontal ground surface, whereas this paper designs a PDC flowing into an inclined one-dimensional channel with a slope angle based on the experimental setting of PELE (see Fig. 1a). The basic equations of the upper and lower layers and the deposit are described below (see Shimizu et al. (2019) for the numerical procedures).
The upper layer is modeled as a suspended load consisting of solid particles, volcanic gas, and entrained ambient air. The basic equations of the suspended load with thickness h(x, t) , velocity u(x, t) , density (x, t) , solid mass fraction n s (x, t) , volcanic gas mass fraction n v (x, t) , air mass fraction n a (x, t) , temperature T(x, t) , and specific heat at constant pressure C p (x, t) are as follows. Conservation of flow mass: Conservation of entrained air mass: Conservation of solid particle mass: (3) t n s h + x n s uh = −n s W s cos

Conservation of flow momentum:
Conservation of flow enthalpy:

Equation of state:
Here t is the time, a (≡ p∕(R a T a ) ) is the density of ambient air, E is the entrainment coefficient (see Eq. (A.1) of Shimizu et al. 2019), W s is the settling speed of solid particles at the bottom of the suspended load, g is the gravitational acceleration, z c is the height of the basal contact, C dc is the basal-drag coefficient of the suspended load, u H is the lower-layer velocity, C pa is the specific heat of air at constant pressure, T a is the temperature of ambient air, C s is the specific heat of solid particles, s is the mean density of solid particles, R a is the gas constant of air, R v is the gas constant of volcanic gas, and p is the pressure. The mass fractions satisfy the condition of n s + n v + n a = 1 . The specific heat of the suspended load at constant pressure is given by C p = n s C s + n a C pa + n v C pv . The above depth-averaged equations provide a good approximation of the dynamics of currents where turbulent mixing is sufficiently intense to maintain vertically uniform concentration (cf. Bonnecaze et al. 1993). To describe realistic suspended-load dynamics, a balance between the buoyancy pressure driving the flow front and the resistance pressure caused by the acceleration of the ambient air at the front (i.e., the front condition): is taken into account, where the subscript N denotes the front and F N is the imposed frontal non-dimensional parameter. Equation (7) is applicable to the gravity currents for a wide range of density differences between the current and the ambient fluid (Ungarish 2009;Shimizu et al. 2017).
The lower layer is modeled as a homogeneous bed load (i.e., a homogeneous flow with a higher particle concentration than that of the suspended load) consisting of solid particles and air. The basic equations of the bed load with thickness h H (x, t) and velocity u H (x, t) are as follows. Conservation of solid particle mass: Conservation of solid particle momentum: Here the subscript H denotes the higher particle concentration flow (i.e., the bed load), sH is the particle volume fraction in the bed load, sD is the particle volume fraction in the deposit, W sH is the settling speed of particles at the bottom of the bed load, z b is the height of the contact between the bed load and the deposit, and C db is the basal-drag coefficient of the bed load. The bed load is assumed to have a constant bulk density ) is the density of the gas phase in the bed load and T 1 is the initial temperature of the upper suspended load. To reproduce the experimentally observed fluid dynamical features of the bed load with particle volume fractions less than sD ( ∼ 10 −2 ), we set the particle volume fraction in the lower layer ( sH ) as a model parameter in the present two-layer model.
The deposit progressively aggrades upward from the bottom of the two-layer PDC. The aggradation rate of material in the deposit is expressed by the deposition speed, D , as Here, the value of D is calculated from the mass balance at the deposition surface as where s * is the particle volume fraction in the bed load ( sH ) or that in the suspended load ( s ≡ n s ∕ s ) and W s * is the effective settling speed of particles from the bed load ( W sH ) or that from the suspended load ( W s ). The aggradation for the suspended load occurs when the particle-settling rate from the suspended load is lower than that from the bed load at the position where the bed load is absent (i.e., the two conditions that the right-hand side of Eq. (8) < 0 and h H = 0 are simultaneously satisfied).

Suspended load
The source of the upper suspended load in the simulations is modeled as a supply of homogeneous mixture at a constant mass flow rate from the inlet boundary x = 0 (Fig. 1). The values of the inlet boundary conditions for the suspended load (i.e., thickness h 1 , velocity u 1 , solid mass fraction n s1 , volcanic gas mass fraction n v1 , temperature T 1 , and mean solid density s ) are obtained from the experimental data at Profile 1 ( x = 0 ) (see Table 1). The values of h 1 , u 1 , n s1 , n v1 , and T 1 (Table 1) are estimated by depth-and time-averaging experimental data of flow velocity, temperature, and particle volume fraction as a function of time t and height z at Profile 1 (see Supplementary Information S1 in ESM 1). As the densities of particles in the experiment depend on their particle sizes (ESM 2), the value of s is estimated by the mass-weighted average where n s,i and s,i are the depth-and time-averaged solid mass fraction and solid density of the i-th particle class obtained from the experimental data at Profile 1 (see ESM 2 for details of the estimation). The inlet mass flow rate (per unit length) of the suspended load is given as 1 u 1 h 1 for time t = 0-4 s and as 0 for t > 4 s (Fig. 1c), where 1 represents the density of the suspended load at x = 0 and is estimated by the equation of state (Eq. (6)). The time interval of 4 s refers to the duration of high-speed camera footage available to obtain the velocity data at the static observer location in the experiment.
The flow dynamics of the suspended load is dependent on three factors: the imposed frontal non-dimensional parameter ( F N ), the basal-drag coefficient ( C dc ), and the settling speed of solid particles into the bed load (or the deposit) ( W s ). The values of F N and C dc (Table 1) are estimated from existing models; F N is based on the theoretical model for steady-state inviscid gravity currents (Benjamin 1968) and C dc is estimated on the basis of the empirical formula of Hager (1988) (cf. Hogg and Pritchard 2004). The value of W s for the polydisperse system (Table 1) is calculated from the mean terminal velocity (see the "Effective values of settling speeds ( W s and W sH )" section for details). These values of F N , C dc , and W s are assessed by comparison of the numerical results of the suspended load with the experimental data at x > 0.

Bed load
In the simulations, a bed load is supplied at the inlet boundary x = 0 at a constant mass flow rate (Fig. 1). The inlet boundary conditions of the bed load (i.e., thickness h H1 , velocity u H1 , and solid volume fraction sH ) are obtained from the experimental data (see Table 1 and Supplementary Information S1 in ESM 1). The value of h H1 (Table 1) is based on the time-series of the bed-load thickness at Profile 1 ( x = 0 ). The value of u H1 (Table 1) is estimated by depth averaging the time-averaged height-dependent fitting function of flow velocity at Profile 1 between heights 0 and h H1 . The value of sH (Table 1) is based on the observation that the bed load has solid volume fractions of ∼ 10 −2 (Brosch and Lube 2020). The inlet mass flow rate (per unit length) of the bed load is given as H u H1 h H1 for time t = 0-4 s and as 0 for t > 4 s (Fig. 1c). The duration of t = 0-4 s roughly represents the time interval within which the bed load existed at Profile 1 in the experiment (see Fig. S1a in ESM 1).
The flow dynamics of the bed load is controlled by two major factors: the basal-drag coefficient ( C db ) and the settling speed of particles at the bottom of the bed load ( W sH ). The value of C db (Table 1) is based on the empirical formula of Hager (1988) (i.e., C db = 0.025Re −0.2 H1 ; cf. Hogg and Pritchard 2004), where Re H1 (≡ H u H1 h H1 ∕ H ) represents the Reynolds number of the bed load at x = 0 , and the bulk dynamic viscosity of the bed load ( H ) is set to 10 −5 Pa s. The effective value of W sH is evaluated below.

Effective values of settling speeds ( W and W )
Among the above parameters, the effective values of settling speeds ( W s and W sH ) play a key role. Previously, the

0.02
Volume fraction of solid particles in bed load sedimentation process from particle-laden density currents has been explained by a simple settling model based on the terminal velocity W T (e.g., Bonnecaze et al. 1993;Choux and Druitt 2002;Girolami et al. 2010). We use the mean terminal velocity, W T = ∑ � n s,i W T,i � , as a parameter for the polydisperse system (see ESM 2 for details of the estimation). Here the terminal velocity of each i-th particle class ( W T,i ) is estimated by Kunii and Levenspiel (1969): where d i is the particle diameter, g (= 1.2 kg/m 3 ) is the density of the gas phase, g (= 10 −5 Pa s) is the dynamic viscosity of the gas phase, and Re s,i (≡ g d i W T,i ∕ g ) is the particle Reynolds number.
Generally, the settling speed of particles in particle-fluid mixture deviates from W T due to the effects of fine (wellcoupled) particles on the mixture density and viscosity as well as that of counter flow associated with particle settling (Marble 1970). These effects can be taken into account by a formula (Richardson and Zaki 1954) as (referred to as the hindered-settling model), where m(= 2 -12 ) is the empirical exponent depending on the particle size. Equation (13) is applicable for a wide range of the particle volume fraction from dilute suspension to fluidized bed (Khan and Richardson 1989;Druitt et al. 2007). It implies that the hindered-settling effects diminish for s * ≪ 1 and that the effective value of settling speed in the suspended load ( W s ) is approximated by W T . Consequently, we use W s = W T for the suspended load (Table 1), which means that the particle size distribution in the suspended load is assumed to stay constant through time and space. We also use W sH = 1 − sH m W T (i.e., Eq. (13)) for the bed load in the reference simulation (referred to as "Run Reference"). We consider that care should be taken in the application of Eq. (13) for the bed load, because the settling speed for the bed load is thought to depend strongly on unknown physical processes in the bed load. Equation (13) does not consider some effects that can play a role at the bottom of particle-laden density currents such as erosion of deposits and saltating/rolling of particles. It also ignores the effect of particle cluster, observed in recent experiments (e.g., Lube et al. 2015;Breard et al. 2016;Weit et al. 2018). For these reasons, in simulations other than Run Reference, we set W sH as a tuning parameter and estimate its value on the basis of fitting the numerical results to the experimental data at x > 0 Table 1 for the value of W sH in the best fitted simulation (referred to as "Run Best-fit")). We also discuss the sedimentation process in the bed load on the basis of the difference between Runs Reference and Best-fit.

Suspended load
The results of the suspended loads in the numerical simulations are almost unaffected by the characteristics of the bed load regardless of W sH , because the bed loads have a negligible effect on the dynamics of the suspended load (e.g., the interfacial drag between the two layers). Here we describe the flow dynamical features of the suspended load in the results reproducing the behavior of the experimental bed load (i.e., Run Best-fit with W sH = 1.25 × 10 −2 m/s).
In the simulation, a suspended load is generated from the inlet boundary ( x = 0 ) and flows into the run-out sections ( x > 0 ) (see Fig. 2 and Supplementary Movie 1 (ESM 3)). The numerical results successfully reproduce the qualitative features of the suspended load observed in the experiment. They develop a typical gravity current structure comprising a leading thick gravity current "head" and a trailing gravity current "body" (cf. Brosch and Lube 2020). The results also agree well quantitatively with the main components of the experimental data. They reproduce the time evolution of the front position in the experiment (Fig. 3a). The flow thickness in the simulations ( h(x, t) ) is consistent with the time evolution of the flow thickness (particularly the body thickness) at Profiles 2 and 3 ( x = 2.65 and 7.78 m) in the experiment (Fig. 3 b and c). A brief description of the numerical suspended load is as follows: it has the flow velocity u ∼ 4-5 m/s, flow density ∼ 1.7-2.3 kg/m 3 , and Richardson number Ri ∼ 0.1 -0.4 (i.e., entrainment coefficient E ∼ 0.008-0.04 ) at x = 0-7.78 m. Although the thickness of the head in the simulation does not perfectly agree with the experimental data ( Fig. 3 b and c), this difference is acceptably small, given the fact that the frontal shape in the simulation is determined by the mechanical balance of the depth-averaged model (Eq. (7)) and that the other fluid dynamical effects are not considered.
We performed a parametric study for a wide range of input parameters (i.e., W s = 0.3-3 m/s, F N = 1-√ 2 , and C dc = 10 −4 -10 −2 ) to assess the values of W s , F N , and C dc in Table 1. The results of the parametric study indicate that the dynamical features of the suspended load are insensitive to C dc . On the other hand, they are highly sensitive to W s and F N . The evolution of the body thickness is primarily affected by the flow density, which in turn depends on W s . The front position strongly depends not only on W s but also on F N . This is because the time evolution of the flow front, as well as the head and body structure, is controlled primarily by the resistance pressure caused by the acceleration of the ambient air at the flow front (i.e., the front condition (Eq. (7)); Shimizu et al. 2017). The above good agreement of the body thickness between the numerical and experimental results (Fig. 3 b and c) implies that the mass-weighted average of terminal velocities (i.e., W T = ∑ � n s,i W T,i � ) used in the present model can explain the effective value of W s in the experiment. The agreement for the front position (Fig. 3a) implies that the theoretical model of F N for steady-state inviscid gravity currents (Benjamin 1968) used in the present model explains the mechanical balance at the flow front in the experiment.

Bed load and deposit
In the simulations, a bed load is generated from the inlet boundary ( x = 0 ) and flows into the run-out sections ( x > 0 ) (see Fig. 2 and Supplementary Movie 1 (ESM 3) for the results of Run Best-fit with the values of the input parameters in Table 1). The bed load obtains the mass and momentum of particles settling from the upper suspended load. The front of the bed load stops when its frontal parallel mass flux becomes zero owing to basal deposition (cf. Shimizu et al. 2019). The deposits progressively aggrade upward from the bottom of the bed load in the proximal area and directly from the bottom of the suspended load in the distal area where the bed load is absent (cf. Regime 2a of Shimizu et al. 2019).
When the setting speed at the bottom of the bed load ( W sH ) is set to the value based on the hindered-settling model (Run Reference with W sH = 5.9 × 10 −1 m/s; i.e., Eq. (13) with m = 12 ), the numerical results do not reproduce the experimental results for the bed load and deposit. The bed load is hardly observed in the simulation: its runout distance is extremely short ( ∼ 0.1 m; see Supplementary Movie 2 (ESM 4)). Furthermore, the numerical result generates an unrealistically thick ( ∼ 8 cm) deposit near x = 0 , which was not observed in the experiment. Fig. 2 Representative numerical results of a two-layer PDC at times t = a 3.0, b 4.5, c 6.0, and d 6.5 s from the beginning of propagation in Run Best-fit. Thicknesses of the suspended load ( h(x, t) ; red), bed load ( h H (x, t) ; blue), and deposit The agreement between the numerical and experimental results for the bed load and deposit is improved when W sH is set to 1.25 × 10 −2 m/s (Table 1; i.e., Run Best-fit). The numerical results reproduce the time-averaged data of the bed-load thickness at Profiles 2 and 3 ( x = 2.65 and 7.78 m) in the experiment (Fig. 4 a and b). The deposit mass in the simulation explains the spatial average of the final deposit mass per unit area in the area of Profiles 1-3 in the experiment (Fig. 4c); the unrealistically thick deposit near x = 0 observed in Run Reference disappears. The fact that the numerical results do not reproduce the exponential decay of the experimental deposit mass with distance (see Fig. 4c) is due to that the deposition speed at the bottom of the bed load becomes constant value for given W sH and sH in the present model (see Eqs. (10) and (11)).
The results of the numerical simulations for a wide range of W sH indicate that, as W sH increases, the slope of the bed-load thickness (i.e., h H ∕ x ) decreases, and the deposit mass (per unit area) derived from the bed load increases. These results allow us to estimate the possible range of W sH from the following experimental observations: (1) the bed-load thickness has almost the same value (i.e., 0.005-0.02 m) at Profiles 1-3 (see Fig. 4 a and b; see Fig. S1a in ESM 1), and (2) the deposit derived from the bed load has a mass of 0.8-5 kg/m 2 at Profiles 1-3 (see Fig. 4c). The estimated value of W sH depends on the uncertainties of other parameters for the bed load such as h H1 , u H1 , sH , and C db ; these uncertainties are caused mainly by depth-and/or time-averaging procedures for the experimental data at Profile 1 (see Supplementary Information S1 in ESM 1 for details of the estimations of these uncertainties).
The results of the sensitivity analyses indicate that the bedload thickness (i.e., Fig. 4 a and b) is dependent on the above four parameters as well as W sH . When the uncertainties of h H1 , u H1 , sH , and C db are taken into account, the results of our model with W sH = 0-5.4 × 10 −2 m/s agree with the experimental data of the bed-load thickness at Profiles 1-3 (i.e., 0.005-0.02 m). On the other hand, our two-layer model indicates that the final deposit mass derived from the bed load is given by sD s sH sD − sH W sH Δt cos (see Eqs. (10) and (11)), where Δt (= 3.30-3.77 s) is the time interval within which the bed load exists at x = 0-7.78 m (see Figs. 4a, b and S1a). In other words, the final deposit mass from the bed load at specific points depends on sH and W sH for given Δt . Accordingly, when the uncertainty of sH is taken into account, the experimental deposit mass (i.e., 0.8-5 kg/m 2 at x = 0-7.78 m) is explained by W sH = 2.1 × 10 −3 -6.0 × 10 −2 m/s. Unifying these two constrains from the bed-load thickness and the deposit mass, we conclude that the range of W sH best explaining the experimental data of both the bed-load thickness and deposit mass is 2.1 × 10 −3 -5.4 × 10 −2 m/s. The circle with error bars in Fig. 5 shows the possible range of W sH (normalized by the mean terminal velocity W T , i.e., W sH ∕W T = 2.8 × 10 −3 -7.2 × 10 −2 ) estimated from the experimental observation for sH = 0.01-0.05 . The estimated range of W sH ∕W T is two orders of magnitude smaller than that predicted by the hindered-settling model (the gray region in Fig. 5); Eq. (13) predicts W sH ∕W T = 5.4 × 10 −1 -9.8 × 10 −1 for sH = 0.01-0.05.

Discussion
In the preceding section, we have estimated the possible range of the effective settling speed in the experimental bed load ( W sH = 2.1 × 10 −3 -5.4 × 10 −2 m/s), which substantially deviates from the value based on the hindered-settling model (Fig. 5). We discuss below the origin of this deviation.
The present two-layer model includes a number of parameters, whereas W sH is the only tuning parameter in the above comparison with the PELE experiment. The deviation in Fig. 5 may result from the propagation of the uncertainties of unknown factors. To eliminate this possibility, we validate our numerical model for the lower layer using a supplementary comparison with an experiment of initially fluidized granular flows, where the hindered settling successfully explains its sedimentation process (Girolami et al. 2010). The results in Supplementary Information S2 (in ESM 1) show that the numerical results of a basal-layer model with the hinderedsettling model (i.e., Eq. (13)) agree well with the experimental data reported by Girolami et al. (2008). Accordingly, the deviation observed in Fig. 5 suggests that a certain mechanism other than hindered settling plays a role in the sedimentation process in the experimental bed load.
There are at least two physical processes related to W sH that are not considered in Eq. (13): particle cluster and erosion. According to experimental studies, particle clusters develop in the mixture with intermediate particle volume fractions ( 10 −2 -10 −1 ), which can substantially increase the effective values of W sH (Weit et al. 2018;Brosch and Lube 2020). In the experiments of particle-laden density currents, on the other hand, the deposition speed at the base of the currents is affected by a complex combination of deposition and erosion processes, and the erosion process can significantly decrease the deposition speed, D (Andrews and Manga 2012;Brosch and Lube 2020). As shown in Eq. (11), W sH is directly  related to the deposition speed D . The decrease in D leads to a small effective value of W sH . In the present case, the results in Fig. 5 (i.e., the fact that W sH estimated from the experiment is two orders of magnitude smaller than the value based on the hindered-settling model) indicate that the erosion effect is greater than the particle-cluster effect.
In the experiment, the deposition speed varied temporally and spatially, being accompanied by saltating/rolling particles and shifting/breaking sandwaves in the bed load (Brosch and Lube 2020). On the other hand, the present model assumes that sH and W sH (hence D ) have constant values regardless of the time and the distance from the source. The experimental results also suggest that the complex particletransport due to the strong shear in the transient region between the suspended and bed loads significantly influences the effective value of W sH (see Brosch and Lube 2020 for details), which cannot be accounted for by the present two-layer model. The disagreement between the numerical and experimental results for the profile of the deposit mass ( Fig. 4c) is considered to result from those spatiotemporal variations in sH and W sH (hence D ) due to the complex combination of deposition and erosion processes.
In this paper, we have focused on the dynamics of stratified PDCs with a bed load generated at the base. As mentioned in the "Introduction" section, the lower region in a stratified PDC flows as either bed load ( sH ∼ 10 −2 ) or dense underflow ( sH ∼ 0.4 ), and this difference in the lower region changes the flow and sedimentation of PDCs (e.g., Lube et al. 2015Lube et al. , 2020Breard et al. 2018). Our two-layer model can predict the behavior of stratified PDCs only when the inlet source conditions of the lower region (i.e., h H1 , u H1 , and sH ) as well as the effective vertical mass fluxes of particles (i.e., W s and W sH ) and the effective basal friction (i.e., C db ) are provided. Future works will attempt to develop additional models to determine these parameters for cases where dense underflow or bed load develops at the source (e.g., Breard et al. 2018;Lube et al. 2019;Valentine 2020).

Summary
Numerical results of a two-layer depth-averaged model of PDCs with stratification of particle concentrations were compared with an experimental dilute stratified PDC generated at PELE. In the numerical simulations of the present two-layer model, the stratification in PDCs is modeled as a voluminous suspendedload layer with low particle volume fractions ( ≲ 10 −3 ) and a thin lower bed-load layer with higher particle volume fractions ( ∼ 10 −2 ) on the basis of the experimental source conditions. By fitting the numerical results to the experimental data for the bed load and deposit, the settling speed at the bottom of the bed load ( W sH ) has been estimated to be two orders of magnitude smaller than that predicted by the hindered-settling model. The small effective value of W sH is considered to result from the erosion process accompanied by saltating/rolling of particles at the deposition surface. Further understanding of PDC dynamics based on the two-layer model would require similar comparisons under various source conditions (e.g., those where a dense underflow develops).

Appendix
We provide the list of the mathematical symbols used in this paper (Table 2).