© B. Baldoni et al., Published by EDP Sciences, 2026

1 Introduction

Granular materials are ubiquitous in both natural and industrial contexts [1]. In environmental settings, they play a critical role in processes such as erosion and sedimentation [1,2], while in industrial applications, they are involved in transport and storage operations [3,4]. Understanding and predicting the behavior of granular flows is therefore of significant importance. However, studying these flows poses substantial challenges due to the complex interplay of numerous interactions and parameters that influence grain behavior, such as electrostatic forces, van der Waals interactions, humidity, and grain size [5–7]. To investigate granular flows, several modelling approaches can be employed. The most intuitive is the discrete approach, which involves simulating each individual particle within a granular assembly [8]. While this method offers detailed insights, it is computationally intensive and thus limited in its ability to model large-scale or realistic configurations. These challenges can be addressed by adopting a continuum approach, in which granular flows are modeled using modified Navier–Stokes equations, incorporating a constitutive rheology specifically designed to capture the behavior of granular media. Among these rheologies, the well-established μ(I) model [9] provides an effective framework for describing dry granular flows. It is based on the idea that, in systems dominated by frictional grain interactions, the granular material can be modeled as a viscoplastic fluid whose effective viscosity depends on a dynamic friction coefficient. The friction coefficient is expressed as a function of a dimensionless quantity known as the inertial number, I, defined as follows:

$$I=\frac{d\sqrt{2}{D}_2}{\sqrt{\left|p\right|/{\rho }_g}}\mathrm{,}$$(1)

where p is the local pressure in the granular media, d and ρ_g are respectively the radius and the density of the grain particles and $D_2=\sqrt{D:D}$ is the second invariant of the symmetric part of the velocity gradient $D=(\nabla u+\nabla u^{†})/2$. The friction coefficient denoted by μ(I) is expressed as:

$$\mu (I)={\mu }_s+\frac{{\mu }_2-{\mu }_s}{I_0/I+1}\mathrm{,}$$(2)

where μ_s, μ₂ and I₀ are material dependant parameters. While μ(I) provides a robust theoretical foundation, its predictive capabilities hinge on the accurate determination of these parameters, which are often difficult to measure directly in practice.

To address this limitation, we turn to Physics-Informed Neural Networks (PINNs) as a tool for identifying the unknown rheological parameters within the continuum model. PINNs [10] are a class of artificial neural networks that belong to the broader field of physics-informed machine learning [11]. The core principle of PINNs is to embed known physical laws into the neural network training process, thereby guiding the model toward physically consistent predictions. There are several ways to incorporate such physical constraints into a neural network architecture. According to [12], PINNs correspond to a “ANN–differential equation–regularizing” framework: they rely on conventional artificial neural networks (ANNs), and the governing physics is introduced through partial differential equations, whose residuals are minimized during training. PINNs have already shown their effectiveness in solving a wide range of scientific problems [13], particularly in reconstructing flow fields [14,15] or in identifying parameters [16].

In this work, a PINN model is used to reconstruct granular flows and identify the underlying rheological parameters governing their behavior, using a simple experimental configuration: the quasi-two-dimensional collapse of a granular column (Fig. 1).

This paper is organized as follow. Section 2 introduces the granular flow configuration studied in this work and presents the simulation data used to train the model. Section 3 describes the PINN architecture, detailing its structure, the incorporation of physical laws, and the formulation of the loss function. In Section 4, we present the results of the model, focusing on the identification of the rheological parameters μ_s and μ₂, as well as the reconstruction of the pressure field. Finally, Section 5 concludes the paper by summarizing the main outcomes and suggesting directions for future work.

Fig. 1 Collapse of a granular column: image from the laboratory of a quasi-2D configuration from [17].

2 The granular collapse

Among the granular flow configurations for which the μ(I) rheology has been shown to be applicable, the granular collapse is particularly attractive due to its experimental accessibility and simplicity. It constitutes a well-controlled flow configuration that can be reproduced in laboratory experiments [17,18] or simulated numerically [19]. It consists of an initially confined column of grains that is suddenly released and allowed to spread under gravity over a horizontal surface (Fig. 1).

In this study, the flow is modeled in a two-dimensional (2D) configuration (Fig. 2), as it is assumed to be uniform along the z direction. The training data for the PINN model were generated using the two-phase flow solver Basilisk [20], which employs a Volume-of-Fluid (VOF) method to track the interface between air and the granular material. A series of numerical simulations were conducted in which the rheological parameters μ_s and μ₂ were varied within prescribed ranges (${\mu }_s\in [\mathrm{0.2,0.6}]$ and ${\mu }_2\in [\mathrm{0.54,0.84}]$). The flow variables are computed on a hierarchical octree mesh and subsequently interpolated onto a uniform 128 × 128 grid for database generation. Each point in the grid that lies within the granular material, i.e. where the VOF function c is strictly equal to 1, is associated with a spatial location (x_i, y_i) and a time coordinate t_i, for which the velocity components (u, v) are known. The set of all such spatio-temporal points defines the velocity training dataset, denoted as $X^{\text{data}}=\{(x_i\mathrm{,}{y}_i\mathrm{,}{t}_i{)\}}_{1\le i\le N_{\text{data}}}$.

For the pressure field, only the values at the interface between the granular material and the surrounding air are retained. Interface points are identified using the VOF function, namely where c = 0.5. Such values would typically be accessible in an experimental setting. Pressure values are assigned exclusively to these interface points during training. This subset of points constitutes the training dataset for pressure, denoted as X^surf.

Fig. 2 Sketch of the 2D configuration of the granular column collapse at: a) the initial state; b) the final state.

3 The PINN model

The PINN model (Fig. 3) is built upon a standard artificial neural network (ANN) that takes as input two spatial coordinates (x, y) and a temporal coordinate t, and outputs the corresponding components of the velocity field (u, v) and the pressure p.

To incorporate physical constraints, the derivatives of the predicted fields are computed via automatic differentiation [21], and used to evaluate the residuals of the governing equations, the Navier–Stokes equations, and the continuity equation:

$${\rho }_g\left(\frac{\partial v}{\partial t}+v\cdot \nabla v\right)+\nabla p-\nabla \cdot \left(2{\eta }_{g}D\right)-{\rho }_{g}g=0,$$(3)

$$\nabla \cdot v=0,$$(4)

where g = (0,−1) denotes the reduced gravity field and η_g the equivalent viscosity given by the μ(I)-rheology:

$${\eta }_g=\frac{\mu (I)p}{\sqrt{2}{D}_2}\mathrm{.}$$(5)

Thus, unlike a conventional neural network, the loss function in a PINN does not rely solely on the discrepancies between the predicted and observed outputs (e₀). Instead, it incorporates both a data mismatch term and a physics-based residual term. The total loss function is defined as follows:

$${\text{L}}_{tot}={}_{}{\text{L}}_{\text{data}}+{\text{L}}_{\text{res}}\mathrm{,}$$(6)

The equation residuals are included in the loss term 𝐿_res, which penalizes deviations from the governing equations at a set of collocation points X^res. While these points can in principle be chosen arbitrarily, they are taken to coincide with the data points in this study, i.e., $X^{\text{res}}=X^{\text{data}}$. The loss term is defined as:

$${\text{L}}_{\text{res}}=\frac{1}{2\left|X^{\text{res}}\right|}{\displaystyle \sum _{x\in X^{\text{res}}}}\left(e_u{(x)}^2+e_v{(x)}^2+e_{\text{cont}}{(x)}^2\right)\mathrm{,}$$(7)

where e_u, e_v, and e_cont denote the residuals of the governing equations, i.e., the left-hand sides of equations (3) and (4).

The term 𝐿_data in equation (6) accounts for the difference between the network predictions and the available training data, following the structure of a standard supervised loss function:

$$\begin{array}{l}{\text{L}}_{\text{data}}=\frac{1}{2\left|X^{\text{data}}\right|}{\displaystyle \sum _{x\in X^{\text{data}}}}{\left|u_{\text{N}}(x)-u_x\right|}^2+\frac{1}{2\left|X^{\text{data}}\right|}{\displaystyle \sum _{x\in X^{\text{data}}}}{\left|v_{\text{N}}(x)-v_x\right|}^2\\ +\frac{1}{2\left|X^{\text{surf}}\right|}{\displaystyle \sum _{x\in X^{\text{surf}}}}{\left|p_{\text{surf}\mathrm{,}\text{N}}(x)-p_{\text{surf}\mathrm{,}x}\right|}^2\mathrm{,}\end{array}$$(8)

where u_𝒩(x), v_𝒩(x) and p_{surf, 𝒩}(x) denotes the predictions at point x in X^data or X^surf, the training datasets defined in Section 2, and u_x, v_x and p_{surf, x} are the corresponding reference values at point x.

During training, the neural network parameters are updated to minimize the total loss function 𝐿_tot. In order to identify the rheological parameters μ_s and μ₂ of the μ(I) model, two additional trainable variables are introduced alongside the neural network’s weights and biases (Fig. 3). The overall optimization problem can therefore be formulated as follows:

$${\text{Θ}}^{*}\mathrm{,}{\mu }_s^{*}\mathrm{,}{\mu }_2^{*}=\underset{\Theta \mathrm{,}{\mu }_s\mathrm{,}{\mu }_2}{argmin}({\text{L}}_{\text{res}}+{}_{}{\text{L}}_{\text{data}}),$$(9)

where Θ^* denotes the optimized neural network’s parameters, and μ_s^*, μ₂^* the identified values of the rheological parameters.

Fig. 3 Representation of the PINN architecture used in this study. The neural network 𝒩 typically consists of three hidden layers with 250 neurons each. Automatic differentiation is used to compute derivatives of the predicted fields, which are then used to evaluate the residuals of the governing equations. Two additional trainable variables are initialized alongside the neural network’s parameters in order to identify the values of μs and μ2.

4 Results

In this section, we present the results obtained using the PINN model applied to the granular column collapse configuration. We assess the model’s ability to accurately identify the rheological parameters μ_s and μ₂ of the μ(I) rheology and to reconstruct the pressure field.

To this end, a set of neural networks was trained, with one network trained for each flow configuration. Training was performed for 2000 epochs on an Nvidia RTX A5000 GPU with a batch size of 512, while the learning rate was gradually decreased from 5 × 10⁻⁴ to 1 × 10⁻⁵. All networks shared the same architecture, consisting of three hidden layers of 250 neurons. The optimizer was Adam [22] and the swish activation function [23] was used in all layers except the output layer, which had no activation function.

To quantitatively evaluate the accuracy of the predicted fields and the identified parameters, we compute the relative L₂ norm between the PINN predictions and the reference data. For the pressure field p, the relative L₂ norm is defined as:

$${\epsilon }_p=100\times \sqrt{\frac{{\displaystyle \sum _{x\in X}}\left|\right|p_{\text{N}}(x)-p_x{\left|\right|}^2}{{\displaystyle \sum _{x\in X}}\left|\right|p_x{\left|\right|}^2}}\mathrm{,}$$(10)

where p_𝒩(x) is the value predicted by the model at point x, and p_x is the corresponding reference value. This error is expressed as a percentage to facilitate interpretation. The same metric is also employed to assess the accuracy of the identified scalar parameters, μ_s and μ₂.

4.1 Identification of μ_s and μ₂

In order to obtain the best possible parameter identification results, we investigated the influence of the choice of observation points X^data used to train the PINN. In particular, we tested different distributions of X^data to assess whether concentrating observations in regions where the flow is active and strain-rate gradients are large could improve the identification of the rheological parameters. Figure 4 shows the evolution of the identified parameters μ_s (left) and μ₂ (right) over the training epochs for several choices of observation point distributions. The points are either taken from the full uniform computational grid defined in Section 2 or sampled according to a probability law weighted by the maximum absolute value of D₂.

Although all sampling strategies ultimately lead to comparable final values of μ_s and μ₂, slight differences appear in the convergence behavior. In particular, μ₂ typically requires more iterations to converge, which makes it important to ensure the proper convergence of both parameters. In this context, the case using the full uniform grid exhibits the fastest and most stable convergence. In addition, using the full grid leads to a significantly better reconstruction of the pressure field, which is particularly important in the multitask setting considered here, where both parameter identification and pressure reconstruction must be accurately achieved. Consequently, the results presented in the following are obtained from models trained using observation points distributed over the full uniform grid.

Table 1 summarizes the results obtained from the PINN-based parameter identification procedure for various combinations of the true values of μ_s and μ₂. Overall, the model demonstrates robust performance in retrieving both parameters, with relative errors consistently below 2% for μ_s, and typically below 5% for μ₂. The identification of μ_s appears particularly robust, with errors rarely exceeding 1.5%, even across a broad range of values from 0.2 to 0.6. The identification of μ₂ shows slightly more variation, with errors ranging from approximately 1% to 7%, but remains well within acceptable bounds for practical applications. This confirms the ability of PINNs to accurately infer rheological properties from partial flow data, even when both parameters are learned simultaneously.

Fig. 4 Evolution of the identified rheological parameters μs (left) and μ2 (right) as a function of the training epochs for different numbers of observation points Xdata. Observation points are either taken from the full uniform computational grid (in blue) or sampled according to a probability distribution weighted by the maximum absolute value of D2. The dashed lines indicate the reference values of the identified parameters.

4.2 Reconstruction of the pressure field

The last column of Table 1 reports the relative L₂ norm for the reconstructed pressure field ε_p computed over the whole domain, for each test case. These errors are typically below 4%, indicating a high level of agreement between the PINN predictions and the reference pressure fields. The only exceptions occur for the largest values of μ_s and μ₂, where the error increases to 7–8%, which may be attributed to the reduced sensitivity of the flow dynamics to the pressure distribution in these regimes or to numerical stiffness effects. Despite this slight increase, the accuracy remains acceptable, and the results confirm that the pressure field can be reliably reconstructed from velocity and limited interface pressure data using the PINN framework.

Figure 5 illustrates the comparison between the pressure field obtained from the CFD simulation using Basilisk (a) and the one reconstructed by the PINN model (b) at time t = 4.13 for μ_s = 0.3 and μ₂ = 0.54. The two fields exhibit a good agreement, indicating that the PINN successfully captures the spatial structure of the pressure field. In panel (c), three pressure profiles along the x-direction are plotted for different fixed y-positions, highlighted by colored lines in panels (a) and (b). The profiles predicted by the PINN (dashed lines) closely follow those from the CFD data (solid lines). The overlap is particularly strong in the flowing regions of the granular material. Slight deviations are observed in the static zones, where the granular material does not flow. These discrepancies are expected, as the underlying rheological model and governing equations used by the PINN are not valid in these non-flowing regions.

Fig. 5 At t = 4.13, for μs = 0.3 and μ2 = 0.54: a) the pressure field from Basilisk; b) the pressure field predicted by the PINN model; c) comparaison of three pressure profiles along the x-direction, the corresponding y coordinantes for the profiles are indicated in a) and b).

5 Conclusion

We presented a first application of PINNs in the context of continuum modeling of granular materials, using the well-established μ(I)-rheology framework. By training the model on synthetic data of a simple granular column collapse, we showed that PINNs can accurately reconstruct the pressure field and identify key rheological parameters μ_s and μ₂ with high precision. A key advantage of this approach lies in its ability to infer hard-to-measure physical quantities using only limited, easily accessible data, such as velocity fields and surface pressure. Our results suggest that PINNs provide a promising, flexible tool for rheological characterization of granular materials, especially in contexts where experimental measurements are limited or costly. Future work will focus on extending this approach to real experimental data, as well as to more complex granular configurations.

Funding

A part of this study is conducted in the framework of the Agence Nationale de la Recherche (Project ‘Pinn’terfaces’ ANR-23-CE23-0020).

Conflicts of interest

The authors have nothing to declare.

Data availabily statement

Data will be made available on request.

Author contributions statement

Barbara Baldoni: Conceptualization, Data Curation, Investigation, Methodology, Software, Validation, Visualization, Writing - Original Draft. Mickaël Delcey: Methodology, Software. Yoann Cheny: Conceptualization, Funding acquisition, Methodology, Supervision, Writing - Original Draft, Writing & Review & Editing. Adrien Gans: Data Curation, Methodology, Writing & Review & Editing. Mathieu Jenny: Writing & Review & Editing. Sébastien Kiesgen de Richter: Writing & Review & Editing.

References

All Tables

All Figures

	Fig. 1 Collapse of a granular column: image from the laboratory of a quasi-2D configuration from [17].
In the text

	Fig. 2 Sketch of the 2D configuration of the granular column collapse at: a) the initial state; b) the final state.
In the text

Fig. 3 Representation of the PINN architecture used in this study. The neural network 𝒩 typically consists of three hidden layers with 250 neurons each. Automatic differentiation is used to compute derivatives of the predicted fields, which are then used to evaluate the residuals of the governing equations. Two additional trainable variables are initialized alongside the neural network’s parameters in order to identify the values of μs and μ2.

In the text

Fig. 4 Evolution of the identified rheological parameters μs (left) and μ2 (right) as a function of the training epochs for different numbers of observation points Xdata. Observation points are either taken from the full uniform computational grid (in blue) or sampled according to a probability distribution weighted by the maximum absolute value of D2. The dashed lines indicate the reference values of the identified parameters.

In the text

	Fig. 5 At t = 4.13, for μs = 0.3 and μ2 = 0.54: a) the pressure field from Basilisk; b) the pressure field predicted by the PINN model; c) comparaison of three pressure profiles along the x-direction, the corresponding y coordinantes for the profiles are indicated in a) and b).
In the text