Combination of Discrete Element Method and Artiﬁcial Neural Network for Predicting Porosity of Gravel-Bed River

: In gravel-bed rivers, monitoring porosity is vital for ﬂuvial geomorphology assessment as well as in river ecosystem management. Conventional porosity prediction methods are restricting in terms of the number of considered factors and are also time-consuming. We present a framework, the combination of the Discrete Element Method (DEM) and Artiﬁcial Neural Network (ANN), to study the relationship between porosity and the grain size distribution. DEM was applied to simulate the 3D structure of the packing gravel-bed and ﬁne sediment inﬁltration processes under various forces. The results of the DEM simulations were veriﬁed with the experimental data of porosity and ﬁne sediment distribution. Further, an algorithm was developed for calculating high-resolution results of porosity and grain size distribution in vertical and horizontal directions from the DEM results, which were applied to develop a Feed Forward Neural Network (FNN) to predict bed porosity based on grain size distribution. The reliable results of DEM simulation and FNN prediction conﬁrm that our framework is successful in predicting porosity change of gravel-bed.


Introduction
The significantly different size of fine sediment and coarse bed material is the reason for the change in the void space of gravel-bed rivers. If finer sediments occupy the interstitial spaces of coarser bed materials, the void space of the bed material decreases. If finer sediments are supplied at a small rate to a riverbed covered by a completely developed armor coat, the fine sediments can infiltrate into the interstitial spaces of the coarser bed material and move into void spaces. As a result of fine sediment exchange in void space, the ratio of void space over the total amount of gravel-bed and the porosity change [1,2]. The study of porosity variation plays an important role for fluvial geomorphology assessment and in river ecosystem management. From a river management point of view, the amount of fine sediment stored in the void space of the gravel-bed up to 22% may be neglected and porosity can be considered constant [1]. The impact the void spaces of gravel-bed have on habitats for fish and aquatic species are substantial, and are of important considerations in assessing changes in the void structure of bed materials [3]. Morphologically, the porosity strongly influences the rate of bed level changes [4,5]. The reduction of porosity due to the excess of fine sediment prevents hydrologic exchange surface and underground water [6]. Hence, the variation of the bed porosity should be considered in hydromorphological models for gravel-bed rivers and also in defining hydrological conductivity of modeling of ground water [7,8].
Porosity predictions for fluvial gravel-bed mixtures are classified into two types: Empirical prediction and theoretical prediction. The empirical prediction focuses mainly on the relationship change in the bed elevation is calculated based on the mass conservation in the active layer. The grain sorting in the active layer occurs under the flow interaction and exchange processes occurring between the active layer and the active stratum layer. The grain sorting of the active stratum layer is only caused by the exchange process between these two layers. Further, the exchange rate for the fine fraction of size class can be quantified, based on the empirical equation of Toro-Escobar et al. [26]. After defining the grain size distribution, the bed porosity is calculated using the semi-empirical equations proposed by Reference [16]. This equation required an empirical coefficient and was verified for binary mixture. Therefore, the porosity calculation model needs to reduce computational time as well as increase accuracy by taking into account the change of multi size fraction of bed material.
In this study, a framework combining the Discrete Element Method (DEM) and Artificial Neural Network (ANN) was introduced to predict the porosity of gravel-bed rivers. Firstly, DEM was applied to simulate the 3D bed structure formed by fine sediment infiltrating into gravel packing bed. Then, the results of porosity and sediment distributions obtained by the DEM were compared with experimental results to confirm the capacity of the DEM model. An algorithm was developed for calculating porosity and grain size distribution by the depth of flumes and along the flume. Finally, datasets obtained by the DEM model were used to design an ANN model, called Feed Forward Neural Network (FNN), for predicting the bed porosity of gravel-bed river.

Methodology
In this study, we weave multiple techniques together into one platform, which is shown in Figure 1. The detailed implementation of each component is presented in the following subsections.

Discrete Element Method (DEM)
The Discrete Element Method (DEM) was initially suggested by Cundall and Strack [27] to model the mechanical behavior of granular flows and to simulate the forces acting on each particle and its motion. Typically, a particle can be classified into two types of motion in DEM: Translation and rotation. Momentum and energy of particles are exchanged during collisions with their neighbors or a boundary wall (contact forces), and particle-fluid interactions, as well as gravity. Through the application of Newton's second law of motion, we can determine the trajectory of each i-particle (including its acceleration, velocity, and position) from the following equations:

Discrete Element Method (DEM)
The Discrete Element Method (DEM) was initially suggested by Cundall and Strack [27] to model the mechanical behavior of granular flows and to simulate the forces acting on each particle and its motion. Typically, a particle can be classified into two types of motion in DEM: Translation and rotation. Momentum and energy of particles are exchanged during collisions with their neighbors or a boundary wall (contact forces), and particle-fluid interactions, as well as gravity. Through the application of Newton's second law of motion, we can determine the trajectory of each i-particle (including its acceleration, velocity, and position) from the following equations: where m i = the mass of a particle i; → u i = the velocity of a particle; → g = Gravity acceleration; → f i,k = interaction force between particle i and particle k (contact force); → f i,f = interaction force between the particle i and the fluid; I = moment of inertia; → ω i = angular velocity; d i = diameter of the grain i; and → n i,k = directional contact = vector connecting the center of grains i and k. We use a contact force model based on the principle of spring-dashpot as well as suggestions of Hertz-Mindlin [28]. The contact force is obtained from a force analysis method; the stiffness and damping factors are analyzed in two directions: Orthogonal and tangent of the contact surface between the two grains ( Figure 2): ( Water 2019, 10, x FOR PEER REVIEW 4 of 21 where mi = the mass of a particle i; u ⃗ = the velocity of a particle; g ⃗ = Gravity acceleration; f ⃗ , = interaction force between particle i and particle k (contact force); f ⃗ , = interaction force between the particle i and the fluid ; I = moment of inertia; ω ⃗ = angular velocity; di = diameter of the grain i; and n ⃗ , = directional contact = vector connecting the center of grains i and k.
We use a contact force model based on the principle of spring-dashpot as well as suggestions of Hertz-Mindlin [28]. The contact force is obtained from a force analysis method; the stiffness and damping factors are analyzed in two directions: Orthogonal and tangent of the contact surface between the two grains ( Figure 2):  (3) (n) and (τ) are known as two components of contact force in normal and tangential directions; ki = stiffness of grain i; δi,k = the characteristic of the contact and displacement (also called the length of the springs in the two directions above); αi = damping coefficient; and Δui = relative velocity of grain at the moment of collision. Following Coulomb, the value of tangential friction is determined by the product of the friction coefficient μ and the orthogonal force component. In the nonlinear contact force, Hertz-Mindlin model, the tangential force component will increase until the ratio (f (τ) /f (n) ) reaches a value of μ, and it retains the maximum value until the particles are no longer in contact with each other. A detail of the force models, as well as the method for determining the relevant coefficients, can be found in Reference [28].
After calculating all forces acting on the sediment particles as well as the velocity and the position of the particle at a previous time step, we can determine the current velocity and position of grain by solving Equations (3) and (4). The grain size distribution, as well as the bed porosity for whole the domain, can be defined afterward. As a result, we can also estimate the exchange rate of the fine fraction between different bed layers.
The DEM simulations begin with defining the system geometry. This comprises boundary conditions, particle coordinates, and material properties by identifying the contact model parameters, such as the friction and stiffness coefficients. How loading or deforming occurs within the system can be determined by the user through adding loads, deformations, or settlements. The simulation begins as either a transient or dynamic analysis and runs until the completion of a defined number of time steps. An overlap check procedure starts after particles are inserted into the simulation box, which is conducted based on the geometry and coordinates of the particles. Upon the simulation of motions starting, particles that physically encounter each other are detected, and the contact forces are then calculated at each time step. The magnitude of particle forces is related to the distance between the each of the contacting particles. From this data, the resultant force including, body forces, external forces, and moment acting on each particle can be calculated.
Moreover, two sets of equations for the dynamic equilibrium of the particles are computed in the case when particle rotation is blocked. Each particle translational movement is derived from the resultant applied force and each particle rotational movement is formulated from the resultant (n) and (τ) are known as two components of contact force in normal and tangential directions; k i = stiffness of grain i; δ i,k = the characteristic of the contact and displacement (also called the length of the springs in the two directions above); α i = damping coefficient; and ∆u i = relative velocity of grain at the moment of collision. Following Coulomb, the value of tangential friction is determined by the product of the friction coefficient µ and the orthogonal force component. In the nonlinear contact force, Hertz-Mindlin model, the tangential force component will increase until the ratio (f (τ) /f (n) ) reaches a value of µ, and it retains the maximum value until the particles are no longer in contact with each other. A detail of the force models, as well as the method for determining the relevant coefficients, can be found in Reference [28].
After calculating all forces acting on the sediment particles as well as the velocity and the position of the particle at a previous time step, we can determine the current velocity and position of grain by solving Equations (3) and (4). The grain size distribution, as well as the bed porosity for whole the domain, can be defined afterward. As a result, we can also estimate the exchange rate of the fine fraction between different bed layers.
The DEM simulations begin with defining the system geometry. This comprises boundary conditions, particle coordinates, and material properties by identifying the contact model parameters, such as the friction and stiffness coefficients. How loading or deforming occurs within the system can be determined by the user through adding loads, deformations, or settlements. The simulation begins as either a transient or dynamic analysis and runs until the completion of a defined number of time steps. An overlap check procedure starts after particles are inserted into the simulation box, which is conducted based on the geometry and coordinates of the particles. Upon the simulation of motions starting, particles that physically encounter each other are detected, and the contact forces are then calculated at each time step. The magnitude of particle forces is related to the distance between the each of the contacting particles. From this data, the resultant force including, body forces, external forces, and moment acting on each particle can be calculated.
Moreover, two sets of equations for the dynamic equilibrium of the particles are computed in the case when particle rotation is blocked. Each particle translational movement is derived from the resultant applied force and each particle rotational movement is formulated from the resultant applied moment. By knowing the inertia of the particles, particle translational, and rotational accelerations can be calculated. After new contract forces are determined, the particle positions and orientations are updated and ready for the next time step and will be repeated for all time steps. While this system seems to respond in an almost static manner, the Discrete Element Method is a transient or dynamic analysis. Figure 3 shows the calculation series that occur within a given time step. Particle velocities and incremental displacements are the first to be calculated. Here, the equilibrium of each particle in the sequence is considered. In the second series of calculations, upon the system geometry being updated, the forces at each contact in the whole system are then calculated. The particle rotational moment is produced from the normal contact force, as well as the tangential component of the contact force. As the output of these calculated moments and forces, the new particle position is generated for the next time step, and the series of calculation begins again. For every particle-based DEM simulation, the following fundamental assumptions are accepted. The first consideration is that particles are rigid, each possessing a finite inertia that can be described analytically. Moreover, the particles can translate and rotate independently of each other. The detections of new particle contacts are automatically completed by a geometry check algorithm. Physical contacts of particles normally happen over an infinitesimally small area based on the allowed overlapping and consists of only two particles. Particles that interact in DEM simulations are authorized to overlap slightly at the contact point, where the magnitude of the overlap is required to be small. The compressive inter-particle forces can be calculated from the particles overlapping value. Tensile and compression forces can be transferred at particle contact points and is normal in the direction of contact, as well as a tangential force orthogonal to the normal contact force. Furthermore, there is distance between two separating particles where the tensile inter-particle forces are calculated. When particles collide, this force is its maximum value, and then the particles move away from each other, which also means that the contact area diminishing to zero and is no longer used in contact force calculations. The last key assumption is that clusters of the rigid base particles can be used to represent a single particle. A measurable deformation of the composite particles is caused by the relative motion of the base particles within the cluster. These particle agglomerates may also be rigid themselves.
Water 2019, 10, x FOR PEER REVIEW 5 of 21 applied moment. By knowing the inertia of the particles, particle translational, and rotational accelerations can be calculated. After new contract forces are determined, the particle positions and orientations are updated and ready for the next time step and will be repeated for all time steps. While this system seems to respond in an almost static manner, the Discrete Element Method is a transient or dynamic analysis. Figure 3 shows the calculation series that occur within a given time step. Particle velocities and incremental displacements are the first to be calculated. Here, the equilibrium of each particle in the sequence is considered. In the second series of calculations, upon the system geometry being updated, the forces at each contact in the whole system are then calculated. The particle rotational moment is produced from the normal contact force, as well as the tangential component of the contact force. As the output of these calculated moments and forces, the new particle position is generated for the next time step, and the series of calculation begins again. For every particle-based DEM simulation, the following fundamental assumptions are accepted. The first consideration is that particles are rigid, each possessing a finite inertia that can be described analytically. Moreover, the particles can translate and rotate independently of each other. The detections of new particle contacts are automatically completed by a geometry check algorithm. Physical contacts of particles normally happen over an infinitesimally small area based on the allowed overlapping and consists of only two particles. Particles that interact in DEM simulations are authorized to overlap slightly at the contact point, where the magnitude of the overlap is required to be small. The compressive inter-particle forces can be calculated from the particles overlapping value. Tensile and compression forces can be transferred at particle contact points and is normal in the direction of contact, as well as a tangential force orthogonal to the normal contact force. Furthermore, there is distance between two separating particles where the tensile inter-particle forces are calculated. When particles collide, this force is its maximum value, and then the particles move away from each other, which also means that the contact area diminishing to zero and is no longer used in contact force calculations. The last key assumption is that clusters of the rigid base particles can be used to represent a single particle. A measurable deformation of the composite particles is caused by the relative motion of the base particles within the cluster. These particle agglomerates may also be rigid themselves. In this study, we used the open Source Discrete Element Method Particle Simulation (LIGGGHTS) implemented a new Hertz-Mindlin granular contact model [28,30,31], where grains are In this study, we used the open Source Discrete Element Method Particle Simulation (LIGGGHTS) implemented a new Hertz-Mindlin granular contact model [28,30,31], where grains are modeled as compressible spheres with a diameter d that interacts when in contact via the Hertz-Mindlin model [28,30,31]. An algorithm was developed to calculate grain size distribution and porosity from the calculated results of location and diameter of grains.
Defining a simulation time step is one of many essential steps in setting up the DEM. Sufficiently short time steps ensure stability of the system and enable stimulation of the real processes. According to Johnson [27,28], disturbances that occur during motion of particles in a granular system propagate following the Rayleigh waves form along the surface of solid. The simulation time step is included in the Rayleigh time, which is the time the energy wave takes to transverse the smallest element in the system. The simulation time step should be small enough so that any disturbance of a particle's motion only propagates to its nearest neighbors. Velocity and acceleration are assumed to be constant during the time step. Moreover, the time step duration should be smaller than the critical time increment evaluated from theory. Several equations have been proposed for calculating a critical time step [27]. In this study, we applied a time step of 0.00001 s, which is smaller than 20 percent of the Rayleigh time.

Algorithms for Calculating Grain Size Distribution and Porosity of a Cross Section from DEM Results
The results obtained from LIGGGHTS, an open source Discrete Element Method particle simulation software, contains 3D locations and diameter of grains. To calculate the porosity and grain size distribution, a simple algorithm was developed. We used the K different planes with elevations z k (k = 0, . . . , K), which intersect the spherical grain matrix. The diameter of generated circle i (i = 1, . . . , n k ) is dependent on the spherical diameter and the relative position between the k-plane and grain i ( Figure 4).
Water 2019, 10, x FOR PEER REVIEW 6 of 21 modeled as compressible spheres with a diameter d that interacts when in contact via the Hertz-Mindlin model [28,30,31]. An algorithm was developed to calculate grain size distribution and porosity from the calculated results of location and diameter of grains. Defining a simulation time step is one of many essential steps in setting up the DEM. Sufficiently short time steps ensure stability of the system and enable stimulation of the real processes. According to Johnson [27,28], disturbances that occur during motion of particles in a granular system propagate following the Rayleigh waves form along the surface of solid. The simulation time step is included in the Rayleigh time, which is the time the energy wave takes to transverse the smallest element in the system. The simulation time step should be small enough so that any disturbance of a particle's motion only propagates to its nearest neighbors. Velocity and acceleration are assumed to be constant during the time step. Moreover, the time step duration should be smaller than the critical time increment evaluated from theory. Several equations have been proposed for calculating a critical time step [27]. In this study, we applied a time step of 0.00001 s, which is smaller than 20 percent of the Rayleigh time.

Algorithms for Calculating Grain Size Distribution and Porosity of A Cross Section from DEM Results
The results obtained from LIGGGHTS, an open source Discrete Element Method particle simulation software, contains 3D locations and diameter of grains. To calculate the porosity and grain size distribution, a simple algorithm was developed. We used the K different planes with elevations zk (k = 0, …, K), which intersect the spherical grain matrix. The diameter of generated circle i (i = 1, …, nk) is dependent on the spherical diameter and the relative position between the k-plane and grain i ( Figure 4). The diameter of each circle created by the intersection between plane k and grain ith is calculated as: Total solid area (As,k) of all nk grains in plane k is determined: The total area At is calculated based on the shape generated by the plane k cut across the grain matrix, whereby porosity of cross section k is calculated by the following equation: The diameter of each circle created by the intersection between plane k and grain ith is calculated as: Total solid area (A s,k ) of all n k grains in plane k is determined: The total area A t is calculated based on the shape generated by the plane k cut across the grain matrix, whereby porosity of cross section k is calculated by the following equation: To calculate the grain size distribution, the grains in cross-section k are divided into m k size fractions with characteristic grain size D j (j = 1, . . . , m k ) and D j < D j+1 , then the area of each fraction is calculated by The fraction of class j in cross-section k is calculated by the following equation:

Feed Forward Neural Network (FNN)
Artificial Neural Network (ANN) is a general term encompassing many different network architectures. A Feedforward Neural Network (FNN) is an artificial neural network where connections between nodes do not form a cycle [32]. FNN is the first and simplest type of artificial neural network developed. Information of an FNN travels in only one direction, forward, from input nodes, through hidden nodes, then to the output nodes. Further, the most widely used FNN is a multilayer perceptron (MLP). An MLP model contains several artificial neurons otherwise known as processing elements or nodes. A neuron is a mathematical expression that filters signals traveling through the net. An individual neuron receives its weighted inputs from the connected neurons of the previous layer, which are normally aggregated along with a bias unit. The bias unit is purposed to scale the input to a useful range to improve the convergence properties of the neural network. The combined summation is delivered through a transfer function to generate the neuron output. Weighted connections modify the output as it is passed to neurons in the next layer, where the process is repeated. The weight vectors that connect the different network nodes are discovered through the so-called error back-propagation method. During training, these parameters values are varied in order for the FNN output to align with the measured output of a known dataset [33,34]. Changing the connections' weights in the network, according to an error minimization criterion, achieves a trained response. Overfitting is avoided if a validation process is implemented during the training. When the network has been sufficiently trained to simulate the best response to input data, the network configuration is fixed and a test process is conducted to evaluate the performance of the FNN as a predictive tool [22].
In feed-forward networks ( Figure 5), messages are passed forward only. A network with L layers has a parameter and a differentiable function f (l) : R d 1 → R d 1 corresponding to the lth layer. Given an input x ∈ R d o , the network outputs: where each a (l) ∈ R d l is defined and is defined recursively from the base case a (0) x as follows: The training process minimizes a loss function l : R d L ×d L → R over labeled examples (x, y). The gradient of the squared loss on (x, y) with respect to W(L) is The form mirrors the delta rule because a (L) = f (L) W (L) a (L−1) T where a (L−1) does not involve W (L) . By defining the "error term", we can simplify Equation (14) as δ (L) a (L−1) T . Similarly, the gradient with respect to W (l) for l < L can be verified to be δ (l) a (l−1) T where Computing all gradients in a multi-layer network in this manner is commonly known as "backpropagation", which is just a special case of automatic differentiation. For concreteness, here is the backpropagation algorithm for an L-layer feedforward network with the squared loss Output: Feedforward phase Set a (0) ← x , and for l = 1, . . . , L compute: Backpropagation phase Set δ (L) ← a (L) − y f (L) z (L) , and for l = L−1, . . . , 1 compute: Set W (l) ← δ (l) a (l−1) T , and for l = 1, . . . , L.
Water 2019, 10, x FOR PEER REVIEW 8 of 21 The form mirrors the delta rule because a ( ) = f ( ) W ( ) a ( ) where a ( ) does not involve W ( ) . By defining the "error term", we can simplify Equation (14) as δ ( ) a ( ) . Similarly, the gradient with respect to W ( ) for l < L can be verified to be δ ( ) a ( ) where Computing all gradients in a multi-layer network in this manner is commonly known as "backpropagation", which is just a special case of automatic differentiation. For concreteness, here is the backpropagation algorithm for an L-layer feedforward network with the squared loss Input labeled example (x, y) = R × R parameters W ( ) Output: Feedforward phase Set a ( ) ← , and for l = 1, …, L compute: Backpropagation phase Set δ ( ) ← a ( ) − y ⨀f ( ) z ( ) , and for l = L−1, …, 1 compute: Set W ( ) ← δ ( ) a ( ) , and for l = 1, …, L. The optimization algorithm (or optimizer) is the main approach used for training a machine-learning model to minimize its error rate. There are two metrics to determine the efficacy of an optimizer: Speed of convergence (the process of reaching a global optimum for gradient descent); and generalization (the model's performance on new data). Popular algorithms such as Adaptive The optimization algorithm (or optimizer) is the main approach used for training a machine-learning model to minimize its error rate. There are two metrics to determine the efficacy of an optimizer: Speed of convergence (the process of reaching a global optimum for gradient descent); and generalization (the model's performance on new data). Popular algorithms such as Adaptive Moment Estimation (Adam) or Stochastic Gradient Descent (SGD) can capably cover one or the other metric. The Adam optimizer, presented by Kingma and Ba [35], is extensively used for deep learning models requiring first-order gradient-based descent with small memory and the ability to compute adaptive learning rates for different parameters [36]. This method is computationally efficient, easy to implement, and has proven to perform better than the RMSprop and Rprop optimizers [37]. Gradient rescaling is reliant on the magnitudes of parameter updates. The Adam optimizer does not require a stationary object and can work with more sparse gradients. We calculate the decaying averages of past and past squared gradients m t and v t , respectively, as follows: m t and v t are estimates of the first moment (the mean) and the second moment (the uncentered variance) of the gradients, respectively. m t and v t are initialized as vectors of 0, the authors of Adam noticed that they are biased towards zero, particularly during the initial time steps and during smaller decay rates (i.e., β 1 and β 2 are close to 1).
Bias-corrected first and second moment estimates are computed to counteract these biases: Parameters are then updated by: The default value in this study: β 1 = 0.9, β 2 = 0.999 and = 10 −8 with learning rate = 0.001. More detail about this method is available in Reference [35] To create the FNN architecture, one must first determine the number of layers of each type and the number of nodes in each of these layers. In an FNN, one or more hidden layer of sigmoid neurons are often found, subsequently followed by an output layer containing linear neurons or nodes. This is completed because by having multiple layers of neurons with nonlinear activation functions, it allows the network to learn nonlinear relationships that exist between input and output vectors [38]. There is debate surrounding if the performance of FNN improves from the addition of more hidden layers. It has been found that the instances where performance improves with a second (or third, etc.) hidden layer are very few. Thus, one hidden layer is claimed as adequate for most problems FNN aims to solve. Let it be known that the number of neurons in the input layer is equal to the number of input features in the data set. The output layer contains only a single node namely the bed porosity. The optimal size of the hidden layer is normally in the range of the size of the input and output layers [39]. In our study, an FNN, designed for the largest number of inputs 10, with three layers is created with number neural nodes in the input layer of 10, the hidden layer of 8, and the output layers of 1.

Evaluation of the Model Performance
The measures of correlation coefficient (R), root mean square error (RMSE), mean absolute error (MAE) are used to evaluate the performance of these two models, and are formulated in equations: where n is the number of measurements, x i is calculated value ith, y i represents the measured value ith.

Input Parameters for DEM
For modeling porosity, we considered the following two cases: Case-1: Numerical simulations are carried out for a cylinder container with a diameter of 0.18 m. The container is filled with uniform coarse grains with size D = 10 mm and uniform fine grains with size d = 0.14 mm, which are similar to the experiments conducted by McGeary [40]. We simulate bed porosity for 13 different grain size distributions by varying the fraction of fine grains. The average height of the sediment layer is about 0.50 m.
Case  [42]. These experiments were used to investigate fine infiltration processes into gravel-bed, which were conducted in a flume that was 26 m long, 0.9 m wide, and 0.33 m deep (10 cm thick layer of gravel) [42]. A very slow water flowrate was used in experiment No.2 and No.3 (total 7 experiments). Hence, the effects of water flow on infiltration and packing processes can be neglected in the numerical model. The grain and water densities with four model parameters used in the DEM can be found in Table 1. Table 1. Parameters for numerical simulation.

DEM Verification for Porosity
For 13 samples of Case-1, the diameter ratio of fine grain diameter to coarse grain diameter (d/D) is 0.14. In the DEM model, the vibration force with a wiggle amplitude of 0.001 m and a period of 0.06 s was applied to adjust the porosity. Figure 6a,b shows the 3D simulated results of 2/13 cylinder samples for the case without fine sediment and fine sediment fraction 0.4. Figure 6c shows the comparison between the measured and the calculated results for these 13 samples. We can realize that there is a difference between our simulation and measurement: Our simulation results (diamond markers)-the line shape is not as sharp as the measurement line (circle marker). This can be explained because the fine grain in our simulation did not completely fill in the void structure of coarse gravel due to the high friction coefficient of grain, and because of the short time and small amplitude of vibration.
Another cause for the tolerance of the porosity in the simulation is the convex and concave surface of the cylinder, which may lead to an increase of porosity and therefore an increase in the total volume of the cylinder. However, in general, the agreement between our simulation and McGeary [40] depicted that DEM performs well for porosity simulation (Table 2).
Water 2019, 10, x FOR PEER REVIEW 11 of 21 between our simulation and McGeary [40] depicted that DEM performs well for porosity simulation (Table 2).  Figure 7 shows the simulation result of the flume (Case-2) filled by uniform gravel with D = 8 mm. As can be seen in Figure 7b, our DEM model provided a value of 0.47 for bulk porosity along the flume, which is slightly smaller than the measured porosity in 'Run 1′ (0.48) conducted by Navaratnam et at. [41]. Figure 7b shows the comparison of porosity variation by the depth between our simulation and two experiments, Run 1 and Run 2 [41], where Run 2 has been carried out in a larger flume. At the surface and bottom of the flume, our porosity results agree with the measurement data. In the middle of flume elevation, our porosity results did not change as dramatically as the experimental results due to the absolute uniformity of diameters in simulation, which is very difficult to mimic in experiment. In addition, in the middle of flume elevation, our result is significantly smaller than the experimental porosity results. This can be explained by the fact that the uniform spherical grains used in our simulation, required for grain close packing are different from the irregular shapes of gravel used in the experiment for loose packing. In general, from the statistical performance (Table 2), DEM simulations are in a good agreement with the experimental porosity measurement in gravel-bed flume.    Figure 7 shows the simulation result of the flume (Case-2) filled by uniform gravel with D = 8 mm. As can be seen in Figure 7b, our DEM model provided a value of 0.47 for bulk porosity along the flume, which is slightly smaller than the measured porosity in 'Run 1 (0.48) conducted by Navaratnam et at. [41]. Figure 7b shows the comparison of porosity variation by the depth between our simulation and two experiments, Run 1 and Run 2 [41], where Run 2 has been carried out in a larger flume. At the surface and bottom of the flume, our porosity results agree with the measurement data. In the middle of flume elevation, our porosity results did not change as dramatically as the experimental results due to the absolute uniformity of diameters in simulation, which is very difficult to mimic in experiment. In addition, in the middle of flume elevation, our result is significantly smaller than the experimental porosity results. This can be explained by the fact that the uniform spherical grains used in our simulation, required for grain close packing are different from the irregular shapes of gravel used in the experiment for loose packing. In general, from the statistical performance (Table 2), DEM simulations are in a good agreement with the experimental porosity measurement in gravel-bed flume.
which is very difficult to mimic in experiment. In addition, in the middle of flume elevation, our result is significantly smaller than the experimental porosity results. This can be explained by the fact that the uniform spherical grains used in our simulation, required for grain close packing are different from the irregular shapes of gravel used in the experiment for loose packing. In general, from the statistical performance (Table 2), DEM simulations are in a good agreement with the experimental porosity measurement in gravel-bed flume.

DEM Verification for Infiltration
Numerical simulations were carried out for Case-3 and Case-4. We tested the DEM model with two small windows with an edge of 0.15 m. In the first 41,000 time-steps, 1878 gravel packings were generated and reached a stable state. In the case of bridging, from time-step 41,000th to 320,000th, the 120,000th fine grains sediment were inserted. The time it took for fine sediment to settle was 60,000 time-steps. In the case of percolation, from the time-step 41,000 to 480,000, fine sediment (255,976 grains) was fed into the generated gravel-bed. The time for fine sediment to settle was 120,000 time-steps. Real calculation times (in second) for each process extracted from the log file can be seen in the Table 3.  Figure 8 presents the structure of the gravel-bed as well as the distribution of fine sediments for Case-3 and Case-4 at the end of the simulation. Figure 8a shows the 3D structure of the gravel-bed and fine sediment distribution in the bridging case. The infiltration process was stopped when the top gravel layer was filled. Figure 8c shows exemplarily bed materials at the middle x cross-section at the end of the simulation, where the clogging of fine sediment occurs at the surface. The formed fine sediment layer prevented the upper sediment from filling to the sublayer, where almost all void space remained empty. As can be seen in Figure 8(c1), although the diameters of fine particles are significantly smaller than the void space, they connect to build the 'bridge' across gravel called 'cake filtration', which depends on the size ratio of gravel and the vertical fine sediment rate [43][44][45]. Figure 8b shows the 3D structure of a gravel-bed and fine sediment distribution in the percolation case. The ratio of mean diameter gravel and fine sediment (D m /d m ) in our simulation (52,69) is in the range of percolation . In this ratio, fine particles are easy to infill to gravel, a fact consistent with what has been found in previous studies [46]. Figure 8d shows the middle x axis cross-section with most of its void space filled by sediment, however, not all void space was entirely filled. In the bottom of the simulation domain (Figure 8(d1)), fine sediment could not move down because of the bottom walls effect, leading to a sudden increase of fine fraction near the flume bed. This phenomenon usually occurred in flume experiments with gravel and fine sediment [10,42,47]. While there are some limitations in the time and scale of the simulations, it can be said that DEM is suitable for simulating the realistic 3D structure of fine sediment and gravel.  Figure 9 shows the comparison between the predicted results of fine sediment distribution and Gibson's [42] experiments in the bridging and percolation cases. In the top layer-bridging case, the fine sediment fraction reached its highest value (0.6) and values decrease with depth and a final large increase at the bottom. Opposition results found in the percolation case in Figure 9b presents a larger amount of fine sediment that was stored in the sublayer (average fine fraction 0.22) and at the surface layer (average fine fraction 0.19). From the bottom, we can observe a wave-form of fine fraction variation due to bottom wall effect and the interactions between particles. The amplitude of the wave is reduced with the elevation because of the influence of the wall effect, resulting in a reduction in chaotically stacked particles in the upper layer. To evaluate the performance of the model, the correlation coefficient (R), root mean square error (RMSE), and mean absolute error (MAE) are used-details of these equations are introduced in the next part. We reduced the resolution of simulated results from 500 points to 10 points because of the low resolution of the experiment results as well as the collected averaged seven measurements data at six different flume positions (5.5, 7.8, 9.8, 12.5, 16.5, 18.5 m, and 'the still water') [42]. It needs to be emphasized that the experiments and measurements have been conducted in a large flume with no sediment transport, while due to high computational requirements, our model only considered a small window, 0.15 m wide and 0.5 m long with quiescent water. There are small differences between our results and measurement data of Reference [42] due to the rescaling of the experiment. However, we obtained a good agreement between the experimental and numerical results ( Table 4). The validated results of  Figure 9 shows the comparison between the predicted results of fine sediment distribution and Gibson's [42] experiments in the bridging and percolation cases. In the top layer-bridging case, the fine sediment fraction reached its highest value (0.6) and values decrease with depth and a final large increase at the bottom. Opposition results found in the percolation case in Figure 9b presents a larger amount of fine sediment that was stored in the sublayer (average fine fraction 0.22) and at the surface layer (average fine fraction 0.19). From the bottom, we can observe a wave-form of fine fraction variation due to bottom wall effect and the interactions between particles. The amplitude of the wave is reduced with the elevation because of the influence of the wall effect, resulting in a reduction in chaotically stacked particles in the upper layer. To evaluate the performance of the model, the correlation coefficient (R), root mean square error (RMSE), and mean absolute error (MAE) are used-details of these equations are introduced in the next part. We reduced the resolution of simulated results from 500 points to 10 points because of the low resolution of the experiment results as well as the collected averaged seven measurements data at six different flume positions (5.5, 7.8, 9.8, 12.5, 16.5, 18.5 m, and 'the still water') [42]. It needs to be emphasized that the experiments and measurements have been conducted in a large flume with no sediment transport, while due to high computational requirements, our model only considered a small window, 0.15 m wide and 0.5 m long with quiescent water. There are small differences between our results and measurement data of Reference [42] due to the rescaling of the experiment. However, we obtained a good agreement between the experimental and numerical results ( Table 4). The validated results of the DEM method for simulation with fine sediment infiltration into gravel-bed are used to generate data for the FNN model, which is introduced in the next part. the DEM method for simulation with fine sediment infiltration into gravel-bed are used to generate data for the FNN model, which is introduced in the next part.
Data-classification-2: The mixture is characterized by two grain-sizes. Inputs parameters included location of sample (l), fraction of fine grain (with d < 2 mm), and fraction of coarse grain (with D ≥ 2 mm).
Based on the DEM results, grain size distribution and porosity data were generated for 500 different cross sections along the depth (z-direction) and for 800 different cross sections along the flume (x-direction). That means, we created totally 8 datasets: the 1st dataset using Data-classification-1 for bridging case and in z-direction (called Dataset-1), the 2nd dataset using Data-classification-2 for bridging case and in z-direction (called Dataset-2), the 3rd dataset using Data-classification-1 for percolation case and in z-direction (called Dataset-3), the 4th dataset using Data-classification-2 for percolation case and in z-direction (called Dataset-4),the 5th dataset using Data-classification-1 for bridging case and in x-direction (called Dataset-5), the 6th dataset using Data-classification-2 for bridging case and in x-direction (called Dataset-6), the 7th dataset using Data-classification-1 for percolation case and in x-direction (called Dataset-7), and the 8th dataset using Data-classification-2 for percolation case and in x-direction (called Dataset-8). These datasets
Data-classification-2: The mixture is characterized by two grain-sizes. Inputs parameters included location of sample (l), fraction of fine grain (with d < 2 mm), and fraction of coarse grain (with D ≥ 2 mm).
Based on the DEM results, grain size distribution and porosity data were generated for 500 different cross sections along the depth (z-direction) and for 800 different cross sections along the flume (x-direction). That means, we created totally 8 datasets: the 1st dataset using Data-classification-1 for bridging case and in z-direction (called Dataset-1), the 2nd dataset using Data-classification-2 for bridging case and in z-direction (called Dataset-2), the 3rd dataset using Data-classification-1 for percolation case and in z-direction (called Dataset-3), the 4th dataset using Data-classification-2 for percolation case and in z-direction (called Dataset-4),the 5th dataset using Data-classification-1 for bridging case and in x-direction (called Dataset-5), the 6th dataset using Data-classification-2 for bridging case and in x-direction (called Dataset-6), the 7th dataset using Data-classification-1 for percolation case and in x-direction (called Dataset-7), and the 8th dataset using Data-classification-2 for percolation case and in x-direction (called Dataset-8). These datasets contain also the information of cross section location.
Each dataset in x-direction is randomly divided into two subsets of data, namely 80% (400 data) for training and 20% (100 data) for testing purposes. Similarly, we randomly split 800 samples in x-direction of each dataset into two subsets: 80% (640 data) for training and 20% (160 data) testing purposes. Figure 10 shows the exemplary cumulative distribution at 10 different cross sections and grain distribution at cross-section 480th in the z-direction. x-direction of each dataset into two subsets: 80% (640 data) for training and 20% (160 data) testing purposes. Figure 10 shows the exemplary cumulative distribution at 10 different cross sections and grain distribution at cross-section 480th in the z-direction.

Porosity Prediction Based on FNN Model
Porosity depends on pressure and grain size distribution. In our DEM model, the effects of pressure on porosity, the consequence of forces acting on the grain matrix included gravity, buoyancy, grain friction, and contact force were considered. The influence of these factors contributed to the final location of grain obtained from DEM simulations. Grain diameter and the fraction of each size class are used to represent the characteristic of the grain size distribution.
The statistical indices (R, RMSE, and MAE) of the FNN model performance for four porosity predictions in the z-direction (Datasets 1, 2, 3, and 4) are presented in Table 5. In the bridging case, the prediction based on Data-classification-2 (Dataset-2) performed significantly better than Data-classification-1 (Dataset-1). Similarly, in the percolation case, the prediction using Dataset-4 is slightly better than Dataset-3. The prediction for percolation case achieved higher quality of result than the bridging case. It can be explained that in the bridging case, the fine fraction stored in void space of gravel-bed is smaller than in the percolation case. This lead to the prediction model being able to easily capture the effect fine sediment has on porosity output results. Furthermore, in the percolation case, the distinction between the coarse and fine groups is clear because of the large ratio of coarse gravel to fine sediment. As a result, the errors in the predicted results also reached minimum values (RMSE = 0.005753, MAE = 0.003155) from the model using Data-classification-2 of percolation (Dataset-4).

Porosity Prediction Based on FNN Model
Porosity depends on pressure and grain size distribution. In our DEM model, the effects of pressure on porosity, the consequence of forces acting on the grain matrix included gravity, buoyancy, grain friction, and contact force were considered. The influence of these factors contributed to the final location of grain obtained from DEM simulations. Grain diameter and the fraction of each size class are used to represent the characteristic of the grain size distribution.
The statistical indices (R, RMSE, and MAE) of the FNN model performance for four porosity predictions in the z-direction (Datasets 1, 2, 3, and 4) are presented in Table 5. In the bridging case, the prediction based on Data-classification-2 (Dataset-2) performed significantly better than Data-classification-1 (Dataset-1). Similarly, in the percolation case, the prediction using Dataset-4 is slightly better than Dataset-3. The prediction for percolation case achieved higher quality of result than the bridging case. It can be explained that in the bridging case, the fine fraction stored in void space of gravel-bed is smaller than in the percolation case. This lead to the prediction model being able to easily capture the effect fine sediment has on porosity output results. Furthermore, in the percolation case, the distinction between the coarse and fine groups is clear because of the large ratio of coarse gravel to fine sediment. As a result, the errors in the predicted results also reached minimum values (RMSE = 0.005753, MAE = 0.003155) from the model using Data-classification-2 of percolation (Dataset-4). Figure 11 shows the comparison between the FNN based porosity and the data obtained from DEM. In Figure 11a,d, the predicted porosity variation along the depth is compared with porosity based on Dataset-1 and Dataset-2. The dark magenta dotted line represents the DEM based data, where the minimum porosity reaches the value of 0.24 at the surface layer and increases with depth to reach a maximum value 0.6 near the bottom. Near the box bottom, porosity fluctuated widely. Gravel packing created two clear connection areas: a coarse gravel layer with the bottom box and with the upper layer that suddenly increased void spaces seen at z = 0.00 m and z = 0.012 m in Figure 10a,d. This also confirms the problem with laboratory porosity experiments, specifically how the disturbance of packing near the wall of the container causes pores size near the walls to be consistently larger than near the center of the container (also discussed in Reference [48]).   Figure 12 shows the performance of FNN models using four datasets along the horizontal x-direction for bridging and percolation cases in comparison with the DEM based data. The dark magenta dot-dash lines show the DEM based porosity along the horizontal direction. The porosity values changed from 0.37 to 0.51 in the bridging case-upper panel (Figure 12a) and from 0.22 to 0.40 for percolation-lower panel (Figure 12b). The average oscillation amplitude of porosity for these four cases (0.18) near the sidewall varies significantly more than inside the domain (0.06). This  Figure 11b,c,e,f show the performance of the FNN model and the scatter of porosity for the test dataset. As can be seen in Figure 11b, the FNN prediction using Dataset-1 in bridging case overestimated significantly the high peak (0.58, 0.01), (0.47, 0.018), and (0.43, 0.044) in comparison with the DEM based data. A light overestimation also occurs with the FNN prediction based on Dataset-3 in the percolation case (Figure 11d,e). A point worth noting is that both Dataset-1 and Dataset-3 used Data-classification-1 with nine sizes of grains. While, the performance of the FNN model is very good for Dataset-2 and Dataset-4, which are based on Data-classification-2 with two sizes of grains (Figure 11c,f). Overall, as can be seen in Figure 11, FNN models provide good results for porosity prediction. This suggests that the data-driven method based on the grain size distribution is suitable for porosity prediction along the depth in a gravel-bed. Figure 12 shows the performance of FNN models using four datasets along the horizontal x-direction for bridging and percolation cases in comparison with the DEM based data. The dark magenta dot-dash lines show the DEM based porosity along the horizontal direction. The porosity values changed from 0.37 to 0.51 in the bridging case-upper panel (Figure 12a) and from 0.22 to 0.40 for percolation-lower panel (Figure 12b). The average oscillation amplitude of porosity for these four cases (0.18) near the sidewall varies significantly more than inside the domain (0.06). This can be partly explained because, when the distance is equal to one medium radius (0.007 m) from the wall, the center of coarse particles stack along a vertical plane parallel to the wall, where the highest density of material is reached and reduces with distance from the center of the grain. The effect of the sidewall on the increasing porosity was also discussed in the previous experimental study [41,48,49].   In Figure 12a,b, the FNN prediction based on four Datasets (5, 6, 7, and 8) gave poor results at profiles, x = 0.010 and 0.140 and some other profiles. The model did not agree with the curve of DEM based porosity distribution (e.g., x = 0.095-0.135 m). The model tends to over fit a few points with sudden increases and decreases of fine sediment due to the two sidewalls. Inversely, inside the flume when the porosity did not fluctuate significantly, the model was found to under fit. However, in general, the performances of models are in good agreement with the DEM based data. Table 6 shows relatively good results of porosity prediction and in x-direction, FNNs using Dataset-6, and Dataset-8 perform better than using Dataset-5 and Dataset-7, respectively. Eight datasets were used to train the FNN networks. Regarding the influence of the grain classifications on the efficiency of FNN models, we observed slight differences between Data-classification-1 and Data-classification-2. As can be seen in the Tables 5 and 6, the statistical parameters of datasets based on Data-classification-2 are better than datasets based on Data-classification-2. Interestingly, with more detailed classification, the porosity prediction is not as good as in the less detailed classification. This suggest that, with the large size ratio of coarse gravel to fine sediment (D/d greater than 6.4), the detailed classification contained the little information groups (usually the coarse groups with the diameter larger than 2 mm) that may cause some inaccuracies in predicting porosity. This is consistent with the former study conducted by Bui et al. [7] in that the variation of porosity of gravel-bed is mainly caused by the variation in fine sediment rather than the effect of the rearrangement of coarse gravel. The redundancy of usefulness information increased the capability of FNN.
While the time required for training was significantly long (approximation 4 hours for one dataset), we have to train the model only once for each dataset. After the training process, we obtained an explicit relationship between porosity and the inputs, which can then be used to calculate the porosity with other datasets of inputs. Thus, the time needed to calculate porosity is reduced significantly by employing a FNN method in comparison with using sole DEM. It is needed to emphasize that a data-driven method can entirely replace the DEM in calculating porosity. Of cause, this replacement is strictly applied for the same variable range of the grain size distribution. Nonetheless, the reduced simulation time does not only save the computer resources, but also makes the connection between FNN and conventional hydro-morphodynamics model more robust.

Conclusions
The void of bed material plays an important role in fluvial geomorphology, exchange processes between river and groundwater and river ecosystem. Thus, predicting the variation of void space in gravel-beds plays a crucial role in eco-hydraulic management and fine sediment budget. In this research, we developed a model framework, which combined Discrete Element Method and Artificial Neural Network for porosity prediction. DEM realistically simulated porosity samples and fine sediment infiltration into 3D gravel-bed structure. An algorithm was developed to extract the simulated DEM results to calculate grain size distribution and porosity. Eight different datasets were generated based on the DEM results and applied to design several Artificial Neural Networks to find the relationship between grain size distribution and porosity. The results demonstrated that the combination of DEM and ANN is successful in simulating porosity and the infiltration process of fine sediment into the gravel-bed.
The validity of the presented model in the case of flowing water and sediment transport has not yet been verified, but it is believed that this model framework has a good performance for the analysis of bed and porosity changes. As a next step, to increase the quality of the training data, the effect that flow velocity has on sediment packing will be considered by coupling with OpenFOAM 5.0.0, open-source computational fluid dynamics (CFD) software, to study fluid-particles interactions. This model concept can be applied to calculate porosity change in a gravel-bed river for bed-porosity variation models as well as to define the exchange rate of fine sediment in gravel-bed rivers.
Author Contributions: V.H.B. designed the study, processed and analyzed the data, developed the models, interpreted the results and wrote the paper. The study has been carried out under the supervision of M.D.B. and P.R., who contributed to the model development stage with theoretical considerations and practical guidance, assisted in the interpretations and integration of the results and helped in preparation of this paper with proof reading and corrections.
Funding: This research was funded by the Vietnam International Education Development (VIED) Scholarship.