Simulation of the Mzymta River Flow Using the Regularized Physics-Informed Neural Networks

1. Introduction

Climate change is currently exerting a significant influence on water resources management. Specifically, alterations in the hydrological cycle complicate management practices, while the increasing frequency and intensity of extreme events place considerable strain on hydraulic structures and associated infrastructure [1].

Climate change is transforming water resources management from a routine task into a domain of continuous risk analysis and strategic planning. The development of flexible and resilient strategies has become a key priority for ensuring water security [2].

Mathematically, the study of flood phenomena can be formulated as the solution of forward and inverse hydrological problems by applying equations derived from the theory of shallow water. Physics-Informed Neural Networks (PINNs) have been increasingly applied in digital twin frameworks that combine physical models with data-driven methods [3,4,5]. In these models, the residuals of the partial differential equations are directly embedded into the loss function of the neural network. PINNs have been used for both forward solution and inverse parameter-identification problems [6,7,8]. The modern theory of inverse problems for hyperbolic equations is presented in a number of works [9,10]. Early studies in the field of numerical methods for solving inverse hydrological problems [11,12] applied various mathematical approaches to solve the Saint-Venant equations [10,13,14].

Early studies focused on assessing the capabilities of PINNs for solving hydrological problems, particularly for approximating one-dimensional shallow water equations. In [15,16], PINNs were used to assimilate observational data while solving the one-dimensional Saint-Venant equations. The performance of PINNs was evaluated by comparing the results with analytical solutions and numerical simulation data.

Numerical solutions of the one-dimensional equations were used to predict water level profiles in rivers, as well as to estimate the so-called roughness parameter [17]. PINNs models were implemented within the TensorFlow machine learning library using two neural networks. Different activation functions, numbers of layers, and numbers of neurons in the hidden layers were tested.

The aim of the dissertation work [18] was to investigate the applicability of PINNs for approximating the shallow water equations (SWE), which represent a system of PDEs used for modeling free-surface flows and the two-dimensional dam-break problem.

PINNs have also been applied to flow and transport problems in porous media, including the Darcy and Richards equations. For instance, the PINNs method was introduced for solving a coupled system (the Darcy equation and the advection–dispersion equation) and was tested over a range of Péclet numbers. The authors of [19] incorporated finite volume method-based numerical solutions as a reference to evaluate the PINNs results when determining the steady-state hydraulic head and concentration in the coupled system characterized by spatially varying conductivity and velocity fields. The PINNs approach enables the generation of numerical data in real time, which can then be used for modeling riverbed siltation. The prediction of key indicators plays a crucial role in reducing river pollution caused by siltation.

In [20], the researchers conducted a comparative analysis of the PINNs approach and traditional numerical methods for estimating flow velocity, pressure and density in water systems. The PINNs architecture was applied to solve the partial differential equation describing water infiltration into soil, enabling a numerical analysis of the infiltration process. The modeling results for vertical infiltration in the water-soil system demonstrated high accuracy and lower error levels compared to classical methods. Among the studied soils, water infiltrated fastest in light soil, followed by heavy soil, with the lowest infiltration rate observed in medium soil [21]. Several key soil parameters are required to reproduce water processes in soils using the Richardson–Richards equation (RRE). These include water retention curves (WRC) and hydraulic conductivity functions (HCF). Commonly used methods for obtaining WRC and HCF are often inadequate for such modeling.

In [22,23], it was shown that the PINNs approach can overcome this limitation and yield an inverse solution for the RRE, enabling the estimation of WRC and HCF from volumetric water content measurements without requiring initial or boundary conditions.

In [24], PINNs were employed to solve the inverse problem of unsaturated groundwater flow by incorporating observed volumetric water content data directly into the network loss function. Specifically, the Richards equation was coupled with the van Genuchten constitutive model, and PINNs were trained to recover both the spatially distributed solutions of the Richards equation and the unknown parameters of the van Genuchten model from sparse measurement data.

Berardi et al. [25] extended this approach to an adaptive inverse PINNs architecture for various transport models in porous media, including diffusion, advection–diffusion–reaction systems, and mobile–immobile phase formulations. In their method, transport parameters such as diffusion coefficients were treated as trainable variables, and the loss function components (data mismatch, initial and boundary conditions and PDE residuals) were adaptively weighted during training to ensure robust convergence.

In [26], the authors enhanced the PINNs methodology by proposing gradient-enhanced PINNs (gPINNs), in which the loss function includes not only PDE residual terms but also gradient information. Numerical experiments demonstrated that gPINNs achieve higher accuracy and faster convergence compared to standard PINNs, while requiring fewer collocation points for both forward and inverse problems.

Finally, Difonzo et al. [27] applied PINNs to estimate the horizon parameter in peridynamic models of linear microelasticity. The authors investigated several kernel functions (hat, distributed, and V-shaped kernels) and conducted numerical experiments in one- and two-dimensional settings. By treating the horizon size as an additional trainable parameter, they demonstrated that PINNs can accurately recover the nonlocal interaction scale directly from simulated data.

A specific application of PINNs was proposed in [28] for solving the shallow water equations (SWE) on a sphere. The trained PINNs outperformed conventional numerical methods such as finite difference, finite volume, and spectral methods in solving the differential equations.

Recent approaches in the application of PINNs and hybrid neural network architectures for studying hydrological processes in rivers are presented in [29,30,31,32,33,34]. In recent studies, PINNs and other machine learning methods have been actively applied to model climate-driven processes described by partial differential equations. Such approaches combine physical constraints with data-driven learning [3,4,5,35,36]. The development of PINNs architectures in fluid dynamics has focused on evolutionary deep neural networks [37,38,39,40]. A review of modern approaches in the application of PINNs to fluid dynamics is given in [41,42,43,44,45,46]. The work [47] presents an overview of the application of PINNs and PIKANs in physics-informed machine learning. The authors of [48] explore the use of an energy transformer approach to reconstruct three different flow fields from sparse measurement data.

Currently, there exist other approaches to constructing PINNs architectures based on operator neural networks (DeepONet) [49] and Fourier neural operators (FNO) [50]. However, their software implementation is more complex.

For solving inverse problems using PINNs, a careful selection of optimizers is required to compute the gradient of the loss function and find the optimal parameters of the neural network. A review of various optimizers, including GradNorm, PCGrad, Multi-Objective Optimization, and Self-Adaptive PINNs, was conducted by Italian authors in [25]. They developed a novel Adaptive-Inverse PINNs method, in which the loss function weights are computed at each iteration, and applied it to solve a one-dimensional inverse transport problem in porous media.

In the last five years, graph neural networks have been actively developed for studying hydrological processes in rivers. Notable research groups in this area include the team from Delft University of Technology, Delft, the Netherlands, with the mSWE-GNN project; the group from Southern University of Science and Technology, Shenzhen, China, working on graph models (ChebNet, GAT, GraphSAGE, GCNII, GCN, and GIN) as well as Entity-Aware LSTM; the team from Utah State University, Logan, UT, USA, developing a spatio-temporal graph neural network (STGNN); the collaboration between the US Army Corps of Engineers Vicksburg District, Mississippi, USA and the University of New Orleans, Louisiana, USA, working on the PINN-HEC-RAS project; the group from the University of Virginia, Virginia, USA, in collaboration with Nvidia, developing HydroGraphNet, which employs a GNN-KAN model; and the team from Rice University, Houston, TX, USA, with the FloodGNN-GRU project.

The implementation of digital twins for rivers enables a proactive approach to water resources management, allowing problems to be anticipated and mitigated in near real-time. Furthermore, digital twins can provide forecasts under scenario-based changes that cannot be reproduced in field conditions. River digital twins are often built upon hydrological models for runoff generation (e.g., SWAT, HEC-HMS), as well as hydrodynamic models (e.g., Delft3D, HEC-RAS, TELEMAC, MIKE 21). At the same time, artificial intelligence methods are developing rapidly, including PINNs, Graph Neural Networks (GNN), transformers for time series, DeepONet, and Gaussian Processes for improving predictions and uncertainty quantification. High-performance computing is widely employed, including GPUs, distributed computing, and cloud platforms. Modern hybrid digital twins combine physics-based laws with machine learning techniques.

The study by [15] employed PINNs to downscale a large-scale river model to a subgrid resolution. The proposed PINNs framework successfully solved the unsteady one-dimensional Saint-Venant equations, incorporated periodic tidal boundary conditions via Fourier feature embeddings, assimilated heterogeneous data sources, and was coupled with a large-scale river model to produce downscaled solutions.

The study [51] introduced FNOCL, a deep learning neural network for modeling unsteady natural river flows. FNOCL is based on the FNO by incorporating Convolutional Long Short-Term Memory (ConvLSTM) layers, combining robust performance with high efficiency in capturing spatiotemporal patterns. It is claimed that the network effectively handles the complexities of spatial and temporal variations. FNOCL accounts for complex river cross-sectional profiles, enabling accurate flow dynamics modeling in complex river morphologies. In this study, FNOCL is applied to two distinct river networks with varying morphological complexity. Furthermore, FNOCL demonstrates a significant speedup compared to HEC-RAS, achieving, in particular, a 50-fold increase in a complex river network.

The study [52] investigates the potential of using the FNO as a physics-informed surrogate for one-dimensional debris flow modeling, applied to the Rendinara–Morino system in central Italy. To address the challenges of limited observational data and the computational cost of high-fidelity solvers, a validated finite-volume shallow-water model is combined with operator learning to create a site-specific digital twin of debris flow dynamics. Overall, the study demonstrates that integrating physics-based synthetic data with operator learning architectures provides a viable path for constructing accurate, computationally efficient, and site-specific approximations for debris flow modeling under data-scarce conditions.

In this study, we focus on modeling the flow of the Mzymta River in the vicinity of Sochi, Russia. The research is based on the authors’ previous work [53], which was devoted to the inverse problem for the determination of the Manning coefficient using PINNs. The Mzymta is a typical mountain river, characterized by steep slopes, rapidly varying discharges, wave-dominated processes, and complex geometry. For such rivers, classical 1D models often struggle with unknown roughness, incomplete data, and sparse measurements. PINNs can employ physics as a regularizer and extract information even from a limited number of observations.

The paper is organized as follows. Section 1 contains an introduction that provides a brief overview of the research area. Section 2 covers the mathematical model for hydrological problems: Saint-Venant shallow water equations (Section 2.1), initial and boundary conditions (Section 2.2), non-dimensionalization of the SWE equations (Section 2.3) and the stationary problems statement peculiarities (Section 2.4). Section 3 covers the method for solving the problem: Fully-Connected Neural Network (Section 3.1), Physics-Informed Neural Networks (Section 3.2), smart boundary conditions (Section 3.3), hydraulic jump regularization (Section 3.4), Balancing the weights of the loss function components (Section 3.5). Section 4 covers the description of the considered Mzymta River and the data for its simulation. Section 5 covers the main results: full scale flow of the Mzymta River, stationary (Section 5.1) and non-stationary (Section 5.2); short flattening reach with transition from supercritical to subcritical flow, stationary (Section 5.3) and non-stationary (Section 5.4). Section 6 covers the discussion, where both the results obtained and the hardware implementation are discussed. Section 7 provides the conclusions and outlines future research plans.

2. Mathematical Model

To model the flow of the Mzymta River, we employ a one-dimensional unsteady Saint-Venant model. One-dimensional models are applicable when the flow exhibits a dominant direction. The Mzymta River flows within a single main channel, which makes the 1D approach physically justified. The Saint-Venant Shallow Water Equations (SWE) assume hydrostatic pressure distribution, which is valid when vertical accelerations are negligible. This condition is satisfied when the bed slope is less than 10% and the flow scale is much larger than the river depth. The slope of the Mzymta River does not exceed this threshold at the full-river scale. 1D models require fewer computational resources and less input data than 2D models, making them suitable for rapid scenario assessment. The 1D model is optimal when the flow is primarily confined to the channel and there is no extensive floodplain inundation that would require a 2D approach – which is exactly the case considered in this study.

2.1. Saint-Venant Shallow Water Equations

Consider a one-dimensional problem of water flow in a channel of constant width W [m] and length L [m] (Figure 1). The coordinate along the channel is denoted as $x\in [0,L]$, and the time variable as t. The channel has a bottom elevation profile z [m] above sea level. At each point, the flow has a depth h [m] and a longitudinal cross-sectionally averaged velocity u [m/s]. The water discharge through a cross-section of the river is q [${\mathrm{m}}^2$/s] $=hu$. The flow is accelerated by gravity g [m/${\mathrm{s}}^2$] and decelerated by bottom friction, characterized by the bottom roughness coefficient n [s/${\mathrm{m}}^{1/3}$] according to the Manning formula [54].

The flow in the shallow channel is described by the Saint-Venant Shallow Water Equations. The continuity equation takes the following form:

$${\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial \left(hu\right)}{\partial x}}=0,$$

(1)

and the momentum equation is written as follows

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}+g{\displaystyle \frac{\partial h}{\partial x}}+g{\displaystyle \frac{\partial z}{\partial x}}+g{n}^2{\displaystyle \frac{u^2}{h^{4/3}}}=0.\hfill \end{array}$$

(2)

We assume that the bed slope ${\partial }_{x}z$ is negative and the flow velocity u is positive, which reflects the physical condition that the flow is directed downstream and cannot go uphill.

2.2. Initial and Boundary Conditions

The system (1)-(2) is of mixed hyperbolic-parabolic type. To close it, initial and boundary conditions must be specified. We seek the solution in the domain $\mathsf{\Omega }:=[0,T]\times [0,L]$. The initial conditions are written as:

$$q(0,x)=q_0\left(x\right),h(0,x)=h_0\left(x\right).$$

Since the system consists of two first-order equations, two boundary conditions must be specified: one for the discharge and one for the depth. The location of each boundary condition (inlet or outlet) depends on the specific conditions of the characteristics at the boundaries of the domain. Since water always flows into the domain at the inlet boundary, the boundary condition for discharge is imposed at the inlet cross-section:

$$q(t,0)=q_{\mathrm{in}}\left(t\right).$$

For the depth, the imposition of the boundary condition depends on the flow regime. In subcritical flow, the characteristics propagate upstream, and the depth is specified at the outlet. In supercritical flow, conversely, the characteristics propagate downstream, and the depth is specified at the inlet:

$$h(t,L)=h_{\mathrm{out}}\left(t\right)\mathrm{or}h(t,0)=h_{\mathrm{in}}\left(t\right).$$

However, in the problem under consideration, mixed flow regimes may occur (subcritical in one region and supercritical in another). In this case, the boundary conditions become more complicated. For example, a situation may arise where the depth h is prescribed both at the inlet (supercritical regime) and at the outlet (subcritical regime). Alternatively, the depth may not be prescribed at either the inlet or the outlet boundaries (if the flow is subcritical at the inlet and supercritical at the outlet).

2.3. Non-Dimensionalization of the SWE

To solve the problem using PINNs, it is necessary to non-dimensionalize the equations and the boundary conditions. This is because neural networks can effectively capture the nonlinearity of the problem when the input data lie in the vicinity of unity, since it is precisely in this range that the activation function $\sigma$ exhibits substantially nonlinear behavior.

For non-dimensionalization of the equations let us introduce the characteristic scale parameters: time T, length L, depth H, and velocity U. We write down the relations between the dimensional and dimensionless (denoted by the symbol $\tilde{\phantom{x}}$) variables as follows:

$$\begin{array}{cc}\hfill & t=T\tilde{t},x=L\tilde{x},h=H\tilde{h},u=U\tilde{u},q=HU\tilde{q},\hfill \\ \hfill & s={\partial }_{x}z=-\mathrm{Sl}\tilde{s},\hfill \\ \hfill & s_{\mathrm{f}}=n^2{\displaystyle \frac{u^2}{h^{4/3}}}=\left({\displaystyle \frac{n^2{U}^2}{H^{4/3}}}\right)\left({\displaystyle \frac{{\tilde{u}}^2}{{\tilde{h}}^{4/3}}}\right)=\mathrm{Sl}{\tilde{s}}_{\mathrm{f}}\hfill \end{array}$$

We decompose the bottom slope ${\partial }_{x}z$ into the product of a characteristic slope Sl and a normalized slope $\tilde{s}$. Note that ${\partial }_{x}z$ and ${\partial }_x\tilde{z}$ take negative values, while Sl and $\tilde{s}$ are positive.

We introduce the dimensionless complexes characterizing the system:

• The Strouhal number $\mathrm{St}={\displaystyle \frac{UT}{L}}$ which characterizes the time scale of the problem and its unsteadiness;
• The Froude number $\mathrm{Fr}={\displaystyle \frac{U}{\sqrt{gH}}}$, which characterizes the gravitational component of acceleration;
• The Manning number $\mathrm{Sl}={\displaystyle \frac{n^2{U}^2}{H^{4/3}}}$, which characterizes the gravitational acceleration from the slope balanced by the Manning friction force;
• The aspect ratio $L/H$.

Note that the parameters $\mathrm{Sl},U,H$ are chosen in such a way that the slope and friction are characterized by a single characteristic parameter Sl, corresponding to uniform steady Manning flow.

Then the continuity equation is rewritten as

$${\displaystyle \frac{1}{\mathrm{St}}}{\displaystyle \frac{\partial \tilde{h}}{\partial \tilde{t}}}+{\displaystyle \frac{\partial \left(\tilde{h}\tilde{u}\right)}{\partial \tilde{x}}}=0.$$

The momentum equation is rewritten as

$${\displaystyle \frac{1}{\mathrm{St}}}{\displaystyle \frac{\partial \tilde{u}}{\partial \tilde{t}}}+\tilde{u}{\displaystyle \frac{\partial \tilde{u}}{\partial \tilde{x}}}+{\displaystyle \frac{1}{{\mathrm{Fr}}^2}}\left({\displaystyle \frac{\partial \tilde{h}}{\partial \tilde{x}}}+{\displaystyle \frac{L}{H}}\mathrm{Sl}({\tilde{s}}_{\mathrm{f}}-\tilde{s})\right)=0.$$

The scaled initial conditions are rewritten as follows:

$$\tilde{q}(0,\tilde{x})={\tilde{q}}_0\left(\tilde{x}\right),\tilde{h}(0,\tilde{x})={\tilde{h}}_0\left(\tilde{x}\right).$$

The scaled boundary conditions are rewritten as follows:

$$\begin{array}{cc}\hfill & \tilde{q}(\tilde{t},0)={\tilde{q}}_{\mathrm{in}}\left(\tilde{t}\right);\hfill \\ \hfill & \tilde{h}(\tilde{t},1)={\tilde{h}}_{\mathrm{out}}\left(\tilde{t}\right),\hfill & \hfill \mathrm{if}\mathrm{fr}(\tilde{t},1)<1;\\ \hfill & \tilde{h}(\tilde{t},0)={\tilde{h}}_{\mathrm{in}}\left(\tilde{t}\right),\hfill & \hfill \mathrm{if}\mathrm{fr}(\tilde{t},0)>1.\end{array}$$

(3)

Here $\mathrm{fr}=u/\sqrt{gh}=\mathrm{Fr}\tilde{u}/\sqrt{\tilde{h}}$ is the local Froude number, characterizing the flow at the precise point in the vicinity of the inlet or outlet boundary.

2.4. Stationary Problem

The steady-state solution of the SWE (continuity (1) and momentum (2)) reduces to the backwater equation [55,56]:

$${\displaystyle \frac{\partial h}{\partial x}}={\displaystyle \frac{s-s_{\mathrm{f}}}{1-{\mathrm{fr}}^2}}\mathrm{or}{\displaystyle \frac{\partial \tilde{h}}{\partial \tilde{x}}}{\displaystyle \frac{H}{L}}=\mathrm{Sl}{\displaystyle \frac{\tilde{s}-{\tilde{s}}_{\mathrm{f}}}{1-{\mathrm{Fr}}^2{\tilde{u}}^2/\tilde{h}}}$$

(4)

From equation (4), it follows that the sign of the derivative ${\partial }_{x}h$ is determined by the sign of the balance between bed slope and friction $(s-s_{\mathrm{f}})$ and the sign of the denominator $(1-{\mathrm{fr}}^2)$. Based on the combinations of these signs, various steady-state flow regimes are identified. The regime where the water depth increases in the downstream direction (${\partial }_{x}h>0$) is called a backwater curve (or rising profile), whereas the regime where the depth decreases (${\partial }_{x}h<0$) is called a drawdown curve (or falling profile).

A change in the sign of either of these terms results in a change of the flow regime. The transition between the regimes is defined by the normal depth and critical depth.

The normal depth $h_n$ is defined as the depth at which the bed slope equals the friction slope ($s=s_{\mathrm{f}}$). Under this condition, the gravitational force is in equilibrium with the boundary friction, and the water surface is parallel to the bed. The normal depth $h_n$ is obtained from the Manning equation:

$$h_n={\left({\displaystyle \frac{nq}{\sqrt{s}}}\right)}^{3/5}.$$

(5)

The critical depth $h_c$ is defined as the depth at which the Froude number equals unity ($\mathrm{Fr}=1$). It represents the transition between subcritical ($\mathrm{Fr}<1$) and supercritical ($\mathrm{Fr}>1$) flow. The critical depth $h_c$ is given by:

$$h_c={\left({\displaystyle \frac{q^2}{g}}\right)}^{1/3}.$$

(6)

The relative magnitudes of the actual depth h, the normal depth $h_n$, and the critical depth $h_c$ determine the type of water surface profile. If $h_n>h_c$, the bed slope is classified as mild (M); if $h_n<h_c$, it is steep (S); if $h_n=h_c$, it is critical [57,58].

A change in the flow regime can occur either smoothly or abruptly. A transition through the normal depth is smooth and is characterized by a water surface parallel to the bed (${\partial }_{x}h=0$). A transition through the critical depth is typically abrupt (${\partial }_{x}h\to \infty$), manifesting as either a hydraulic jump (see Figure 1) or a waterfall. An exception occurs when the normal depth coincides with the critical depth. In this case, both the numerator and the denominator in Eq. (4) tend to zero and the transition may be smooth.

3. Method

To solve the Shallow Water Equations, the Physics-Informed Neural Networks method is employed in this work. The core idea of the approach is to approximate the solution of the governing equations using a Fully Connected Neural Network (FCNN).

3.1. Fully Connected Neural Network

A fully connected neural network represents an approximation of some mapping $G:\mathbf{x}\to \mathbf{u}$ of the input data vector $x\in {\mathbb{R}}^n$ (e.g., $\mathbf{x}=\{x,t\}$ in our problem) to the output data vector $\mathbf{u}\in {\mathbb{R}}^m$ (e.g., $\mathbf{u}=\{u,h\}$) as a sequence of several nonlinear operators (so-called layers).

Consider, as an example, a neural network with two hidden layers $\mathit{l}=\left\{l_i\right\}$ and $\mathbf{m}=\left\{m_j\right\}$.

$$\begin{array}{cc}\hfill & \mathbf{u}\left(\mathbf{x}\right)={\mathrm{FNN}}_{\mathbf{u}}\left(\mathbf{x}\right)=\mathbf{u}\left(\mathbf{m}\left(\mathit{l}\left(\mathbf{x}\right)\right)\right),\hfill \\ \hfill & \mathit{l}\left(\mathbf{x}\right)=\sigma \left({\mathbf{W}}^{\mathbf{x}\mathit{l}}·\mathbf{x}+{\mathbf{b}}^{\mathit{l}}\right)\hfill \\ \hfill & \mathbf{m}\left(\mathit{l}\right)=\sigma \left({\mathbf{W}}^{\mathit{l}\mathbf{m}}·\mathit{l}+{\mathbf{b}}^{\mathbf{m}}\right)\hfill \\ \hfill & \mathbf{u}\left(\mathbf{m}\right)=\sigma \left({\mathbf{W}}^{\mathbf{m}\mathbf{u}}·\mathbf{m}+{\mathbf{b}}^{\mathbf{u}}\right)\hfill \end{array}$$

(7)

Here ${\mathbf{W}}^{\mathit{l}}$, ${\mathbf{W}}^{\mathbf{m}}$, ${\mathbf{W}}^{\mathbf{u}}$ are linear transition operators from layer to layer, mapping vectors of different dimensions; ${\mathit{l}}_0$, ${\mathbf{m}}_0$, ${\mathbf{u}}_0$ are bias vectors; and $\sigma$ is a nonlinear activation function. The first universal approximation theorem [59] required the activation function to be bounded, continuous, and sigmoidal (i.e., monotone and non-constant). However, subsequent theoretical developments have demonstrated that these conditions can be substantially relaxed. In particular, it has been shown that the activation function must be non-polynomial [60], piecewise continuous, and may be unbounded [61] – a property that encompasses modern activation functions such as ReLU, SeLU, and SiLU/Swish.

A neural network can be viewed as a computational graph. In this graph the vertices store numerical values (such as inputs x, intermediate activations l and m, or outputs u), and the edges represent simple arithmetic operations (e.g., weighted summation W, addition of bias b, and application of a nonlinear activation function $\sigma$).

To ensure that the neural network provides a good approximation of the target mapping, its coefficients (weights) must be trained using gradient-based optimization of the loss function, which measures the discrepancy between the network output $\mathbf{u}\left(\mathbf{x}\right)$ and the training data $\tilde{\mathbf{u}}\left(\mathbf{x}\right)$.

$$\mathcal{L}={\mathcal{L}}_{\mathrm{DATA}}=\parallel \mathbf{u}\left(\mathbf{x}\right)-\tilde{\mathbf{u}}\left(\mathbf{x}\right)\parallel .$$

(8)

The gradient-based optimization process is described by a gradient flow equation

$${\displaystyle \frac{\partial \mathbf{W}}{\partial \tau }}=-{\displaystyle \frac{\partial \mathcal{L}\left(\mathbf{W}\right(\tau \left)\right)}{\partial \mathbf{W}}}.$$

(9)

where $\tau$ is an iteration variable, $\mathbf{W}\left(\tau \right)$ are the network parameters.

Advanced gradient-based methods can introduce additional complexity to the governing equation (9). In particular, the Adam optimizer [62], employed in the present study, leads to a coupled system of equations through its adaptive moment estimation mechanism [63]:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \mathbf{W}}{\partial \tau }}& =-{\displaystyle \frac{\mathbf{m}\left(\tau \right)}{\sqrt{\mathbf{v}\left(\tau \right)}+\epsilon }},\hfill \\ \hfill {\displaystyle \frac{\partial \mathbf{m}}{\partial \tau }}& =d\left(\tau \right)\left[{\displaystyle \frac{\partial \mathcal{L}\left(\mathbf{W}\right)}{\partial \mathbf{W}}}\right]-r\left(\tau \right)\mathbf{m}\left(\tau \right),\hfill \\ \hfill {\displaystyle \frac{\partial \mathbf{v}}{\partial \tau }}& =p\left(\tau \right){\left[{\displaystyle \frac{\partial \mathcal{L}\left(\mathbf{W}\right)}{\partial \mathbf{W}}}\right]}^{\circ 2}-q\left(\tau \right)\mathbf{v}\left(\tau \right),\hfill \end{array}$$

(10)

where $\mathbf{m}\left(\tau \right)$ and $\mathbf{v}\left(\tau \right)$ are the first and second moment estimates, respectively, $\epsilon$ is a small constant for numerical stability, and $d\left(\tau \right)$, $r\left(\tau \right)$, $p\left(\tau \right)$, $q\left(\tau \right)$ are time-dependent coefficients that depend on the hyperparameters of the algorithm. The notation ${[·]}^{\circ 2}$ denotes element-wise squaring, which produces a vector of the same dimension as the input. This is in contrast to the norm operator, which aggregates the components into a single scalar value.

Optimization methods and their corresponding continuous-time formulations may involve not only first-order derivatives but also second-order derivatives. For example, Nesterov accelerated gradient method [64] includes a second-order derivative with respect to the iteration variable [65]. Newton methods use the second-order derivatives (the Hessian matrix), but this leads to a significant increase in computational complexity. Quasi-Newton methods, such as BFGS and its limited-memory variant L-BFGS [66,67,68,69,70], do not explicitly compute second-order derivatives. Instead, they construct an approximation of the Hessian matrix (the matrix of second-order derivatives) or its inverse using only gradient information from successive iterations.

To calculate the derivatives of the loss function with respect to network parameters automatic differentiation [71,72,73,74] is used. This procedure consists of two phases: a forward pass and backpropagation. The forward pass consists of the sequential evaluation of the vertices of the graph, propagating information from the input layer through the hidden layers to the output layer. Conversely, the backward pass (backpropagation) computes the gradients of the loss function with respect to both the input variables x and the trainable parameters (weights W and biases b) by traversing the graph in reverse order and applying the chain rule:

$${\displaystyle \frac{\partial \mathcal{L}\left(\mathbf{u}\right(\mathbf{m}\left(\mathit{l}\right)\left)\right)}{\partial W}}=\sum _{i,k,s}{\displaystyle \frac{\partial \mathcal{L}}{\partial u_i}}{\displaystyle \frac{\partial u_i}{\partial m_k}}{\displaystyle \frac{\partial m_k}{\partial l_s}}{\displaystyle \frac{\partial l_s}{\partial W}}.$$

(11)

In summary, we present the algorithm for minimizing the loss function. The training procedure consists of the following steps:

1.

Initial guess for weight matrices W, biases b, batch of input vectors ${\mathbf{x}}^{\left(i\right)}$;

2.

Forward Pass: calculate layers l and m, output u using formula (7);

3.

Calculate Loss Function using formula (8);

4.

Back Propagation: calculate derivatives of Loss function with respect to network weights using chain rule (11).

5.

Recalculate weights using Adam method (10).

3.2. Physics-Informed Neural Networks

The neural network method allows one to approximate an unknown mapping using a sufficient amount of data ${\mathbf{u}}_{\mathrm{DATA}}$. These data can originate from field observations, laboratory experiments, or numerical simulations performed using conventional models (e.g., finite difference or finite volume methods).

However, neural networks can leverage not only pure data for training but also indirect information, such as the differential equations governing the underlying process. This is the core idea of the PINNs [3,6]: the governing partial differential equations (PDE) are embedded into the loss function as residual terms, allowing the network to learn solutions that satisfy the physical laws even in regions where no direct data are available. This approach effectively regularizes the training process and enables accurate predictions with sparse or noisy observations.

The loss function $ℒ$ (8) is enriched by additional terms describing the physics of the process:

$$\mathcal{L}={\displaystyle \sum _{i=1}^I}{\lambda }_i{\mathcal{L}}_i={\lambda }_{\mathrm{DATA}}{\mathcal{L}}_{\mathrm{DATA}}+{\lambda }_{\mathrm{PDE}}{\mathcal{L}}_{\mathrm{PDE}}+{\lambda }_{\mathrm{IC}}{\mathcal{L}}_{\mathrm{IC}}+{\lambda }_{\mathrm{BC}}{\mathcal{L}}_{\mathrm{BC}}+{\lambda }_{\mathrm{REG}}{\mathcal{L}}_{\mathrm{REG}}.$$

These terms include the residual of the PDE

$${\mathcal{L}}_{\mathrm{PDE}}=\parallel \mathrm{PDE}(u,h,u_t,u_x,h_t,h_x)\parallel ,$$

initial conditions

$${\mathcal{L}}_{\mathrm{IC}}=\parallel \mathbf{u}(x,0)-{\mathbf{u}}_0\left(x\right)\parallel .$$

boundary conditions

$${\mathcal{L}}_{\mathrm{BC}}={\lambda }_{u_{\mathrm{in}}}\parallel u(0,t)-u_{\mathrm{in}}\left(t\right)\parallel +{\lambda }_{h_{\mathrm{in}}}\parallel h(0,t)-h_{\mathrm{in}}\left(t\right)\parallel +{\lambda }_{h_{\mathrm{out}}}\parallel h(0,t)-h_{\mathrm{out}}\left(t\right)\parallel .$$

and some regularization terms

$${\mathcal{L}}_{\mathrm{REG}}=\parallel \mathrm{REG}(u,h,u_t,u_x,h_t,h_x)\parallel .$$

To compute $\mathrm{PDE}(u,h,u_t,u_x,h_t,h_x)$ and $\mathrm{REG}(u,h,u_t,u_x,h_t,h_x)$ it is necessary to compute the derivatives ${\partial }_x\tilde{u}$, ${\partial }_t\tilde{u}$, ${\partial }_x\tilde{h}$, ${\partial }_t\tilde{h}$ using the automatic differentiation and the chain rule (11).

3.3. Smart Boundary Conditions

To address the problem of imposing boundary conditions (3) depending on the Froude number in our model, we employed smart boundary conditions with adaptive weights for the boundary condition components:

$${\lambda }_{h_{\mathrm{in}}}={\displaystyle \frac{1+tanh\left(k\right(\mathrm{fr}-1\left)\right)}{2}},{\lambda }_{h_{\mathrm{out}}}={\displaystyle \frac{1-tanh\left(k\right(\mathrm{fr}-1\left)\right)}{2}}.$$

When $\mathrm{fr}\gg 1$ (supercritical flow), $tanh\left(k\right(\mathrm{fr}-1\left)\right)\to 1$, so that ${\lambda }_{h_{\mathrm{in}}}\to 1$ and ${\lambda }_{h_{\mathrm{out}}}\to 0$. Thus, the inlet depth boundary condition is enforced, which is consistent with the characteristic propagation direction for supercritical flow. Conversely, when $\mathrm{fr}\ll 1$ (subcritical flow), $tanh\left(k\right(\mathrm{fr}-1\left)\right)\to -1$, yielding ${\lambda }_{h_{\mathrm{in}}}\to 0$ and ${\lambda }_{h_{\mathrm{out}}}\to 1$, enforcing the outlet depth boundary condition. The parameter $k>0$ controls the sharpness of the transition between the two regimes.

3.4. Hydraulic Jump Regularization

To adaptively control the regularization in the vicinity of the hydraulic jump (shock), we introduce several variants of a shock detector $S\left(x\right)$.

The first idea is to identify the shock using the local gradient of the solution:

$$S_{h_x}(x,t)=1-{\displaystyle \frac{{\left({\partial }_{x}h\right)}^2}{\underset{x,t}{max}{\left({\partial }_{x}h\right)}^2}}.$$

This quantity measures the relative magnitude of the local gradients with respect to the global maximum gradients. In smooth regions, $S\left(x\right)$ is close to 1, while near shocks or steep gradients, $S_{h_x}\left(x\right)$ approaches 0.

The second shock identification criterion relies on the local Froude number fr, which serves as an indicator of the region where a sudden flow regime transition is likely to take place:

$$S_{\mathrm{fr}}(x,t)=ln\left({\mathrm{fr}}^2(x,t)\right).$$

The logarithmic formulation ensures consistent scaling across supercritical and subcritical flow regimes.

The adaptive switching coefficient can be computed using an exponential function

$$\lambda (x,t)=1-exp\left({\displaystyle \frac{-S(x,t)}{\sigma }}\right),$$

or a hyperbolic tangent function:

$$\lambda (x,t)={\displaystyle \frac{1}{2}}\left(1+tanh\left({\displaystyle \frac{S(x,t)-S_0}{\sigma }}\right)\right),$$

where $S\left(x\right)$ is the shock detector, the coefficient $\sigma$ controls the sharpness of the function and $S_0$ is the threshold of function transition. To avoid completely disabling the equation, we use a weighted ${\lambda }_{\mathrm{REG}}$:

$${\lambda }_{\mathrm{REG}}(x,t)={\lambda }_{min}+(1-{\lambda }_{min})\lambda (x,t),$$

where ${\lambda }_{min}$ is a small positive constant that ensures the equation is never fully turned off. When $S\left(x\right)\approx 0$ (in the vicinity of a shock), ${\lambda }_{\mathrm{REG}}(x,t)\approx {\lambda }_{min}$; when $S(x,t)\approx 1$ (in smooth regions), ${\lambda }_{\mathrm{REG}}(x,t)\approx 1$. This allows us to apply regularization selectively, e.g., by diminishing the contribution of the momentum equation (2) to the loss while avoiding over-smoothing in the rest of the domain:

$${\mathcal{L}}_{\mathrm{PDE}\mathrm{REG}}(x,t)={\lambda }_{\mathrm{REG}}(x,t){\mathcal{L}}_{\mathrm{PDE}}(x,t).$$

To sharpen the discontinuity, we introduce a penalty term that minimizes the area where regularization is applied, thereby confining the transition to a narrow region:

$${\mathcal{L}}_{\mathrm{SHOCK}}(x,t)=(1-{\lambda }_{\mathrm{REG}}(x,t)).$$

In addition, a total variation (TV) regularization term is applied within the regularized region:

$${\mathcal{L}}_{\mathrm{TV}}=(1-{\lambda }_{\mathrm{REG}}(x,t)){\parallel {\partial }_{x}h\parallel }_{L_1}.$$

This approach is motivated by the Rudin–Osher–Fatemi (ROF) algorithm for noise removal [75]. The ROF model is well-known for its edge-preserving properties, as it penalizes the $L_1$ norm of the gradient rather than the $L_2$ norm, thereby allowing for sharp discontinuities in the solution.

3.5. Balancing the Weights of the Loss Function Components

In PINNs, one of the critical challenges is balancing the weights ${\lambda }_i$ of the loss function components ${\mathcal{L}}_i$.

Several approaches exist for determining the basis on which the weights should be balanced. In the first (and simplest) approach, it is postulated that the components of the loss function should be balanced:

$${\lambda }_i\left(\tau \right){\mathcal{L}}_i\left(\tau \right)\to {\mathcal{L}}_{\mathrm{TARGET}}.$$

(12)

However, during training, it is not the magnitudes of the loss components themselves that matter most, but rather their rates of convergence. Therefore, in the second approach, it is proposed to balance the learning rates of the individual components:

$${\lambda }_i\left(\tau \right){\displaystyle \frac{\partial {\mathcal{L}}_i\left(\tau \right)}{\partial \tau }}\to {\displaystyle \frac{\partial {\mathcal{L}}_{\mathrm{TARGET}}}{\partial \tau }}.$$

(13)

This method is more computationally expensive, as it requires either tracking the history of the loss evolution or explicitly computing the rate of change of the loss with respect to the network weights.

In our work, we employ the method described in Eq. (12) with a time-dependent target ${\mathcal{L}}_{\mathrm{TARGET},i}\left(\tau \right)$. By introducing a dependence on the iteration number $\tau$, the method can alternately focus on the boundary conditions, the continuity equation, or the momentum equation.

In our calculations, we adopt the following strategy. The weights of the boundary conditions are fixed at ${\lambda }_{\mathrm{BC}}=1$, while the weights of the governing equations are adaptively adjusted toward their target values

$${\lambda }_{{\mathrm{PDE}}_i}(\tau +\mathsf{\Delta }\tau )={\lambda }_{{\mathrm{PDE}}_i}\left(\tau \right)·CLIP\left({\displaystyle \frac{{\mathcal{L}}_{{\mathrm{TARGET}}_i}\left(\tau \right)}{{\mathcal{L}}_{{\mathrm{PDE}}_i}\left(\tau \right)}},1/\delta ,\delta \right).$$

where $CLIP(·)$ is the procedure of value clipping, $\delta$ is the limiter of maximal change of ${\lambda }_{{\mathrm{PDE}}_i}$.

The target value ${\mathcal{L}}_{{\mathrm{TARGET}}_i}\left(\tau \right)$ is set by the formula

$${\mathcal{L}}_{{\mathrm{TARGET}}_i}\left(\tau \right)=\underset{j}{max}(k{\mathcal{L}}_{{\mathrm{BC}}_j},{\mathcal{L}}_{{\mathrm{PDE}}_i},{\mathcal{L}}_{ϵ})·{10}^{a+(b-a){\displaystyle \frac{1+sin\left({\omega }_i\tau \right)}{2}}}.$$

(14)

where ${\mathcal{L}}_{ϵ}$ is the lower limit of amplitude base coefficient, k is the fraction coefficient of BC loss, a and b are the lower and upper limits of the amplitude oscillation, ${\omega }_i$ is the frequency of amplitude oscillation.

Figure 2 presents the typical evolution of the weighted loss functions ${\mathcal{L}}_i$ (left) and the corresponding loss weights ${\lambda }_i$ (right). The left panel shows that the boundary conditions converge rapidly to a limit of ${10}^{-7}$. This limit is determined by the learning rate. Therefore monotonic convergence requires reducing the learning rate. In contrast to the boundary conditions, the PDE residuals converge at a slower rate and display low-frequency oscillations within the target interval specified by Eq. (14). The right panel shows the history of the loss weights ${\lambda }_i$. The weights for the BC are constant and equal to unity throughout the entire training process. The weights for the PDE start from a very low level (${10}^{-8}$) and are then adaptively recomputed to satisfy the target value defined in Eq. (13). Note that the average ${\lambda }_{\mathrm{PDE}}$ increases over iterations (right panel), while their weighted product oscillates near a constant level (left panel), implying that the unweighted PDE losses ${\mathcal{L}}_{\mathrm{PDE}}$ decrease accordingly. It is worth noting that this strategy enables the exploration of a wide range of loss weight combinations, thereby promoting a flexible search for the global minimum and reducing the risk of converging to local minima.

4. The Mzymta River: Site Description and Hydrological Characteristics

This work is devoted to the numerical modeling of the flow of the Mzymta River, which flows in the Sochi region, Russia.

The Mzymta River is the longest and most abundant river on the Black Sea coast of Russia. It is located within the Adler District of Sochi city and flows through popular resort areas, including Krasnaya Polyana, Esto-Sadok and Rosa Khutor. The river is approximately 89 km long, and its drainage basin covers an area of 885 ${\mathrm{km}}^2$.

• Flow regime: This is a typical mountain river with a sharp drop in altitude and a fast current. It originates from the high-mountain Lake Kardyvach in the upper reaches of the Caucasus Mountains and flows into the Black Sea near Adler. In its lower reaches, the river becomes calmer and approaches a lowland character. The river has a mixed feeding regime, with spring–summer floods and rainfall-induced flash floods. The mean annual discharge is about 50 ${\mathrm{m}}^3$/s (with a maximum recorded discharge of 764 ${\mathrm{m}}^3$/s).
• Tributaries: The main tributaries include the Kepsha, Chvizhepse, Beshenka, Laura, and Pslukh.
• Nature and landscape: The Mzymta valley is characterized by picturesque mountain scenery. In its upper reaches, there are beautiful lakes, and the river itself flows through gorges and valleys.
• Utilization: The river holds significant touristic value and is an important recreational area in the Krasnaya Polyana and Rosa Khutor region. It is popular among outdoor enthusiasts, particularly for rafting and hiking. A road and a railway run along the river. Its waters are used for water supply and hydroelectric power generation. Near the village of Kazachiy Brod lies the famous Adler trout farm. In the lower reaches, there are major settlements such as Adler and Sirius, featuring sports centers, stadiums, a university, hotels, and shopping centers.

For simulation of the Mzymta River we need to determine its basic model parameters. All parameters are presented in the Table 1.

The Manning roughness coefficient can be estimated using standard hydraulic tables, such as those presented in [58]. Among the channel types listed therein, the most appropriate category for the Mzymta River is “Mountain streams, no vegetation in the channel, banks usually steep, with trees and brush along banks submerged at high stages.” The corresponding range of the Manning roughness coefficient is $0.03$–$0.07$, depending on the degree of boulder content and vegetation cover.

For simulation of the Mzymta River the following data were collected:

• The longitudinal bottom profile $z\left(x\right)$
• Time series of water discharge and water level data for the river over the period from 01.06.2023 to 31.08.2023.

4.1. The Longitudinal Bottom Profile

The longitudinal bottom profile $z\left(x\right)$ for the Mzymta River was constructed using digital topographic map data from Google Earth. Figure 3 presents the river basin map.

Note that the elevation data in Google Earth are based on SRTM (Shuttle Radar Topography Mission) data, whose nominal accuracy is about 16 meters [76]; however, in narrow, deep gorges, the error can reach up to 100 meters [77]. These errors in bed elevation estimates are related to the specific characteristics of the radar surveying technique. In steep terrain, the radar signal interacts primarily with the valley walls rather than the valley floor, leading to a systematic overestimation of the bottom elevation. This effect is particularly pronounced in narrow canyons, where the line-of-sight of the radar beam is partially obstructed, resulting in significant discrepancies in the derived longitudinal bed profile. As a consequence, the slopes obtained from SRTM data in such regions should be interpreted with caution, and the model results may require calibration against in-situ measurements for accurate predictions.

Figure 4 (left) shows the obtained bed elevation profile above sea level along the river course. From this profile, the length of the river was determined to be $L=84744$ m, the elevation of the source above sea level was found to be $z_{max}=1850$ m, and the average slope of the river was calculated as $\mathrm{Sl}=z_{max}/L=0.0218$.

For detailed modeling of the hydraulic jump, we chose a short river reach ($L=12.19$ m), where the Mzymta River transitions from a steep slope (${\partial }_{x}z=-0.06$) to a mild slope (${\partial }_{x}z=-0.01$), as shown in Figure 4 (right).

The errors in the bottom profile can be observed in Figure 4: in the lower reaches of the Mzymta River (60–75 km from the source), in the region of the Akhshtyr and Atskhu canyons, sharp increases in the bed elevation (up to 65 meters) are visible, which are physically unrealistic.

To solve the problem of the non-monotonic profile, we performed data smoothing using the smoothing spline method [78,79]. The constructed spline not only provides a smoothed representation of the river bed profile but also allows us to readily obtain its derivative ${\partial }_{x}z$ along the x coordinate ot the river. This derivative is required in equation (2). The smoothness of the spline ensures that the derivative is continuous and free from unphysical oscillations, which is essential for stable numerical simulations.

4.2. The Time Series of Water Discharge and Water Level

As input boundary conditions for discharge and water level, we used observational data provided by Roshydromet (the Russian Federal Service for Hydrometeorology and Environmental Monitoring). Figure 5 summarizes the time series collected from four hydrological gauging stations located along the Mzymta River. Rosa Khutor, Krasnaya Polyana, Gorki-Gorod are located in the higher reaches at an altitude more than 450 meters above sea level. Kazachiy Brod is located in the lower reaches of the river at an elevation of approximately 68 meters above sea level.

From these data, we estimated the characteristic low-flow discharge $Q=20$ ${\mathrm{m}}^3$/s, the characteristic low-flow depth $H=0.3$ m, and also calculated the characteristic discharge per unit width $q=0.4$ ${\mathrm{m}}^2$/s, the velocity $U=1.32$ m/s, and the Manning coefficient $n=0.05$ s/${\mathrm{m}}^{1/3}$.

5. Results

The simulations were performed using the DeepXDE library, a famous open-source framework for PINNs [36,49]. We used the computational resources of the server in Sirius University with an NVidia GPU A100. The training time for the PINNs typically amounts to 2-12 hours in the unsteady case. For the computations, we used a fully connected neural network with a width of 64 neurons and a depth of up to 6 layers. For smooth problems, 3 layers are sufficient, whereas problems involving hydraulic jumps are better approximated by deeper networks. We used the Swish (SiLU) activation function, as it has been shown to perform well for problems with sharp gradients and discontinuities. For the steady-state formulation, we used 1000 collocation points for the interior and 1 point on each boundary, whereas the unsteady formulation required 40,000 interior points, along with 1000 points on each boundary and 1000 initial condition points. The study involved selecting the optimal number of collocation points. For the steady-state case, the number of points varied from 100 to 2000. For the unsteady-state case, the number of points varied from 10,000 to 100,000. The final choice was based on achieving the required accuracy and reasonable calculation time. The Adam optimizer was chosen for its computational efficiency and robustness, as it is a gradient-based method well-suited for navigating the complex, non-convex loss landscapes typical of PINNs, avoiding local minima and exhibiting resilience to noise and sparse gradients. The learning rate was initially set to ${10}^{-4}$ and later decreased to ${10}^{-6}$ for fine-tuning of the results. The coefficients $\lambda$ for the boundary and initial conditions were set to unity, while the weights for the governing equations were adaptively adjusted across a wide range (spanning from ${10}^{-8}$ to ${10}^1$) to balance the contribution of each term to the total loss. All the numerical parameters of the model are collected in Table 2.

The analysis of the neural network hyperparameters was presented by the authors in their earlier study [53].

5.1. Full-Scale Stationary Flow

In this section, a steady-state flow analysis is performed for the entire Mzymta River, from its source to its mouth (L=84744 m).

For this scale, the friction parameter is large ($L/H\mathrm{Sl}=6145$), and the river tends toward normal flow conditions, i.e. $h\approx h_n$. The only location where a slight deviation from the normal depth profile is observed is the mild-slope reach near the river mouth. This is because in the mild-slope region, the depth increases and the velocity decreases, which in turn reduces the contribution of the term $s-s_{\mathrm{f}}$ in the momentum equation (2).

Under the given flow parameters (Table 1) and the bed slope profile (Figure 4, left), both supercritical and subcritical regimes occur, thus requiring boundary conditions for the water depth h at both the inlet and the outlet. Due to the stiffness of the problem, we set the depth h equal to the normal depth $h_n$.

Figure 6 presents the profiles of the critical depth $h_c$ and the normal depth $h_n$, along with the depth profiles obtained using the PINNs and the explicit fourth-order Runge–Kutta method (RK4) for the backwater equation (4).

For the specified parameters, the profiles of the critical depth $h_c$ and the normal depth $h_n$ intersect, and the river changes its flow regime several times, from subcritical to supercritical and back.

In this regard, solving the backwater equation (4) using classical methods (e.g., RK4) presents a difficulty. On supercritical reaches, the equation is Lyapunov-stable when solved from upstream to downstream, and the normal depth acts as an attractor of the solution. On subcritical reaches, conversely, the equation is stably integrated only in the upstream direction. When attempting to integrate the equation directly, the method remains stable until the critical depth is reached, after which it becomes unstable and highly sensitive to small perturbations in the solution. To compute the depth profile, we performed separate calculations for each reach, starting from the locations where the flow regime changes (Figure 6).

The results of numerical simulations demonstrate excellent agreement between the analytical assessment of the solution, the classical RK4 method, and the PINNs approach (with a maximum relative error of approximately 1%). Moreover, the PINNs operate stably even under multiple transitions between subcritical and supercritical regimes and vice versa.

5.2. Full-Scale Non-Stationary Flow

In this section, we perform modeling of unsteady flow in the Mzymta River. As initial conditions, we take the steady-state flow obtained in Section 5.1 and perturb the river discharge by a sinusoidal perturbation:

$$q_{\mathrm{in}}=q_0+q_{A}sin\left(\omega t\right).$$

(15)

Figure 7 shows the spatial-temporal distributions of the flow depth h and discharge Q, along with the characteristic lines $u+\sqrt{g/h}$, which correspond to the propagation of the inlet perturbation. The depth profile is observed to be advected downstream at the characteristic speed $u+\sqrt{g/h}$, indicating a finite propagation velocity rather than an instantaneous response to the upstream forcing.

Figure 8 shows the distribution of the flow depth h and Froude mask in the $(x,t)$-plane (left) and its projection onto the x-axis (right). The left panel displays oscillations in the water depth induced by the unsteady inlet discharge, as well as multiple transitions from subcritical to supercritical flow regimes. The right panel reveals that the depth profile remains close to the normal depth profile $h_n$ from the steady-state solution, although it is subject to perturbations that propagate downstream as the inlet discharge varies in time.

5.3. Stationary Hydraulic Jump Modeling

In the previous Section 5.1 and Section 5.2, it was shown that at the scale of the entire river, friction smoothes out all perturbations, and the river profile becomes close to normal flow. Therefore, to demonstrate the capability of PINNs to compute a hydraulic jump, we considered a short river reach where the bed slope flattens (Figure 4, right), thus promoting a transition from supercritical to subcritical flow. Due to the flattening of the bed slope, the normal depth increases sharply, leading to the formation of a backwater curve.

Figure 9 presents the depth profiles obtained using RK4 method and the PINNs, along with the critical depth $h_c$ and normal depth $h_n$ profiles. At the inlet and outlet, we imposed depth boundary conditions perturbed from the normal depth: in the upstream reach $\tilde{h}\left(0\right)<{\tilde{h}}_n\left(\tilde{s}\left(0\right)\right)$, and in the downstream reach $\tilde{h}\left(1\right)<{\tilde{h}}_n\left(\tilde{s}\left(1\right)\right)$.

To solve the problem using the Runge–Kutta method, it is necessary to separately integrate the backwater equation (4) downstream in the supercritical region and upstream in the subcritical region. In the middle of the domain, the reaches are matched via a hydraulic jump. The RK4 method requires manual splitting of the domain into subdomains according to the flow regime and manual matching of the corresponding reaches.

The PINNs method without regularization (red dashed curve) automatically produces a smoothed transition between the two regimes. However, applying regularization in the region where fr is close to unity enables the model to capture a sharp hydraulic jump.

It can be observed that in smooth regions, the PINNs and RK4 results are in good agreement.

5.4. Non-Stationary Hydraulic Jump Modeling

In this section we consider the non-stationary sinusoidal perturbation of the hydraulic jump problem, demonstrated in the previous Section 5.3.

Figure 10 presents the distributions of the flow depth h and discharge Q over the coordinates $(x,t)$ in the unsteady problem with a hydraulic jump. The colors indicate the regions of subcritical and supercritical flow.

Figure 11 presents the comparison of non-regularized and regularized solutions. It can be observed that the steady-state and unsteady flow profiles are similar for both the regularized and non-regularized cases. The main visible differences are a steeper hydraulic jump profile and larger amplitude of depth fluctuations in the upstream reach for the regularized solution.

6. Discussion and Future Plans

The main challenges we encountered in modeling hydrological processes using PINNs were:

• Balancing the weights in the loss function. Different stages of training require different weight ratios, and manual balancing becomes tedious. Therefore, we attempted to automate this process.
• Setting initial and boundary conditions, especially in the stiff full-scale formulation, where friction dominates the momentum equation and the system rapidly relaxes to normal flow, while the Saint-Venant equations degenerate into simplified analogs. In this case, specifying boundary and initial conditions inconsistent with normal flow leads to a sharp jump in the solution at the boundary. In this context, PINNs allow us to solve the problem without prescribing initial conditions at all, autonomously determining the initial distributions of depth and discharge that best match the near-normal flow.
• Modeling flow regime transitions and hydraulic jumps presents a significant difficulty for PINNs, as these networks are inherently designed to approximate smooth functions. To capture the discontinuity associated with a hydraulic jump, we introduced regularization terms into the loss function, which effectively inform the network about the presence and nature of the shock.

Future plans for model development:

1.

For the Mzymta River, one of the most challenging problems is the spatial variability of roughness $n=n\left(x\right)$. One can make the parameter $n\left(x\right)$ trainable, using known discharge data $Q(x,t)$, cross-sectional area $A(x,t)$, and minimizing the solution error with PINNs. This would yield a map of the effective channel roughness.

2.

The most promising formulation appears to be the following: a PINNs for the joint identification of $Q(x,t)$, $n\left(x\right)$, and $A(h,x)$ from water level data $h(x,t)$, discharge measurements $Q(x,t)$ and the Saint-Venant equations. In this framework, the network simultaneously recovers the hydrodynamics, roughness, and hydraulic characteristics of the channel. For the Mzymta River, this would no longer be merely a surrogate model, but a tool for extracting hidden parameters of the river system.

3.

If a long-term archive of hydrographs for the Mzymta River is available, instead of a standard PINNs, one could train a Physics-Informed DeepONet: $Q_{\mathrm{in}}\left(t\right)⟶h(x,t),Q(x,t)$. This would yield an operator model capable of instantly predicting the river’s response to a new flood event. For mountain rivers with rapid flood wave formation, this is likely the most promising research direction for the coming years.

4.

A two-dimensional (2D) formulation could also be explored to capture the spatial variability of the flow across the Mzymta River’s channel width. While the 1D model is sufficient for the main channel, a 2D approach would be necessary in regions with significant lateral gradients, such as tributary confluences, anabranches, or in wide floodplain areas. However, this would considerably increase the computational cost and data requirements.

7. Conclusion

The novelty of the work lies in the study of hydrological processes for a mountain river with a given relief profile using the PINNs method.

A dynamic weighting strategy was employed for the loss function in the PINNs framework. The coefficients were updated at every training epoch to ensure a proper balance among the various loss components, thereby preventing any individual term from exerting a dominant influence on the optimization.

The full-scale simulations demonstrated that the flow rapidly approaches a normal-depth profile consistent with the prescribed bed topography, Manning coefficient, and discharge. This behavior is explained by the rapid damping of any local perturbations, which become negligible at the river scale. Non-stationary simulations showed that the flow profile is slightly perturbed relative to the initial state as the inlet discharge varies. The perturbations propagate downstream along the characteristic curves.

A particularly challenging and interesting problem in river hydraulics is the regime transition from supercritical to subcritical flow, which occurs through a hydraulic jump. In the present study, we show that PINNs are capable of accurately capturing such discontinuous flow features. For this purpose, we selected a short reach of the Mzymta River characterized by a flattening of the bed slope, where a backwater profile develops and connects to the upstream supercritical flow through a hydraulic jump.

To accurately simulate hydraulic jumps, we introduced a regularization strategy for the loss function. The approach relies on a shock-detector function, which identifies potential jump regions based on the local Froude number (close to unity) and the gradient of the water depth. Within the detected jump zone, the weight of the momentum equation in the loss is reduced. Furthermore, to sharpen the transition, we add a regularization term that penalizes the spatial extent of the detected jump area.

The PINNs solutions were verified by comparing with the classical Runge–Kutta method. Verification showed that the results obtained using the two methods are in excellent agreement.