A Non-Equilibrium Sediment Transport Model for Dam Break Flow over Moveable Bed Based on Non-Uniform Rectangular Mesh

Gangfeng Wu 1, Zhehao Yang 1, Kefeng Zhang 1,*, Ping Dong 1,2 and Ying-Tien Lin 3 1 School of Civil Engineering and Architecture, Ningbo Institute of Technology, Zhejiang University, Ningbo 315100, China; zjdxwgf@gmail.com (G.W.); haomingzezhe@126.com (Z.Y.); Ping.Dong@liverpool.ac.uk (P.D.) 2 School of Engineering, University of Liverpool, Liverpool L69 3GH, UK 3 Ocean College, Zhejiang University, Hangzhou 310058, China; kevinlin@zju.edu.cn * Correspondence: kfzhang@hotmail.com; Tel.: +86-138-5785-6003


Introduction
Dam break flow is one of the research focuses in hydraulic engineering due to its tremendous destructive consequences. For predicting such flow the depth averaged shallow water models have been widely used in the real-life applications [1]. Early models assume the bed elevation is fixed and are only applicable to dam break flows over fixed beds [2][3][4][5][6]. However, significant sediment movements and morphological changes are commonly observed in some catastrophic flood events caused by dam failures. Therefore, it is necessary to account for the sediment transport and bed evolution when simulating dam break flow. In recent years, some numerical models have been proposed in the literature to simulate dam break flows on movable beds [7][8][9][10][11][12][13][14][15][16][17][18][19][20]. In order to improve the accuracy of simulating dam-break flow-induced sediment transport, the numerical models are evolving from equilibrium description to non-equilibrium description of sediment transport [21,22]. Equilibrium sediment transport models based on the assumption that the bed load or the total load are instantaneously adapted to equilibrium state have been widely used for morphodynamic problems [17][18][19][20]. However, the assumption of local equilibrium is not accurate for fast geomorphic transients and may lead to unrealistic results [23,24]. Compared to equilibrium models, non-equilibrium sediment transport models which take into account spatial and temporal lag effects between transport capacity and local hydrodynamic conditions are more appropriate for strong sediment transport caused by dam break flow. The main challenges in numerical aspects of developing a robust model are preserving the positivity of water depth and satisfying the well-balanced property or the C-property [25,26].
Based on the projection step, the Godunov-type schemes generally fall into two categories: upwind schemes (Riemann solver schemes) [9,12,14,[17][18][19][20]27] and central schemes [6,[28][29][30]. Central schemes do not require any Riemann problem solvers and characteristic decomposition, which is a powerful black-box solver for the hyperbolic partial differential equations (PDEs). Kurganov et al. [31] proposed a numerical scheme, which belongs to Godunov-type Riemann-problem-solver-free central schemes, to solve 1D and 2D shallow water equation system with discontinuous bottom topography. In order to reduce the numerical dissipation, the one-sided local speeds are adopted to estimate of the width of the Riemann fans. So the proposed central scheme has an upwind nature and is called central-upwind scheme. They proved that the central-upwind scheme is well-balanced, is able to preserve the positivity of water depth when the Courant number is less than 0.25, and has key features of the Godunov-type central schemes, i.e., robustness, efficiency and simplicity. Liu et al. [15] further shows that the central-upwind scheme is less expensive than other popular Riemann solvers while still capable of tracking wet/dry fronts accurately. The central-upwind schemes are widely used in numerical models to simulate dam break flow over fixed [6,30] and movable bed [15,21,29,32], and satisfactory simulation results are achieved. However, most of these existing models are based on uniform Cartesian grids or triangular mesh, and limited studies are dedicated to non-uniform rectangular or quadtree grids.
Uniform Cartesian grids are widely used in computational simulation of free surface flows due to it's simplicity. However, since the same resolution has to be adopted in the whole domain, large number of computational nodes and high computational cost are needed when high-resolution grids are used. In fact, in large-scale simulations, the high-gradient regions where fine grids are required are much less than the whole computational domain [33]. In order to improve computational efficiency and at the same time maintain the required accuracy, non-uniform grids have been recently used to simulate dam break flow over movable bed, either in the form of non-uniform triangular mesh [12,15,34] or rectangular mesh [35,36]. Compared to the triangular mesh, non-uniform rectangular mesh is more convenient to achieve high order accuracy in space, while keeps most advantages of the rectangular mesh [37].
In this paper, a coupled hydrodynamic and non-equilibrium sediment transport model on non-uniform rectangular mesh is developed to simulate dam break flow over movable bed. Following Wu [1], the enhanced shallow water and sediment transport equations are adopted, including the terms that account for the mass and momentum exchange between flow phase and sediment phase. The interface flux is computed by the central-upwind scheme with second order accuracy in space. The quantitative performance of the model has been investigated by comparing the results with measurement data in four dam break laboratory experiments.

Hydrodynamic Equations
The two-dimensional shallow water equations, which are derived by integrating the Navier-Stokes equations over depth with hydrostatic pressure assumption, are widely used to simulate dam break flow over fixed bed. However, strong sediment transport and rapid bed changes, which could generate significant effect on water flow, are occurring during dam break flow over movable beds. Thus, the enhanced shallow water equations including the mass and momentum exchange between flow phase and sediment phase to describe the hydrodynamic procedure of the dam break flow over movable bed are required, which are [1]: ∂ρhv ∂t where t = time; x and y = the space coordinate; h = flow depth; u and v = flow velocity in x and y directions, respectively; z b = bed elevation; g = gravity acceleration; ρ = density of water and sediment mixture; C t = the volumetric concentration of total-load sediment; ρ s and ρ w = density of sediment and water, respectively. ρ b = the density of the water and sediment mixture in the bed surface layer, ; p = the porosity of the bed material; n = Manning's coefficient for channel bed.

Morphodynamic Equations
Simulation of the sediment transport processes induced by dam break flow can be achieved by calculating the total-load directly [1]. Therefore, the sediment transport and bed evolution equations are given as: where q t = total-load sediment transport capacity, q t = q b + q s ; q b and q s = the equilibrium bed load and suspended sediment transport rate, respectively; L = the non-equilibrium adaption length.

Non-Equilibrium Adaptation Length
The non-equilibrium adaptation length L is a characteristic distance for sediment to adjust from a non-equilibrium state to the equilibrium state, which has been investigated by many researchers [38][39][40]. The non-equilibrium adaptation length L could be determined by [1,10]: where L b = the adaptation length of bed-load, which is related to the dominant bed form or sediment transport scale; α = the adaptation coefficient of suspended load, whose values mostly range between 0.25 to 5.0. It should be mentioned that these two parameters are usually determined through calibration against experimental data. w s = w s0 (1 − C t ) m = settling velocity of sediment particles; m = exponent against experimental data indicating the hindered influence of high sediment concentrations on settling velocity. w s0 = the settling velocity of a single sediment particle in quiescent water, which is given by [41]: where d = the diameter of the sediment particles; ν = the kinematic viscosity of water.

Sediment Transport Rate
The commonly-used formulas of Wu et al. [42], are adopted to calculate the equilibrium transport rate of bed load q b and suspended load q s : where τ be = effective bed shear stress, τ be = (n /n) 3/2 τ b ; n = Manning's coefficient related to grain roughness, n = d 1/6/ 20; n = Manning's coefficient for channel bed; τ b = bed shear stress, ; τ c = the critical bed shear stress for incipient sediment motion, which is calculated using the relationship of Wu et al. [42]:

Non-Uniform Rectangular Mesh
The non-uniform rectangular mesh, as shown in Figure 1, is used in this study. To obtain this non-uniform rectangular mesh, a straightforward method proposed by Liang [33] is implemented: Firstly, the whole computational domain is discretized by a coarse uniform rectangular grid. Then, a specific subdivision level is allocated to the coarse background cell and the coarse cell will be divided into 2 lev × 2 lev fine elements, where lev = a subdivision level. Finally, regularize the grid to make sure that any cell is at most two times bigger or smaller than its neighbors (2:1 rule), as shown in Figure 1. After regularization of the grid, neighboring cells for an arbitrary cell could be determined by a simple algebraic relationship [33]. Moreover, in order to simplify the computer programing and reduce memory usage, the computational cells are numbered in a 1D sequence and a non-staggered approach is used, which defines all primary variables h, u, v, and C t at cell centers and represent the average values of each cell.

Numerical Method
In this study, a central upwind scheme based on finite volume method is used to solve the equations for the dam break flow problem, which involves mix flow regimes and sediment transport.
Second-order accuracy in space is achieved by the nonnegative water depth reconstruction method proposed by Liang [4]. More details of the numerical method is given in the following paragraphs. In order to present this numerical method clearly and compactly, Equations (1)-(4) is reformulated into a shallow water equation system and written in vector form:

∂Q ∂t
where η = water surface elevation, η = h + z b . Q = a vector of the conserved variables. F and G = the flux vectors in x and y direction, respectively. S = the source term vectors which include bed slope terms and bed friction terms.

Finite Volume Discretization
Integrating of Equation (11) over non-uniform rectangular control element i shown in Figure 2, and applying the Green theorem and the Euler time-integration method, one can get the following discretized equation: where the superscript o = the time level and subscript i = the index of the cell. ∆t = the time step. ∆x and ∆y = cell sizes in x and y directions, respectively. F E , F W , G N and G S = the fluxes through the east, west, north and south cell interfaces, respectively. S i = the source term at the cell center. It should be mentioned that when evaluating the flux through the interface of two cells with different sizes, special attention must be taken to maintain the mass and momentum conservation. As shown in Figure 2, two smaller rectangular cells locate at the right side of the cell i. The flux vector F E through the eastern interface could not be directly calculated, but given by F E = (F E1 + F E2 )/2, where F E1 and F E2 are the fluxes through the western interface of two smaller neighboring cells. This step could preserve the conservative property of the finite volume scheme.

Nonnegative Water Depth Reconstruction
To ensure the C-property and deal with wetting and drying accurately in a uniform rectangular grid Liang [4] proposed a nonnegative water depth reconstruction method. In this section the same procedure has been extended in 2D non-uniform rectangular grids to achieve second order accuracy in space and C-property.
Firstly, the Riemann states, the conserved variables (η, hu, hv) and water depth h at interface, are reconstructed from the available flow information at the cell center. For a cell i the left Riemann states at the cell interface i + 1/2 in x direction (eastern face) are given by: where ϕ i = a slope limiter evaluated separately for considered variables. In this study the minmod limiter is used. The velocity components of interface are calculated as: When the cell is dry (h < 10 −8 m), the velocities are set to be zero and not calculated by Equation (17). The left bed elevation at the cell interface i + 1/2 in x direction is given by The right Riemann states η R i+1/2, x at the cell interface i + 1/2 in x direction can be obtained using the similar method. Then, a single value of bed elevation at the cell interface i + 1/2 is defined following Audusse et al. [43]: In order to preserve the reconstructed water depth to be nonnegative, the face values of water depth are corrected as: Finally, the rest Riemann states can be recalculated by: In general the above reconstruction procedure will not affect the C-property of the numerical scheme. However, further modifications need to be implemented if the bed elevation of the dry cell is higher than the water surface elevation of the neighbor wet cell [4]: This treatment makes the reconstructed bed elevation equal to water level at the interface of the wet cell and the dry cell, which ensures the numerical scheme to be well balanced for the dry-bed application.

Numerical Fluxes Calculation
The central-upwind scheme proposed by Kurganov et al. [31] is applied to calculate mass and momentum flux (the first three components of F and G) at interface. Taking the calculation of flux through the cell interface i + 1/2 in x direction as an example, the mass and momentum flux can be calculated as: where Q L i+1/2,x and Q R i+1/2,x = the reconstructed Riemann states on the left and right side of the cell interface i + 1/2 in x direction, respectively. a + i+1/2,x and a + i+1/2,x = one-sided local wave speeds on the left and right side of the cell interface i + 1/2 in x direction, respectively and are calculated by: The sediment transport flux (the fourth component of F and G) at interface is determined by upwind scheme. Thus, the sediment flux can be calculated as: where (C t ) L i+1/2,x and (C t ) R i+1/2,x = the volumetric sediment concentration in cells located on the left and right hand sides of interface i + 1/2 in x direction, respectively. F i+1/2,x (1) = mass flux calculated by Equation (29).

Treatment of Bed Slope and Friction Source Terms
The bed slope source terms are simply evaluated by central difference method and no special treatments are needed to ensure the well balance property. Taking the bed-slope source term in x direction as an example: Similarly, the bed slope term for y direction can be discretized using the same method.
Since the water depth is present in the denominator of the friction source term, the fully explicit scheme for friction term may lead to instability of the numerical scheme when the water depth becomes very small. Therefore, the semi-implicit scheme, which is more stable than explicit scheme, is used to discretize the friction term:

Bed Deformation Calculation
The bed deformation is updated at each time step using the following explicit equation:

Model Stability
The above numerical model is explicit and its stability is controlled by the Courant-Friedrichs-Lewy (CFL) condition. Therefore, the time step can be expressed as: where Cr = the Courant number specified in the range 0 < Cr < 0.25 and Cr = 0.2 is adopted in all simulations in this study.

Model Validation
In this section, four laboratory experiments are used to test the performance of the proposed model. The first two cases are adopted to test the capability of the developed model to simulate dam break flow over fixed irregular bed with wetting and drying. Then, the model is validated against two benchmark tests of dam break flow over a dry and erodible bed to demonstrate it can simulate sediment transport and bed erosion with reasonable accuracy. For all tests, g = 9.8 m/s 2 and ρ w = 1000 kg/m 3 .

Partial Dam Break Experiment
The partial dam break experiment conducted by Fraccarollo and Toro [2] is used to evaluate the accuracy of the numerical model to simulate dam break flow over a fixed flat bed with wetting and drying. The detailed setup of the experiment is shown in Figure 3, in which the whole domain is 4 m long and 2 m wide. The dam with a breach of 0.4 m is located 1 m away from the upstream end. The initial water surface elevation in the reservoir surrounded by solid walls is 0.6 m, while the downstream flood plain is initially dry and has open boundaries in three sides. The breach in the middle of the dam is removed instantaneously at t = 0 s and the water in the reservoir rushes onto the downstream flood plain. Five gauges are placed at different locations to record the time history of water depth. The coordinates of five gauges are given in Table 1.
In this simulation, the computational domain is firstly discretized by 40 × 20 rectangular cells with a uniform cell size of 0.1 m. The area near the breach with high water surface gradient is subdivided and the subdivision level is set to be 3. Figure 4 presents the non-uniform rectangular grids for the whole domain. The total simulation time is 10 s and the time step is adaptively determined by the CFL condition. Figure 5 presents the comparisons between the calculated and measured time history of water depth at five gauge stations. It can be found that the numerical water depth agree quite well with the measured data, which demonstrate that the model is able to simulate dam break flow over fixed flat bed with wetting and drying accurately.

Simulation of Dam Break Flow over Triangular Obstacle
The laboratory experiment of dam break flow over triangular obstacle, which was performed in the Concerted Action on Dam-Break Modeling (CADAM) project [44] and had been widely used as an benchmark test by many researchers [14,21,36,45] , is considered to validate the current model's ability to simulate dam break flow with wet/dry front over fixed irregular bed. Figure 6 presents the configuration of the experiment, in which the flume is 38 m long and the dam is located 15.5 m away from the upstream end. A symmetric triangular obstacle of 0.4 m high and 6 m long is installed at 13 m downstream of the dam. The initial water level in the reservoir is 0.75 m and the downstream area is a dry bed. The dam suddenly fails at t = 0 s and the water in the reservoir flows to the downstream area. As shown in Figure 6, six gauges are located at 4 m(G4), 8 m(G8), 10 m(G10), 11 m(G11), 13 m(G13) and 20 m(G20) away from dam to record the time history of water depth. The boundary conditions for the upstream and downstream ends are solid wall and free outlet, respectively. Figure 7 presents the computational grid for decomposing the whole domain, in which the cell size of the background grid is 0.5 m and subdivision level 3 is applied to the triangular obstacle. As suggested by Liang and Marche [45], the Manning's coefficient is set to be 0.0125. The simulation is run for 90 s. Figure 8 shows the calculated time history of water depth at six gauges. It can be seen the water depths and arrival time of the wave are well predicted except for gauge G20. Compared to the measurement data, the water depth is overestimated by the present model and the similar discrepancy was also presented in other models [45][46][47], which used different numerical schemes. One possible explanation for this deviation is the depth integrated governing equations is not suitable for the complex water flow phenomenon between triangular obstacle and outlet. In order to further evaluate the differences between the calculated and measured water depths quantitatively, the average discrepancy (L 1 norm) is defined as Overall, the calculated result agree with measurements quite well and the wet/dry fronts are well captured, which confirms the effectiveness and robustness of our model to deal with dam break flow over fixed irregular bed with wet/dry fronts.

Dam-Break Flow in an Erodible Channel with a Sudden Enlargement
In order to validate the ability of the present model to deal with dam-break flow over movable bed, an experiment of dam-break flow in an erodible channel with a sudden enlargement conducted in the laboratory of Université catholique de Louvain (UCL) in Belgium [48,49] , which had been used by many researchers [14,[17][18][19]22,27] to test their models, is simulated here. As shown in Figure 9, the experimental flume is 6 m long with a sudden enlargement of width from 0.25 m to 0.5 m at 4 m downstream of the inlet. The whole flume is covered by a layer of 0.1 m thick, fully saturated sand with a median diameter of 1.82 mm. The density of the sand is 2680 kg/m 3 , and the deposit porosity is 0.47. The dam is installed at 1 m upstream of the sudden expansion and the initial water depth in the upstream and downstream of the dam are 0.25 m and 0 m, respectively. Eight ultrasonic gauges P1-P8 are placed at downstream of channel to record the water level in a period of 12 s, and the coordinates of gauges are given in Table 2. At the end of the experiment, the bed elevations at nine cross sections cs1-cs9 located every 5 cm from 10 to 50 cm after the sudden expansion are measured using a laser sheet imaging technique.
For the simulation, the computational domain is firstly discretized by coarse rectangular grids with dimension of 0.05 m × 0.01 m. Then, the area located between 2.5 m and 4.5 m away from upstream end is refined and the subdivision level is set to be 2. Figure 10 presents the computational grids for the whole domain. As suggested by Wu et al. [27], the calibrated suspended-load adaptation coefficient α is equal to 4.0, the adaptation length of bed-load L b is set to be 0.025 m, and the Manning's coefficient is given as 0.025. Figure 9. The experimental configuration of dam break flow over movable bed in a sudden-expanded channel (dash lines: measurement cross sections; circles: ultrasonic gauges) [48,49].  Figure 11 shows the calculated and measured water level at gauges P1, P2, P3, P6 and P7 in a period of 12 s. It can be found that the calculated water level agree with the measurement data quite well. The average discrepancies L 1 of water level at five gauges (i.e., P1, P2, P3, P6 and P7) are 0.0130 m, 0.0064 m, 0.0068 m, 0.0066 m and 0.0064 m, respectively. Figure 12 presents the comparison of the calculated and measured bed elevation at five cross sections: cs1, cs3, cs5, cs7 and cs9. At cross sections cs1 and cs3, the deposition height along the left side of the channel is underestimated, but the general erosion and deposition patterns are reproduced. At cross sections cs5, cs7 and cs9, the calculated bed elevation profiles are in good agreements with measured values. The average discrepancies L 1 of bed elevation at five cross sections (i.e., cs1, cs3, cs5, cs7 and cs9) are 0.0075 m, 0.0067 m, 0.0052 m, 0.0039 m and 0.0028 m, respectively. Figure 12 also shows the calculated bed profile by Juez et al. [18] and Wu et al. [27]. It can be found that the calculated result by Wu et al. [27] and present model, both of which are non-equilibrium sediment transport model, is better than equilibrium model proposed by Juez et al. [18]. So the spatial and temporal lag effects between transport capacity and local hydrodynamic conditions should be taken into account in strong sediment transport simulation induced by dam break flow. In general, the experiment of dam-break flow in an erodible channel with a sudden enlargement can be simulated with an reasonable accuracy, which demonstrate the ability of the present model to deal with dam-break flow over movable bed.

Partial Dam-Break Flow over Mobile Bed
In this section, the benchmark experiment of partial dam-break flow over mobile bed, which was also performed at UCL-Belgium [16] and has been used to test morphological models by many researchers [17,21,27,35,50], is adopted here to further verify the ability of the model to simulate bed evolution of a mobile bed caused by partial dam-break flow. As shown in Figure 13, the flume is 3.6 m wide and about 36 m long. The 1 m wide gate located between two impervious blocks is lifted immediately to represent partial dam-break. The gate center is set to be the origin of the Cartesian coordinate system. The fixed bed of the flume is paved with an 85 mm thick sand layer extending from 1 m upstream of the gate to 9 m downstream of the gate. The medium diameter, specific gravity and the initial porosity of the fully saturated sand lay are 1.61 mm, 2.63 and 0.42, respectively. The whole duration of the experiment is 20 s, which ensures the dam-break wave cannot reach the outlet and avoids the reflection effect from the downstream end of the flume. The initial water level in the upstream reservoir is 0.47 m and the downstream area is dry bed. Four pairs of ultrasonic gauges G1-G8 are symmetrically placed at downstream of the gate to measure the water level, and the coordinates of gauges are given in Table 3. Bed profiles along three longitudinal sections (y = 0.2 m, 0.7 m, and 1.45 m) are measured at the end of the experiment. Figure 13. The experimental configuration of partial dam break flow over movable bed [16]. As shown in Figure 14, a coarse mesh consisting of 360 × 36 rectangular grids (∆x = ∆y = 0.1 m) in longitudinal and lateral directions is firstly used to cover the whole domain. Then, the domain between 1.0 m upstream of the dam and 4.0 m downstream of the gate is refined and the subdivision level is 2. The Manning's coefficient for the fixed bed and movable bed are equal to 0.01 and 0.0165, respectively [16]. Following the suggestion of Wu et al. [27], the calibrated value of the suspended-load adaptation coefficient α and the adaptation length of bed-load L b are set to be 4.0 and 0.025, respectively. Figure 15 presents the calculated and measured bed elevation at the end of the experiment. Both of the measured and calculated results show that significant erosion occurs around the downstream of the gate which lead to a large scour hole, and some of the eroded sediment deposited along the sidewalls downstream of the hole. Although the erosion depth near the downstream of the dam is overestimated by the present model due to strong 3D features of the flow in this area, the computed and measured main erosion and deposition patterns are rather similar. Further quantitative comparisons between calculated and measured bed profile along three longitudinal sections (y = 0.2 m, 0.7 m and 1.45 m) are presented in Figure 16. Two sets of experimental bed elevation data (Experiment1 and Experiment2) measured using the same experimental setup are included in the figure. Strong variations between two experiments are observed, especially in the area between x = 0 m and x = 2 m. The possible reason for the difference between two repeated experiments is the random process of the turbulent flow. The average discrepancies L 1 of bed elevation along three longitudinal sections compared to the Experiment1 are 0.010 m, 0.012 m and 0.016 m, respectively. Compared to the Experiment2, the average discrepancies L 1 of bed elevation are 0.015 m, 0.011 m and 0.006 m, respectively. It should be noted that the deviations between calculated and measured bed elevation were in a similar magnitude with the discrepancy of the measured bed elevation among different experiment runs. In this respect, the computed bed elevations were generally in reasonable agreement with the measurements.  Figure 17 shows the calculated and measured water level at four pairs of gauges in a period of 20 s. It can be found that the calculated water levels at G1 and G4 are about 40 mm lower than the measurement, which are also reported by other researchers [27,35]. One possible reason is that strong

Conclusions
In this paper, a depth-averaged 2D finite volume model based on non-uniform rectangular mesh has been developed to simulate dam break flow over movable beds. The model uses the enhanced shallow water equations which consider the effects of the mass and momentum exchange between flow phase and sediment phase. The non-equilibrium transport of total load is adopted in the sediment transport model. In the model, the non-uniform rectangular mesh is used to refine the grid in computational domains with complex flow structures, such as dam break locations. The fluxes at interface are calculated by the central-upwind scheme, which belongs to the class of Godunov-type Riemann-problem-solver-free central schemes and is able to capture the moving wet/dry fronts accurately. The nonnegative water depth reconstruction method is used to achieve second-order accuracy in space and C-property.
The developed model is firstly tested using two experimental cases of dam break flow over fixed irregular bed. The water depth and advance front celerity are well predicted by the present model. Then the model is further validated against two laboratory experiments of dam-break flow over movable bed. The comparison with the measurement data shows that the proposed model is able to calculate the water surface and final bed level generally well. In this study, the model is applied in the flume-scale dam-break flood with flow-sediment interaction. Since in reach-scale applications, the flow conditions and bed configurations are more complex, a further study of present model in real-life applications will be addressed.
Author Contributions: G.W. carried out the numerical simulations and wrote the first draft of the manuscript. Z.Y. and Y.L. contributed to the data analysis. K.Z. and P.D. conceived and supervised the study and edited the manuscript. All authors reviewed the manuscript.