<figure><img src="image_1.png" /></figure>

HAL
open science

# Inferring viscoplastic models from velocity fields: a physics-informed neural network approach

Martin Lardy, Sham Tlili, Simon Gsell

► To cite this version:

Martin Lardy, Sham Tlili, Simon Gsell. Inferring viscoplastic models from velocity fields: a physics-informed neural network approach. 2025. hal-05389416

HAL Id: hal-05389416

https://hal.science/hal-05389416v1

Preprint submitted on 30 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Graphical Abstract

## Inferring viscoplastic models from velocity fields: a physics-informed neural network approach

Martin Lardy, Sham Tlili, Simon Gsell

<figure><img src="image_3.png" /></figure>
Inferring viscoplastic models from velocity fields: a
physics-informed neural network approach

Martin Lardy<sup>a,b</sup>, Sham Tlili<sup>b</sup>, Simon Gsell<sup>a</sup>

<sup>a</sup>Aix Marseille Univ, CNRS, Centrale Med, IRPHE, Turing Center for Living Systems, Marseille, France

<sup>b</sup>Aix-Marseille Univ, CNRS, IBDM, Turing Center for Living Systems, Marseille, France

**Abstract**

Fluid-like materials are ubiquitous, spanning from living biological tissues to geological formations, and across scales ranging from micrometers to kilometers. Inferring their rheological properties remains a major challenge, particularly when traditional rheometry fails to capture their complex, three-dimensional, and often heterogeneous behavior. This difficulty is exacerbated by system size, boundary conditions, and other material-specific physical, chemical, or thermal constraints. In this work, we explore whether rheological laws can be inferred directly from flow observations. We propose a physics-informed neural network (PINN) framework designed to learn constitutive viscoplastic laws from velocity field data alone. Our method uses a neural network to interpolate the velocity field, enabling the computation of velocity gradients via automatic differentiation. These gradients are used to estimate the residuals of the governing conservation laws, which implicitly depend on the unknown rheology. We jointly optimize both the constitutive model and the velocity field representation by minimizing the physical residuals and discrepancies from observed data. We validate our approach on synthetic velocity fields generated from numerical simulations using Herschel-Bulkley, Carreau and Panastasiou models under various flow conditions. The algorithm reliably infers rheological parameters, even in the presence of significant noise. We analyze the dependence of inference performance on flow geometry and sampling, highlighting the importance of shear rate distribution in the dataset. Finally, we explore preliminary strategies for model-agnostic inference via embedded model selection, demonstrating the potential of PINNs for identifying the most suitable rheological law from candidate models. This study illustrates how machine learning, and PINNs in particular, can enhance our ability to probe the rheology of complex fluids using velocity field data alone—paving the way for new approaches in computational rheology and material characterization.

Keywords: Rheology inference, Physics-Informed Neural Networks, Machine learning

PACS: 0000, 1111

2000 MSC: 0000, 1111

# 1. Introduction

Many natural and industrial materials exhibit complex rheological behavior, where deformation responds nonlinearly to the stresses driving it. Accurately modeling these rheological properties is crucial for understanding and predicting how materials deform—ranging from
micrometer-scale processes such as tissue morphogenesis [1], to kilometer-scale geodynamic flows in the Earth's mantle [2, 3].

Traditional rheometry involves bringing a sample to the laboratory and subjecting it to controlled stress or strain to characterize its constitutive behavior. This technique has been key to understand the mechanics of complex materials such as granular media [4] and dense suspensions [5, 6]. However, this approach also has several limitations. First, it typically provides only one-dimensional, local measurements that may not generalize to describe large-scale or spatially heterogeneous phenomena involving non-local effects or anisotropy. Second, many materials are difficult or impossible to study using standard rheometers—due to size constraints (from biological tissues to planetary layers), interfacial tribology (e.g., soft materials that do not adhere well to solid surfaces), or extreme physical conditions (e.g., high temperatures in lava, low temperatures in snow).

An alternative to conventional rheometry is to infer rheological properties directly from the way a material flows in situ—by observing its deformation in complex, natural, or engineered configurations. Across scales, material deformation is governed by the same fundamental conservation laws for mass and momentum. This raises the question: Can we invert these equations using experimental flow data to recover the constitutive law of the material?

Machine learning techniques have recently been developed to tackle this challenge by identifying unknown terms in partial differential equations (PDEs) that describe physical systems [7, 8, 9, 10, 11, 12]. In fluid mechanics, such methods have been used to infer hidden features such as boundary conditions, external forces, or material properties, using approaches ranging from Bayesian inference [13, 14] to neural networks [15, 16, 17].

Among these, Physics-Informed Neural Networks (PINNs) have emerged as a promising framework for solving forward and inverse problems in physics and engineering [18, 19]. In PINNs, physical variables such as velocity or stress are approximated using neural networks, and training is driven by minimizing a loss function that combines data fitting with residuals from governing PDEs. Automatic differentiation is used to compute spatial and temporal derivatives efficiently [20], enabling the enforcement of physical constraints during learning. This framework naturally accommodates inverse problems, allowing unknown parameters—including those in constitutive laws—to be inferred as part of the optimization. PINNs have been successfully applied to a variety of fluid mechanics problems [21], including the recovery of missing flow data [22, 23], and inference of power-law [24] or viscoelastic [25] properties from velocity fields, as well as viscoelastoplastic laws from rheometry data [26, 27].

In this work, we address the problem of inferring viscoplastic constitutive laws directly from flow velocity data using a PINN-based approach. Our algorithm employs a single neural network to represent both the velocity and pressure fields. We apply it to synthetic velocity fields generated from simulations of steady flows in confined geometries, inspired by microfluidic setups used to probe the rheology of foams [28], emulsions [29], and living tissues [30, 31]. These flows are computed using established viscoplastic models, including Herschel-Bulkley, Carreau, and Papanastasiou models.

We first demonstrate that, given prior knowledge of the correct constitutive model, the algorithm can accurately infer its parameters using only a single velocity snapshot—also reconstructing the corresponding pressure field. Importantly, the method remains robust in the presence of artificial noise, which is promising for application to experimental data.
We systematically analyze how inference performance depends on the structure of the input data, including the number and spatial distribution of velocity measurements. Our results highlight that a broad distribution of shear rates in the sampled data is critical for successful parameter recovery, and we propose a shear-rate-biased sampling strategy to improve performance. Finally, we present a preliminary strategy for model selection, enabling our algorithm to identify the most appropriate rheological model from a set of candidates before optimizing the associated parameters. Altogether, we believe our work lays a foundation for automated rheology inference from velocity fields, offering a data-driven alternative to classical rheometry that may later be applied to complex or experimentally constrained systems.

## 2. Rheology inference approach

### 2.1. Inverse problem formulation

We focus on two-dimensional incompressible fluid flows in creeping regime, governed by the Stokes equations:

$$ \nabla \cdot \mathbf{U} = 0, \qquad (1a) $$

$$ \nabla P = \nabla \cdot [2\eta S(\mathbf{U})], \qquad (1b) $$

with $\mathbf{U} = (u, v)$ the velocity field, $P$ the pressure and $S(\mathbf{U}) = (\nabla \mathbf{U} + (\nabla \mathbf{U})^T)/2$ the shear-rate tensor. For a viscoplastic fluid, the effective viscosity $\eta$ depends on the local shear rate $\dot{\gamma} = \sqrt{2S : \bar{S}}$. In the following, we aim at inverting equations (1) for a prescribed velocity field $\mathbf{U}$, i.e. finding the viscosity function $\eta$ and the pressure field $P$ such that equations (1a-b) are satisfied. In order to avoid the trivial solution $\eta = 0$ and $P = \text{const.}$, we introduce a characteristic viscosity $\eta_c \neq 0$ such that $\eta = \eta_c\eta^*$, with $\eta^*$ a non-dimensional viscosity. The other flow quantities are also made non-dimensional as follows:

$$ S^* = \frac{S}{\dot{\gamma}_c}, \quad U^* = \frac{U}{U_0}, \quad \mathbf{X}^* = \frac{\mathbf{X}}{D}, \quad P^* = \frac{P}{\eta_c \dot{\gamma}_c}, \qquad (2) $$

with $\mathbf{X} = (x, y)$ the spatial coordinates, $U_0$ a characteristic velocity (for example the mean flow velocity), $D$ a characteristic length and $\dot{\gamma}_c = U_0/D$ the characteristic shear rate. The non-dimensional form of equations (1a-b) reads

$$ \nabla \cdot \mathbf{U}^* = 0, \qquad (3a) $$

$$ \nabla P^* = \nabla \cdot (2\eta^* S^*). \qquad (3b) $$

In the following, we introduce three different rheological models for $\eta(\dot{\gamma}_c)$, namely the Herschel-Bulkley, Carreau and Papanoustasiou [32] models.

### 2.2. Viscoplastic models

We present here the three rheological models considered in our work and their non-dimensional form. The models are summarized in table 1.

**Herschel-Bulkley model.** The Herschel-Bulkley model accounts for shear-thinning/thickening and yield-stress behaviors of the fluid. This model is popular to characterize viscoplastic
<table><caption>Table 1: Summary of the rheological models considered in this work. The non-dimensional rheological parameters $α_i$ are obtained by choosing the characteristic viscosities $η_c = kγ̇_c$ (Herschel-Bulkley), $η_c = η_0$ (Carreau) and $η_c = μ_0$ (Papanastasiou).</caption><thead><tr><th></th><th># of parameters</th><th>α<sub>0</sub></th><th>α<sub>1</sub></th><th>α<sub>2</sub></th></tr></thead><tbody><tr><td>Herschel-Bulkley (5)</td><td>2</td><td>σ<sub>0</sub>/kγ̇<sub>c</sub><sup>n</sup></td><td>n</td><td>-</td></tr><tr><td>Carreau (7)</td><td>3</td><td>η<sub>∞</sub>/η<sub>0</sub></td><td>λγ̇<sub>c</sub></td><td>n</td></tr><tr><td>Papanastasiou (9)</td><td>2</td><td>τ<sub>0</sub>/μ<sub>0</sub></td><td>mγ̇<sub>c</sub></td><td>-</td></tr></tbody></table>

properties, and has been successfully used to describe the steady rheology of soft glassy and jammed materials or polymeric gels [33]. The viscosity function reads

$$ \eta_{HB} = \frac{\sigma_0}{\dot{\gamma}} + k\dot{\gamma}^{n-1}, \qquad (4) $$

with $σ_0$ the yield stress, $k$ the fluid consistency and $n$ the power-law index. We choose $η_c = kγ̇_c^{n-1}$ as the characteristic viscosity, assuming $k \neq 0$. The non-dimensional viscosity becomes

$$ \eta_{HB}^* = \frac{\alpha_0}{\dot{\gamma}^*} + (\dot{\gamma}^*)^{\alpha_1 - 1}, \qquad (5) $$

with $α_0 = \sigma_0/(k\dot{\gamma}_c^n)$ and $α_1 = n$ the two non-dimensional rheological parameters to be inferred. Note that $α_0$ corresponds to the standard definition of the Bingham number for yield-stress fluids.

**Carreau model.** The Carreau model only accounts for shear-thinning / thickening effects, but the power-law behavior of the fluid is localized within a certain shear-rate range. Outside of the power-law region, the fluid recovers a Newtonian behavior. The viscosity reads

$$ \eta_C = \eta_\infty + (\eta_0 - \eta_\infty)(1 + (\lambda\dot{\gamma})^2)^{\frac{n-1}{2}}, \qquad (6) $$

where $η_0$ and $η_\infty$ are the limit viscosities at zero and infinite shear rates, respectively, $λ$ is the characteristic time and $n$ is the power-law index. The Carreau model thus involves one extra parameter compared to the Herschel-Bulkley model. We choose $η_c = η_0$ as the characteristic viscosity, leading to the following non-dimensional viscosity law:

$$ \eta_C^* = \alpha_0 + (1 - \alpha_0)(1 + (\alpha_1 \dot{\gamma}^*)^2)^{\frac{\alpha_2 - 1}{2}}, \qquad (7) $$

with $α_0 = η_∞/η_0$, $α_1 = λγ̇_c$ and $α_2 = n$ the three non-dimensional parameters to be inferred.

**Papanastasiou model.** This model is a regularized version of the Bingham model (i.e. Herschel-Bulkley model with $n = 1$), accounting for yield-stress behavior but avoiding the divergence of the viscosity at zero shear rate. The viscosity model reads

$$ \eta_P = (1 - e^{-m\dot{\gamma}}) \frac{\tau_0}{\dot{\gamma}} + \mu_0, \qquad (8) $$
<figure><img src="image_4.png" /><figcaption>Figure 1: Overview of our PINN rheology inference algorithm. We use a neural network $N_\theta$ to represent the flow quantities $U_\theta^*$ and $P_\theta^*$. The distance between the output velocity $U_\theta^*$ and the velocity data $U^*$ is quantified by the data loss $\mathcal{L}_{data}$ (11). We then use automatic differentiation to determine the spatial gradients of $U_\theta^*$ and $P_\theta^*$, that are then combined with the rheological model $\eta_\alpha^*$ to estimate the physical loss $\mathcal{L}_{phy}$ (12). The total loss consists of the weighted sum of $\mathcal{L}_{data}$ and $\mathcal{L}_{phy}$, and also includes an additional loss $\mathcal{L}_c$ (13) enforcing the positivity of the rheological parameters $\alpha$. The rheology inference is performed by minimizing the total loss over the physical/rheological parameters $\alpha$ and interpolation parameters $\theta$.</figcaption></figure>

with $\mu_0$ the limit viscosity at zero shear rate, $\tau_0$ the yield stress and $m$ the characteristic time. We choose $\eta_c = \mu_0$, giving the following non-dimensional viscosity:

$$ \eta_P^* = (1 - e^{-\alpha_1 \dot{\gamma}^*}) \frac{\alpha_0}{\dot{\gamma}^*} + 1, \qquad (9) $$

with $\alpha_0 = \frac{\tau_0}{\mu_0 \dot{\gamma}_c}$ and $\alpha_1 = m \dot{\gamma}_c$ the non-dimensional rheological parameters to be inferred.

## 2.3. Rheological inference approach

Our inference algorithm is summarized in figure 1. We aim at inverting equations (3a-b) by determining the viscosity function $\eta_\alpha^*$, with $\alpha = \{\alpha_0, \alpha_1, ...\}$ the parameters of one of the three rheological models above. We perform the inverse problem from a set of velocity samples $U_n = (u_n, v_n)$, $n \in \{1, ..., N\}$ with $N$ the total number of samples. Each sample is normalized as

$$ U_n^* = \frac{(u_n - \bar{u}, v_n - \bar{v})}{U_0}, \qquad (10) $$

where $\bar{u}$ and $\bar{v}$ are the mean values of the $u$ and $v$ fields. The position $X_n$ of each velocity sample is also normalized as $X_n^* = X_n/D$.

Assuming that the data $U_n^*$ might be noisy and sparse, we first perform an interpolation of $U_n^*$ with a field $U_\theta^*(X^*)$, where $\theta$ is a set of interpolation parameters. In the PINN approach, we use a neural network $N_\theta$ to represent the function $U_\theta^*(X^*)$, with $\theta$ the bias and weights of the network. In practice, the output of the network is not directly $U_\theta^*$ but the
stream function $\psi_{\theta}$; the velocity is then reconstructed as $u_{\theta} = \partial\psi_{\theta}/\partial y^*$ and $v_{\theta} = -\partial\psi_{\theta}/\partial x^*$. This allows us to enforce the incompressibility condition (3a) for the velocity field $U_{\theta}^*$. In a compressible scenario, then the output would directly be $u_{\theta}$ and $v_{\theta}$. We also include the unknown pressure $P_{\theta}^*$ as an output of the same network $N_{\theta}$. The neural network used here is a fully connected feedforward network composed of 4 hidden layers with 20 neurons each, with a tanh activation function.

In order to train the network $N_{\theta}$, we first introduce the data loss $\mathcal{L}_{data}$ which quantifies the distance between the interpolated velocity field and the data $U_n$,

$$ \mathcal{L}_{data} = \|U_{\theta}^* - U_{n}^*\|, \qquad (11) $$

with $\|\cdot\|$ the 2-norm squared ($\|\cdot\|_2^2 = \sum^N(.)^2$), with N the number of training points. Assuming a guess rheological model $\eta_{\alpha}^*$, we then define the physical loss $\mathcal{L}_{phy}$ based on the momentum balance residual (3b),

$$ \mathcal{L}_{phy} = \| -\kappa \nabla P_{\theta}^* + \nabla \cdot (2\eta_{\alpha}^* S_{\theta}^*) \|, \qquad (12) $$

where $\nabla P_{\theta}^*$ and $S_{\theta}^*$ are determined using automatic differentiation of the network $N_{\theta}$. Even though not necessary, in the following we choose to calculate $\mathcal{L}_{phy}$ at the same colocation points as for $\mathcal{L}_{data}$. In (12) the pressure field can be rescaled by any arbitrary prefactor $\kappa$. As the pressure field magnitude depends on the characteristic viscosity $\eta_c$, we have found that the inference is improved when rescaling the pressure using $\kappa \propto 1/\eta_c$. In practice, we use $\kappa = \dot{\gamma}_c^{n-1}$ for the Herschel-Bulkley model and $\kappa = 1$ for the Carreau and Papanastasiou models. In the case of the Herschel-Bulkley model, this allows to encode the rescaling of the pressure by the power-law index $n$, so that the neural network does not have to learn this rescaling each time the rheological model is updated. In order to enforce the positivity of the rheological parameters, we finally introduce a positivity loss

$$ \mathcal{L}_p = \sum_i \max(0, -\alpha_i). \qquad (13) $$

The total loss defining our regression problem finally reads

$$ \mathcal{L} = \mathcal{L}_{data} + W\mathcal{L}_{pde} + \mathcal{L}_{p}, \qquad (14) $$

where $W$ is a constant weight used to balance the data and physical losses. In the following, we automatically fix the value of $W$ at the beginning of each simulation, by computing the initial ratio between $\mathcal{L}_{data}$ and $\mathcal{L}_{phy}$ such that $[\log_{10}(\mathcal{L}_{data}/\mathcal{L}_{phy})] + 1 = w$ and $W = 10^{w-1}$. During the rheology inference process, we aim to minimize the total loss $\mathcal{L}$ by performing a descent along gradients $\partial\mathcal{L}/\partial\theta$ and $\partial\mathcal{L}/\partial\alpha$ with respect to the interpolation and physical/rheological parameters, respectively. We perform this gradient descent using the ADAM algorithm [34], with a learning rate equal to 0.0025. When the optimal parameters $\alpha$ and $\theta$ are found, the algorithm provides a rheological model $\eta_{\alpha}$ minimizing the physical residual for a velocity field $U_{\theta}$ as close as possible to the velocity data $U_n$. All the Python scripts used in this work are based on the Pytorch library and are available at https://github.com/martinLARD/Pinn.git [35].
## 2.4. Benchmarking approach

We test our rheology inference approach using a series of non-Newtonian flows generated by numerical simulations. The simulations are based on the lattice-Boltzmann and immersed boundary methods [36, 37]. We perform simulations in different configurations similar to those classically used experimentally to characterize foams [28], emulsions [29], or cellular tissues [30, 31], such as a straight channel, different types of non-uniform channels, and a channel including a circular obstacle. The simulations are carried out for different viscoplastic models and various values of the associated parameters. Our lattice-Boltzmann method is based on a generalized-Newtonian approach, i.e. viscoplastic behaviors are purely modeled by locally varying the effective fluid viscosity, as described in detail in [36]. The velocity fields resulting from the simulations are sampled uniformly or with a bias based on the shear-rate distribution (see Section 4). For each inference run, the sampling points remain fixed during the calculation. We quantify the inference accuracy by calculating the relative log-viscosity error between the true viscosity function $\eta^*$ and the inferred viscosity $\eta_\alpha^*$,

$$E_{\eta} = \frac{1}{M} \sum_{i=0}^{M} \frac{|\log(\eta_{i}^{*}) - \log(\eta_{\alpha,i}^{*}(\dot{\gamma}^{*}))|}{|\log(\eta_{i}^{*})|}, \quad (15)$$

where M is the number of $\eta^*$ points available.

# 3. Inference of rheological parameters

## 3.1. Overview of the inference process

Figure 2 shows a series of inference examples applied to Herschel-Bulkley flows in a uniform channel, a channel obstructed by a circular obstacle and a non-uniform channel. In all cases, we perform numerical simulations providing velocity fields (figure 2a) from which we uniformly sample a total of 800 velocity data points provided to the inference algorithm. We will further discuss the effect of the colocation points in section 3.1.1 and section 4. Each configuration leads to a distinct viscosity pattern (figure 2b) that needs to be inferred by the algorithm. We fix the number of epochs/iterations to $6 \times 10^5$.

In all cases, the learned rheological parameters $\alpha_0$ and $\alpha_1$ converge towards their true values (figure 2c) as the total loss decreases during the training process (figure 2d). However, the accuracy of the parameter inference depends on the relative influence of each specific parameter on the resulting viscosity field. We thus further evaluate the inference quality by reconstructing the viscosity law $\eta^*(\dot{\gamma})$ and the viscosity field $\eta(x)$ from the inferred parameters (figures 2e,f). In all cases, the reconstructed viscosity matches with the true one, even when small errors are observed on the rheological parameters (e.g. $\alpha_0$ for the flow past the obstacle).

We plot the output fields of the neural network and compare them to the true fields in figure 3. The overall pressure pattern is well predicted by the PINN algorithm, even though pressure boundary conditions are not provided during the inference process. Some discrepancies are however visible between the original and inferred pressure fields, especially near solid boundaries in the flow past the cylinder and in the wavy channel. In the uniform channel case, pressure error mostly occurs in the central plug region, where the fluid is unyielded and the effective viscosity is very high (figure 2b). Overall, these pressure errors
<figure><img src="image_5.png" /><figcaption>Figure 2: Overview of the inference process by inferring Herschel-Bulkley parameters in different configurations. (a,b) We perform simulations of Herschel-Bulkley flows with different values of the parameters $\alpha_0$ and $\alpha_1$ in a uniform channel (top), a channel obstructed by a circular obstacle (center) and a non-uniform channel (bottom). The velocity field (column (a)) is sampled as the input data for the inference algorithm, while the viscosity field (column (b)) is unknown during the inference. (c,d) As the total loss $\mathcal{L}$ decreases during learning, the algorithm converges to inferred values for $\alpha_0$ and $\alpha_1$ (solid lines) that are close the true values (dashed lines). (e,f) The viscosity functions reconstructed using the inferred values of $\alpha_0$ and $\alpha_1$ match the true viscosity functions, both as a function of the shear rate (column (e)) or in space (column (f)).</figcaption></figure>

however remain smaller than 10% during the inference (figure 3b). As expected, discrepancies are even smaller for the velocity fields (figure 3c,d), due to the minimization of $\mathcal{L}_{data}$. Again, largest errors are observed close to solid boundaries and within the unyielded regions, but overall the velocity error remains smaller than 10% for all configurations.

### 3.1.1. Effect of colocation points

We quantify the impact of the number of colocation points on the inference accuracy in figure 4. We perform a series of 20 inferences with either 100, 400, 800, or 1200 collocation points and $8 \times 10^5$ epochs, for (a) the Herschel-Bulkley and (b) Carreau flows in a uniform channel. We note that the inference accuracy is generally higher for the Carreau model. For both rheological models, the inference accuracy generally increases with the number of points, which can be noted both from the values of the inferred parameters (figure 4a,b) and from the overall viscosity law error (figure 4c). Interestingly, the inference variability also drastically decreases with the number of colocation points. This is particularly visible from the viscosity error in the case of the Herschel-Bulkley rheology (figure 4c). This can be interpreted as a decrease of the inference uncertainty as the amount of input information is increased. However, for 800 colocation points and more, the inference accuracy does not significantly decrease with the number of points anymore; it even starts to increase again in the case of the Carreau model. Similar detrimental effect of a too large amount of colocation
<figure><img src="image_6.png" /><figcaption>Figure 3: Velocity and pressure fields at the output of the neural network, for the three examples shown in figure (2). We introduce the notation $\tilde{A} = \frac{A-\min(A)}{\max(A)-\min(A)}$, for A a given spatial field. The figure shows Herschel-Bulkley examples with parameters ($\alpha_0, \alpha_1$), from top to bottom, equal to (1.2, 0.1), (1, 0.1) and (1, 0.1). (a) Pressure fields. The upper part of the plot shows the output of the neural network ($\tilde{P}_{\theta}^*$); the lower part shows the pressure obtained from numerical simulation ($P$) both normalized. (b) Absolute difference between both pressure fields with a factor $a$ before the normalized pressure tuned to minimize the absolute difference integrated on the entire channel. This parameter is from top to bottom equal to 0.94, 1 and 0.93. (c) Horizontal velocity fields. The upper part of the plot shows the output of the neural network; the lower part shows the velocity obtained from numerical simulation. (d) Absolute difference between both velocity fields, with $b$ from top to bottom equal to 1.07, 0.98 and 0.96.</figcaption></figure>

points has also been reported in previous works on PINNs [38].

## 3.2. Noise robustness

An important challenge for the inference algorithm is to be able to perform on noisy experimental data. To explore this aspect, we add artificial noise to the velocity data provided to our algorithm. As the nature of the noise that needs to be removed from experimental data is generally unknown, we consider two different types of artificial noise, namely a gaussian and a uniform noise. We perform series of inference runs with various noise amplitudes,
<figure><img src="image_7.png" /><figcaption>Figure 4: Impact of the number of colocation points on the inference accuracy. A series of 20 inferences are performed during 800 000 epochs, each conducted with either 100, 400, 800, or 1200 collocation points, for (a) the Herschel-Bulkley and (b) Carreau flows in a uniform channel. (a,b) Mean values of the inferred rheological parameters. The error bars indicate the standard error of the mean. (c) Mean viscosity error (15) for both rheological models. Error bars indicate the error propagated from the rheological parameter values.</figcaption></figure>

corresponding to various standard deviation values for the gaussian noise and various limit values for the uniform noise, both normalized by the mean flow velocity.

Figure 5 shows the inference results in the case of Herschel-Bulkley flow in a uniform channel. For both types of noise, the inference accuracy decreases with the noise amplitude (figure 5a,b). However, even with a 10% noise amplitude, the inferred parameters remain reasonably close to their true values, especially in the case of uniform noise. At the end of the inference process, the interpolated velocity at the output of the neural network is entirely denoised and matches the original velocity data (5(c,d)). This illustrates that the physical loss $\mathcal{L}_{phy}$ acts as a physics-based regularization function during the velocity interpolation.

Interestingly, the presence of noise does not increase the inference uncertainty (quantified by the error bars in figure 5a,b), but rather introduces a consistent bias towards a Newtonian rheology. This is further illustrated in figure 5e, which shows (i) the original velocity profile across the channel, (ii) the sampling points across the same cross-section after addition of a 20% Gaussian noise, and (iii) the inferred velocity at output of the neural network. The reconstructed velocity indeed lacks the non-newtonian features observed in the original data, especially the central plug typical of viscoplastic channel flows. However, it is also clear from the figure that noise amplitude is large compared to the velocity difference between the Newtonian and non-Newtonian flows; capturing such non-Newtonian effects in a very noisy environment thus appears to be challenging for the method.

## 4. Biased sampling in different geometries

Accurately inferring viscosity laws requires to sample data across a broad range of shear rates. However the shear rate distribution may vary from one geometry to another, implying a geometry-bias in the rheology inference. Here, we propose to control this bias by performing a spatially non-uniform sampling, taking into account the shear-rate distribution in the available data. We noticed that in the considered configurations, the shear rate tends to be naturally distributed towards lower values. A first simple approach in order to alter the shear-rate distribution consists in biasing it towards higher values. Given access to the velocity field, we first compute the corresponding shear rate field. By binning this field, we
<figure><img src="image_8.png" /><figcaption>Figure 5: Rheology inference on noisy velocity data. (a,b) Mean values of the infered Herschel-Bulkley parameters in a uniform channel, over a series of 20 inference runs, for (a) Gaussian and (b) uniform noise of various amplitudes. The error bars indicate the standard error of the mean. Snapshots below the figures show the noisy velocity fields. (c) Velocity field rectonstruction in the case of a 5% Gaussian noise (star symbol in (a)): (top) original velocity field, before addition of noise, (bottom) infered velocity field at the output of the neural network. (d) Absolute difference between the original and infered velocity fields. (e) Velocity profile along the vertical direction at the center of the channel in the case of a 20% Gaussian noise: original velocity profile, sample points after addition of noise and infered velocity profile.</figcaption></figure>

determine the distribution of shear rate values across the bins, which provides a distribution denoted by $\Gamma(X)$. We then construct a new sampling probability distribution based on this information:

$$P_s(X) = \frac{\Gamma(X)^c}{\sum_i [\Gamma(X_i)^c]}, \qquad (16)$$

with c a user-defined constant controlling the bias magnitude and ($X_i$) being the points of the field. Equation (16) provides a simple way to promote sampling in high-shear regions.

Before the training process, we sample velocity data according to the probability distribution $P_s$ (16). Figure 6 illustrates this biased sampling in the case of the flow past an obstacle. As the parameter c increases, the spatial distribution of the sampling points gets spatially non-uniform (figure 6a-c) and the statistical distribution of shear-rate values gets
<figure><img src="image_9.png" /><figcaption>Figure 6: Illustration of the shear-biased sampling in the case of the flow past an obstacle. (a-c) Sampling point position for different values of the bias parameter *c*: (a) uniform sampling (*c* = 0), (b) *c* = 1 and (c) *c* = 4. (d) Box plots showing the statistical shear-rate distribution in each case. In each plot the horizontal lines correspond to the first and third quartiles, used to compute the inter-quartile range (IQR). (e) IQR as a function of the bias parameter *c*. (f) Inference error (15) as a function of the bias parameter *c*.</figcaption></figure>

broader (figure 6d). We quantify this aspect using the inter-quartile range (IQR) of the shear-rate distribution. The IQR provides a robust measure of the diversity of information within the distribution, as it is insensitive to extreme values. The IQR increases as the sampling gets more biased (figure 6e), which appears to significantly improve the inference accuracy. We interpret this as the beneficial effect of improving the diversity of shear-rate values over which the viscosity law can be learned.

We apply the same sampling strategy to various cases in different geometries and summarize the results in figure 7. Even though variations are observed from one geometry to another, this confirms that increasing the IQR generally improves the inference quality. This also confirms that the inference process can be applied to a variety of geometries, on a broad range of rheological parameters values. The bias parameter *c* thus provides a useful way to tune the sampling process depending on the shear-rate distribution within the available data, in order to improve the inference quality.

## 5. Automatic rheology model selection

So far we have illustrated the capacity of our PINN algorithm to fit rheological parameters based on velocity data, provided that the correct rheological model is known. We now address the more challenging problem of rheological model selection: is it possible to select the most
<figure><img src="image_10.png" /><figcaption>Figure 7: Effect of the inter-quartile range (IRQ) of the shear-rate distribution on the inference accuracy, for different configurations and different rheological parameters. Solid lines connect points with identical geometry and rheology, but different IQR achieved through biased sampling. Each point is averaged over 20 inference runs performed with a fixed value of the bias parameter c. Horizontal error bars indicate the standard error on the mean IQR value, as data points are re-sampled for each inference run. Vertical error bars represent the error propagated from the standard error of the mean of the rheological parameters.</figcaption></figure>

efficient model among a list of model candidates, while fitting the corresponding parameters? We consider a set of L candidate viscosity laws, $\eta_1^*, \eta_2^*, \dots, \eta_L^*$, each with an associated pressure scaling parameter $\kappa_i$ and their corresponding number of viscosity parameters $n_1, \dots, n_L$ (i.e, size of the set $\alpha$). Each viscosity model is associated to a physical loss,

$$ \mathcal{L}_{phy,i} = \| -\kappa_i \nabla P_{\theta}^* + \nabla \cdot (2\eta_i^* S_{\theta}^*) \|. \quad (17) $$

We build the total physical loss using a function analogous to the Akaike Information Criterion (AIC) used in Bayesian inference [39] or in some sparse regression algorithms using a $l_0$-norm regularization ([40, 41, 11]),

$$ \mathcal{L}_{phy,tot} = \sum_{i} \beta_i (\mathcal{L}_{phy,i} + \|\eta_i\|_0), \quad (18) $$
with $\beta_i \in [0, 1]$ an activation function following

$$ \beta_i(b_i) = \frac{\tanh rb_i) + 1}{2}, \qquad (19) $$

and $b_i$ are model activation parameters that need to be learned during the training process. The parameter $r$ sets the transition slope of the activation function; in the following we set $r = 10$. The total physical loss (18) depends on the number of physical parameters $\|\eta_i\|_0 = n_i$ involved in each model, avoiding a bias towards the selection of the model involving the largest number of fit parameters. We also need our algorithm to promote sparsity of the resulting rheological model by promoting the selection of a single rheological model during the learning process, i.e. setting one $\beta_i$ to 1 and the others to 0. We thus add a sparsity loss analogous to sparsity losses mentioned in [42, 43]

$$ \mathcal{L}_{\beta} = 2 \max_i(n_i) |1 - \sum_i \beta_i|. \qquad (20) $$

The factor $2 \max_i(n_i)$ is included to prevent the scenario where $\beta_i = 0$ for all $i$. Indeed, if $\beta_i = 0$ for all $i$ then $\mathcal{L}_{\beta} = 2 \max_i(n_i)$, which is greater than $n_i$ for all $i$. Therefore,

$$ (\mathcal{L}_{\beta} + \mathcal{L}_{phy,tot})_{\beta=\{0,...,0\}} > (\mathcal{L}_{\beta} + \mathcal{L}_{phy,tot})_{\beta=\{0,...,1,...,0\}}, \qquad (21) $$

which prevents the existence of a local minimum when all $\beta_i = 0$. The sparsity loss (20) also prevents the case where $\beta_i = 1/L$ for all models.

Indeed,

$$ (\mathcal{L}_{phy,tot} + \mathcal{L}_{\beta})_{\beta=\{\frac{1}{L},...,\frac{1}{L}\}} = \frac{1}{N} \left( \sum_{i=1}^{L} n_i + \mathcal{L}_{phy,i} \right) \qquad (22) $$

$$ > \min_i(n_i + \mathcal{L}_{phy,i}), \qquad (23) $$

i.e.

$$ (\mathcal{L}_{phy,tot} + \mathcal{L}_{\beta})_{\beta=\{\frac{1}{L},...,\frac{1}{L}\}} > (\mathcal{L}_{phy,tot} + \mathcal{L}_{\beta})_{\beta=\{0,...,1,...,0\}}. \qquad (24) $$

Taken together, the total model selection loss reads

$$ \mathcal{L}_{MS} = \mathcal{L}_{data} + w\mathcal{L}_{phy,tot} + \mathcal{L}_{p} + \mathcal{L}_{\beta}. \qquad (25) $$

We split the learning process into two phases. First, we perform $3 \times 10^4$ epochs using the model selection loss (25). Then, we stop the model selection by choosing the model with the highest value of $\beta_i$. The associated activation function is set to 1 and the others to 0. We finally perform additional epochs using the single-model loss (14) to fit the rheological parameters.

We illustrate the performance of this approach in figure 8. We provide the algorithm with three model candidates: Herschel-Bulkley, Carreau, and Papanastasiou. It is worth noting that the Papanastasiou model is a regularized version of the Herschel-Bulkley model for $n = 1$, meaning there exists a set of parameters $\alpha$ for which the Papanastasiou and Herschel-Bulkley viscosity laws are equivalent. To facilitate effective model selection, we constrain the Papanastasiou viscosity parameters $\alpha_0$ and $\alpha_1$ to be greater than 1 by adding
<figure><img src="image_11.png" /><figcaption>Figure 8: Rheology model selection. (a) Overview of the various viscosity functions $\eta^*(\dot{\gamma})$ possibly generated by the Herschel-Bulkley, Carreau and Papanastasiou models as their rheological parameters are varied. For each model, the green line indicates the actual case provided to the algorithm to test our model selection approach. With $\tilde{\gamma} = \frac{\dot{\gamma}-\min(\dot{\gamma})}{\max(\dot{\gamma})-\min(\dot{\gamma})}$ (b) Set of parameters with their associated viscosity provided to the algorithm to test our model selection approach. (c) Proportion of selected model over 30 tests when passing a simulation with the same set of parameters as in the column first column. (d) 3 examples of the convergence history of the 3 weights $\beta_i$</figcaption></figure>

a term $\max(0, 1 - \alpha_i)$ in the loss $L_p$. We select three examples in which the algorithm is applied to velocity data generated with one of the three models. In each example, we perform a total of 30 inference runs and show the resulting proportion of selected models in figure 8c. In all cases, the algorithm selects the correct rheological model, except for a few runs in which the Panastasiou and Herschel-Bulkley models are confused, especially when the true
model is the Herschel-Bulkley model. This can be attributed to the similarity between both
models for certain parameter sets (see figure 8). As shown in figure 8d, the dynamics of
the activation functions during training varies from one model to the other, but in all cases
the correct model is rapidly promoted compared to other model candidates. Overall, the
proposed model selection method is effective and easily extendable to a larger number of
candidate models.

6. Discussion

Since their introduction, physics-informed neural networks (PINNs) [18] have seen rapid
development and adoption across a range of physical and engineering disciplines. Their
strength particularly lies in solving inverse problems—inferring physical properties or gov-
erning models from data—by leveraging both observational data and physical laws to con-
strain solutions. This makes them especially well-suited for tasks such as rheological model
inference, where the underlying constitutive behavior is often complex or partially unknown.

Recent studies have begun to explore the application of PINNs to rheological prob-
lems, either by using classical rheometry data [26, 27] or, more recently, by learning from
velocity fields [25, 24]. In this work, we demonstrate that even relatively simple PINN ar-
chitectures—consisting of a single fully connected neural network—are capable of inferring
viscoplastic properties directly from noisy velocity data. Notably, our approach does not
require prior knowledge of boundary conditions, making it particularly promising for ana-
lyzing experimental data, where stress boundary conditions are typically unknown and the
no-slip assumption at solid boundaries may not hold [28, 29, 30, 31].

The methodology we propose is designed to be readily extensible to a wide range of vis-
coplastic fluid types. In addition, we introduce a practical strategy for model selection during
the inference process, enabling automatic identification of the most suitable constitutive law
from a set of candidate models. This adds a significant layer of flexibility and robustness to
the framework, potentially easing its adoption in experimental and industrial settings.

Several directions for future research naturally emerge from this work. A primary next
step is the application of our approach to experimental velocity fields. This would involve
benchmarking the method on systems with well-characterized rheology, followed by its de-
ployment on materials for which traditional characterization techniques are limited or infeasi-
ble. Furthermore, the framework can be extended to accommodate more complex rheological
behaviors, such as viscoelastic or viscoplastoelastic responses. One particularly compelling
objective is to generalize the model selection capability to span these broader rheological
classes, allowing the PINN to autonomously determine the most appropriate constitutive
description for a given dataset. Lastly, there remains a need for deeper theoretical un-
derstanding of PINN training dynamics, particularly in the context of inverse rheological
problems. Issues such as convergence rates, robustness to noise, and identifiability of model
parameters warrant further investigation. In the longer term, we envision that approaches
such as the one developed here will play a central role in advancing computational rheo-
logy and in addressing interdisciplinary challenges, notably in the study of biological and
geophysical materials, where conventional modeling approaches often fall short.
Acknowledgements

We thank Christophe Eloy, Hélène Delanoë-Ayari, Yoël Forterre, Philippe Roudot and Andonis Gerardos for useful discussion. The project leading to this publication has received funding from the “Investissements d’Avenir” French Government program managed by the French National Research Agency (ANR-16-CONV-0001) and from “Excellence Initiative of Aix-Marseille University – A*MIDEX”. We also acknowledge the French National Research Agency, under grant ANR-24-CE45-1959. The Centre de Calcul Intensif d’Aix-Marseille is acknowledged for granting access to its high performance computing resources.

Declaration of generative AI in the writing process

The authors used ChatGPT in order to improve the writing quality of their manuscript.
After using this tool/service, the authors reviewed and edited the content as needed and
take full responsibility for the content of the publication.

References

[1] P.-F. Lenne, V. Trivedi, Sculpting tissues by phase transitions, Nature Communications 13 (1) (2022) 664.

[2] A. Davaille, P. Carrez, P. Cordier, Fat plumes may reflect the complex rheology of the lower mantle, Geophysical Research Letters 45 (3) (2018) 1349–1354.

[3] A. Tommasi, D. Mainprice, G. Canova, Y. Chastel, Viscoplastic self-consistent and equilibrium-based modeling of olivine lattice preferred orientations: Implications for the upper mantle seismic anisotropy, Journal of Geophysical Research: Solid Earth 105 (B4) (2000) 7893–7908.

[4] P. Jop, Y. Forterre, O. Pouliquen, A constitutive law for dense granular flows, Nature 441 (7094) (2006) 727–730.

[5] É. Guazzelli, O. Pouliquen, Rheology of dense granular suspensions, Journal of Fluid Mechanics 852 (2018) P1.

[6] B. Etcheverry, Y. Forterre, B. Metzger, Capillary-stress controlled rheometer reveals the dual rheology of shear-thickening suspensions, Physical Review X 13 (1) (2023) 011024.

[7] S. L. Brunton, S. L. Brunton, J. L. Proctor, J. L. Proctor, J. N. Kutz, J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proceedings of the National Academy of Sciences of the United States of America (2016). doi:10.1073/pnas.1517384113.

[8] J. Bongard, H. Lipson, Automated reverse engineering of nonlinear dynamical systems, Proceedings of the National Academy of Sciences 104 (24) (2007) 9943–9948.

[9] S. H. Rudy, S. L. Brunton, J. L. Proctor, J. N. Kutz, Data-driven discovery of partial differential equations, Science advances 3 (4) (2017) e1602614.
[10] K. Kaheman, J. N. Kutz, S. L. Brunton, Sindy-pi: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics, Proceedings of the Royal Society A 476 (2242) (2020) 20200279.

[11] Z. Chen, Y. Liu, H. Sun, Physics-informed learning of governing equations from scarce data, Nature communications 12 (1) (2021) 6136.

[12] J. Colen, M. Han, R. Zhang, S. A. Redford, L. M. Lemma, L. Morgan, P. V. Ruijgrok, R. Adkins, Z. Bryant, Z. Dogic, et al., Machine learning active-nematic hydrodynamics, Proceedings of the National Academy of Sciences 118 (10) (2021) e2016708118.

[13] J. D. Christopher, N. T. Wimer, C. Lapointe, T. R. Hayden, I. Grooms, G. B. Rieker, P. E. Hamlington, Parameter estimation for complex thermal-fluid flows using approximate bayesian computation, Physical Review Fluids 3 (10) (2018) 104602.

[14] A. Kontogiannis, R. Hodgkinson, E. L. Manchester, Learning rheological parameters of non-newtonian fluids from velocimetry data, arXiv preprint arXiv:2408.02604 (2024).

[15] T. Fan, K. Xu, J. Pathak, E. Darve, Solving inverse problems in steady-state navier-stokes equations using deep neural networks, arXiv preprint arXiv:2008.13074 (2020).

[16] K. Xu, A. M. Tartakovsky, J. Burghardt, E. Darve, Learning viscoelasticity models from indirect data using deep neural networks, Computer Methods in Applied Mechanics and Engineering 387 (2021) 114124.

[17] S. L. Brunton, B. R. Noack, P. Koumoutsakos, Machine learning for fluid mechanics, Annual review of fluid mechanics 52 (1) (2020) 477-508.

[18] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational physics 378 (2019) 686-707.

[19] I. E. Lagaris, A. Likas, D. I. Fotiadis, Artificial neural networks for solving ordinary and partial differential equations, IEEE transactions on neural networks 9 (5) (1998) 987-1000.

[20] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, J. M. Siskind, Automatic differentiation in machine learning: a survey, Journal of machine learning research 18 (153) (2018) 1-43.

[21] S. Cai, Z. Mao, Z. Wang, M. Yin, G. E. Karniadakis, Physics-informed neural networks (pinns) for fluid mechanics: A review, Acta Mechanica Sinica 37 (12) (2021) 1727-1738.

[22] M. Raissi, A. Yazdani, G. E. Karniadakis, Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations, Science 367 (6481) (2020) 1026-1030.

[23] Q. Lou, X. Meng, G. E. Karniadakis, Physics-informed neural networks for solving forward and inverse flow problems via the boltzmann-bgk formulation, Journal of Computational Physics 447 (2021) 110676.
[24] R. Zhai, D. Yin, G. Pang, A deep learning framework for solving forward and inverse problems of power-law fluids, Physics of Fluids 35 (9) (2023).

[25] S. Thakur, M. Raissi, A. M. Ardekani, Viscoelasticnet: A physics informed neural network framework for stress discovery and model selection, Journal of Non-Newtonian Fluid Mechanics (2024) 105265.

[26] M. Mahmoudabadbozchelou, S. Jamali, Rheology-informed neural networks (rhinns) for forward and inverse metamodelling of complex fluids, Scientific reports 11 (1) (2021) 12015.

[27] M. Mahmoudabadbozchelou, K. M. Kamani, S. A. Rogers, S. Jamali, Digital rheometer twins: Learning the hidden rheology of complex fluids through rheology-informed graph neural networks, Proceedings of the National Academy of Sciences 119 (20) (2022) e2202234119.

[28] I. Cheddadi, P. Saramito, B. Dollet, C. Raufaste, F. Graner, Understanding and predicting viscous, elastic, plastic flows, The European Physical Journal E 34 (2011) 1–15.

[29] I. Golovkova, L. Montel, E. Wandersman, T. Bertrand, A. M. Prevost, L.-L. Pontani, Depletion attraction impairs the plasticity of emulsions flowing in a constriction, Soft Matter 16 (13) (2020) 3294–3302.

[30] S. Tlili, M. Durande, C. Gay, B. Ladoux, F. Graner, H. Delanoë-Ayari, Migrating epithelial monolayer flows like a maxwell viscoelastic liquid, Physical Review Letters 125 (8) (2020) 088102.

[31] S. L. Tlili, F. Graner, H. Delanoë-Ayari, A microfluidic platform to investigate the role of mechanical constraints on tissue reorganization 149 (20) dev200774. doi:10.1242/ dev.200774.

[32] T. C. Papanastasiou, Flows of materials with yield, Journal of rheology 31 (5) (1987) 385–404.

[33] D. Bonn, M. M. Denn, L. Berthier, T. Divoux, S. Manneville, Yield stress materials in soft condensed matter, Reviews of Modern Physics 89 (3) (2017) 035005.

[34] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).

[35] https://github.com/martinLARD/Pinn.git.

[36] S. Gsell, U. d'Ortona, J. Favier, Lattice-boltzmann simulation of creeping generalized newtonian flows: theory and guidelines, Journal of Computational Physics 429 (2021) 109943.

[37] S. Gsell, J. Favier, Direct-forcing immersed-boundary method: a simple correction preventing boundary slip error, Journal of Computational Physics 435 (2021) 110265.
[38] F. Chen, Y. Meng, K. Li, C. Yang, J. Yang, Global physics-informed neural networks (gpinns): from local point-wise constraint to global nodal association, arXiv preprint arXiv:2503.06403 (2025).

[39] K. P. Murphy, Machine learning: a probabilistic perspective, MIT press, 2012.

[40] J. Weston, A. Elisseeff, B. Schölkopf, M. Tipping, Use of the zero-norm with linear models and kernel methods, Journal of machine learning research 3 (Mar) (2003) 1439–1461.

[41] T. Pang, F. Nie, J. Han, X. Li, Efficient feature selection via $l_{2,0}$-norm constrained sparse regression, IEEE Transactions on Knowledge and Data Engineering 31 (5) (2018) 880–893.

[42] H. Narasimhan, Learning with complex loss functions and constraints, in: International Conference on Artificial Intelligence and Statistics, PMLR, 2018, pp. 1646–1654.

[43] A. Cotter, H. Jiang, M. Gupta, S. Wang, T. Narayan, S. You, K. Sridharan, Optimization with non-differentiable constraints with applications to fairness, recall, churn, and other goals, Journal of Machine Learning Research 20 (172) (2019) 1–59.