Implicit Finite-Difference Scheme for Two-Dimensional Flood Modelling Using Shallow-Water Equations

1. Introduction

The current stage of development in Earth sciences and engineering hydrology is characterised by a profound transformation in approaches to studying water systems [1,2,3]. Simplified, static assessment methods are being replaced by complex, dynamic models that can reproduce the evolution of natural and anthropogenic processes across space and time. In this context, the shallow water equation system occupies a unique position as a fundamental theoretical basis, ensuring a balance between the physical completeness of the description of the phenomenon and computational feasibility [1,4,5,6]. However, a true revolution in the capabilities of analysis and forecasting is taking place at the intersection of this mathematical theory with geoinformation system technologies. The integration of shallow water equations and geoinformation systems forms a new research paradigm, in which abstract hydrodynamic calculations acquire spatial specificity, transforming them into tools for management decision-making, risk assessment, and engineering design.

The mathematical model of the system in conservative form for modelling the processes investigated in this work takes the form [2,3]:

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial H}{\partial t}}+{\displaystyle \frac{\partial \left(Hu\right)}{\partial x}}+{\displaystyle \frac{\partial \left(Hv\right)}{\partial y}}=0,\\ {\displaystyle \frac{\partial \left(Hu\right)}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left(H{u}^2+{\displaystyle \frac{g}{2}}{H}^2\right)+{\displaystyle \frac{\partial \left(Huv\right)}{\partial y}}=0,\\ {\displaystyle \frac{\partial \left(Hv\right)}{\partial t}}+{\displaystyle \frac{\partial \left(Huv\right)}{\partial x}}+{\displaystyle \frac{\partial }{\partial y}}\left(H{v}^2+{\displaystyle \frac{g}{2}}{H}^2\right)=0,\end{array}\right.$$

(1)

where H(x,y,t) is the total water depth; u(x,y,t) is the velocity along x; v(x,y,t) is the velocity along y; g is the acceleration due to gravity.

The shallow water equations represent a class of hyperbolic partial differential equations that describe the motion of a fluid with a free surface in a gravitational field. Their derivation from the full three-dimensional Navier–Stokes equations via a depth-averaging procedure is based on the key physical assumption of a hydrostatic pressure distribution, which holds for the vast majority of geophysical flows in which the horizontal scales significantly exceed the vertical ones [3,7]. This simplification is not a shortcoming, but rather a methodological tool that allows secondary vertical accelerations to be filtered out, focusing on the dominant horizontal transport of mass and momentum.

The relevance of these equations in modern research stems from their ability to reproduce a wide spectrum of wave motions critical to understanding the dynamics of the oceans, seas and planetary atmospheres. In particular, they accurately describe high-frequency inertial-gravitational waves (Poincare waves), which are responsible for the rapid redistribution of energy within a basin, and low-frequency Rossby waves, whose dynamics determine synoptic variability and large-scale circulation [7,8,9]. It is precisely this ability to model multi-scale processes that makes shallow-water equations and their modifications indispensable tools in oceanography, limnology, and river hydraulics for the study of inundations, floods, and tsunami surges onto coasts (Figure 1).

Studies aimed at accounting for the medium’s variable density and compressibility are of particular importance [7,10,11]. The development of ‘compressible shallow water equations’ models paves the way for modelling fundamentally new classes of phenomena. This concerns not only dust storms on Mars or two-phase flows under terrestrial conditions (for example, the transport of sediments or snow-dust masses), but also fundamental processes in astrophysical discs, where gravitational interaction and compressibility play a decisive role.

Another revolutionary direction is the synergy of classical hydrodynamic models with machine learning methods [12,13,14]. The current trend lies in the creation of so-called “universal” or “data-driven hydrodynamic models” [13,15]. Applying differential programming to shallow-water equations allows not only for solving equations with pre-defined parameters, but also for “training” the model based on available observations. This makes it possible, for example, to identify the spatial distribution of the Manning roughness coefficient with unprecedented accuracy, adapting the model to the local characteristics of the river channel and floodplain, which is critical for reliable flood zone forecasting [15,16,17].

The greatest practical potential of shallow water equations is realised precisely through their close integration with geographic information systems [18]. Such a system ceases to be merely a tool for visualising final results. It becomes an integral component of the entire modelling cycle – from the preparation of input data to the interpretation of the results obtained.

The foundation of any hydrodynamic calculation is a digital elevation model, which serves as the spatial basis for constructing the computational grid [16,19,20]. It is through the use of geographic information systems that hydrologically correct processing of the digital elevation model is carried out – removal of artefacts, “filling” of closed depressions, and calculation of surface runoff directions. Next, using geoprocessing tools, raster and vector layers containing information on spatially distributed parameters—such as substrate types, vegetation cover and soils—are overlaid onto the study area. Each of these types is assigned numerical values for hydraulic resistance (roughness coefficients), which allows the model to respond appropriately to landscape heterogeneity (Figure 2) [20,21].

The model verification process also takes on a new dimension thanks to modern geoinformation systems. Modelling results (velocity fields, depths, free-surface levels) are imported back into the GIS environment for comparison with Earth remote sensing data, ground-based hydrological stations, or historical flood maps [10,16]. Such spatial analysis allows not only a visual assessment of the model’s adequacy, but also the calculation of quantitative accuracy metrics.

It should be noted that the implementation of the shallow water equation system using explicit schemes, despite their high speed, has significant limitations on the time step due to the Courant–Friedrichs–Lewy condition, which is given by [1]:

$$\mathsf{\Delta }t\le CFL\cdot \min \left({\displaystyle \frac{\mathsf{\Delta }x}{\left|u_{i,j}\right|+\sqrt{g{H}_{i,j}}}},{\displaystyle \frac{\mathsf{\Delta }y}{\left|v_{i,j}\right|+\sqrt{g{H}_{i,j}}}}\right).$$

The value of the parameter CFL≈0.5 for 2D explicit schemes; CFL ≤ 1 for explicit 1D schemes and for Lax-Friedrichs-type schemes – typically 0.3–0.5.

The synthesis of hydrodynamic models based on shallow-water equations and the analytical capabilities of geographic information systems creates a powerful toolkit for solving a wide range of applied problems in national security, the economy, and ecology [22,23,24]. In the field of natural disaster risk management, a key task is the creation of flood hazard and risk maps. The application of these equations enables the modelling of the dynamics of flood waves, the passage of catastrophic floods with varying return periods, and storm surges. Integrating these scenarios into geographic information systems enables quantification of the population at risk, identification of critical infrastructure at risk (bridges, hospitals, chemical plants), and calculation of expected economic losses. This provides a scientifically sound basis for spatial planning, insurance and the development of evacuation plans.

In environmental terms, this combination of technologies enables the prediction of pollutant migration in water bodies. By modelling currents using such equations and superimposing contaminant movement trajectories onto them, it is possible to predict with high accuracy the zones of pollutant sedimentation, including heavy metals, petroleum products, and radionuclides [8,14]. This approach is indispensable for assessing the impact of accidental discharges, designing water protection zones, and justifying remediation measures for disturbed areas.

2. Analysis of Shallow Water Equation Modelling Methods

The choice of method depends on the specific physics of the process and engineering requirements. Recent comparative studies allow the following generalisations.

For smooth flows, high-order methods and compact finite-difference methods yield the best results. For flows with discontinuities (e.g., hydraulic jumps), the finite volume method (particularly the counter-flow scheme) is the natural choice due to its conservativeness and ability to localise the wave front accurately. For example, the abstract of one of the studies notes that the developed high-precision scheme significantly outperforms other methods in accuracy, reducing error by several orders of magnitude. A comparison of the above methods has shown that both are reliable for quasi-steady states that account for friction, but require different approaches to approximating the source terms. A critical requirement for modern schemes is the property of balance. The finite element and finite volume methods achieve this through a special approximation of the source term, such as the seabed topography. For the finite volume method, this means that the scheme must accurately solve the given problem in the corresponding mathematical formulation [1].

Traditional finite-difference methods on structured grids are the fastest. Finite-volume methods offer a good balance between speed and accuracy. The finite-element method and the Galerkin method are more computationally expensive, but their advantage in accuracy on complex grids may offset this drawback. A study of methods for the local inertia equation showed that implicit methods have the highest stability and the smallest error when modelling the water-land boundary, whereas the exponential method can produce oscillations. A summary comparison of the flood modelling methods considered that use shallow-water equations is given in Table 1.

Traditional finite-difference schemes on regular grids. These are historically the first and most extensively studied approaches. They are based on approximating derivatives with finite-difference analogues on structured grids. This category includes such well-known schemes as the Lax–Friedrichs, Lax–Wendroff, and McCormack schemes, among others, which are described in detail in the classical literature. The high computational speed of these methods is their main advantage. Thanks to the simplicity of the data structures and the regularity of the computations, they lend themselves excellently to vectorisation and parallelisation, making them the fastest of all methods. The basic schemes are relatively easy to programme, making them attractive for rapid prototyping. As a recent comparative study has shown, one of the traditional finite-difference schemes demonstrated the best accuracy in modelling large-scale ocean currents and storm surges on the shelf.

It is also worth noting the shortcomings of such schemes. One such shortcoming is the poor approximation of complex geometries. The use of rectangular grids yields a ‘stepped’ approximation of the coastline, leading to significant errors, particularly in complex coastal areas. Classic high-order schemes (e.g., Lax–Wendroff) generate undesirable oscillations (discontinuities) near hydraulic jumps. First-order schemes (Lax–Friedrichs) are monotonic but have high scheme viscosity, which ‘smears’ sharp wave profiles. Such schemes require special procedures to ensure the exact conservation of the state of rest over an uneven seabed. To overcome the geometry problem, modifications have been developed to better represent curved boundaries on a rectangular grid. For example, the technique of introducing diagonal segments in addition to vertical and horizontal ones to denote the ‘water-land’ boundary. In this case, significantly higher accuracy is achieved in coastal zones than with the traditional stepwise approach, without the need for excessive grid refinement. This approach yields more accurate results without a significant increase in the number of grid nodes, thereby reducing computational time. Of course, it should be noted that programme implementation is more complex – it requires modifying classical algorithms and implementing more complex programme logic. However, this approach still lags behind unstructured grids in terms of flexibility, which are used in the finite element and finite difference methods.

Mesh-free generalised finite difference methods represent a modern direction in the development of classical difference methods, designed to overcome the limitations of structured meshes. These modifications enable the construction of difference approximations with irregularly spaced nodes, significantly increasing the method’s flexibility. It is worth noting the flexibility of such methods, as they can handle complex geometries without requiring a grid, which is their main advantage over classical finite-difference methods. In such methods, high accuracy can be achieved on smooth solutions. There are spatio-temporal versions of generalised difference methods that demonstrate high efficiency and stability, for example, for the Korteweg–de Vries and regularised long-wave equations, which are models of dispersive waves in shallow water. Of course, the selection of neighbouring nodes for approximation and the construction of system matrices for each point are considerably more computationally expensive than simple indices in a structured grid. The method is considerably more complex to programme than the classical finite-difference method. It is not as well-established a standard for shallow-water equations as the finite-volume method.

Finite difference method for dispersive shallow-water equations. A distinct class is formed by difference schemes developed to solve non-classical, dispersive shallow-water models (of the Businescu type) that describe wave propagation while accounting for frequency dispersion. This method and its modifications allow the modelling of phenomena that are inaccessible to classical, non-dispersive shallow-water equations, such as the transformation of a tsunami in deep water. Research shows that the stability conditions for finite-difference schemes approximating dispersive equations are weaker than for schemes approximating classical shallow-water equations. However, in this case, the analysis of the dispersive and dissipative properties of such schemes is considerably more complex. It should also be noted that, for two-dimensional cases, there is the problem of the scheme’s phase error depending on the direction of wave propagation, which requires the development of special, rotation-invariant schemes.

Simplified finite-difference schemes for approximating the kinematic wave. For very simple configurations, the simplified form of the Saint-Venant equations – the kinematic wave equations – is used instead of the full equations. These are the simplest of all methods, very fast and stable under certain conditions. These methods do not work in the presence of a backwater, in the presence of changes in flow direction, or on complex terrain with steep gradients. They are unsuitable for most engineering problems on rivers. A comparative table of the main known difference methods used in oceanography is given below (Table 2).

3. Proposed Method

To construct an implicit difference scheme for (1), we shall first construct an explicit difference scheme, the principle of which will also be used for the implicit difference scheme. Let us assume we have a structured rectangular grid with nodes of the type $\left(i,j,k\right)$, where

$$x_i=i\cdot h_x,h_x={\displaystyle \frac{b-a}{N_x}},i=\overline{0;N_x+1};$$

$$y_i=j\cdot h_y,h_y={\displaystyle \frac{d-c}{N_y}},j=\overline{0;N_y+1};$$

$$t_k=k\cdot h_t,h_t={\displaystyle \frac{T_{finish}}{N_t}},k=\overline{0;N_t}.$$

All subsequent calculations will form a system of equations in discrete form for the internal nodes of the computational domain, i.e., for$i=\overline{1;N_x},j=\overline{1;N_y},k=\overline{1;N_k}.$ In this case, the discrete analogue of the first equation of the shallow water equation system for any internal node$\left(i,j,k\right)$ can be expressed as follows:

$${\displaystyle \frac{H_{i,j}^{k+1}-H_{i,j}^k}{\mathsf{\Delta }t}}=-{\displaystyle \frac{{\left(Hu\right)}_{i+\frac{1}{2},j}^k-{\left(Hu\right)}_{i-\frac{1}{2},j}^k}{\mathsf{\Delta }x}}-{\displaystyle \frac{{\left(Hv\right)}_{i,j+\frac{1}{2}}^k-{\left(Hv\right)}_{i,j-\frac{1}{2}}^k}{\mathsf{\Delta }y}}.$$

(2)

It is clear that in (2) nodes of type

$$\left(i+0.5,j,k\right),\left(i-0.5,j,k\right),\left(i,j+0.5,k\right),\left(i,j-0.5,k\right)$$

have non-integer indices. In this case, it is possible and even recommended to use an approximate representation via average values, i.e., via the expressions:

$${\left(Hu\right)}_{i+{\displaystyle \frac{1}{2}},j}^k=H_{i+{\displaystyle \frac{1}{2}},j}^k{u}_{i+{\displaystyle \frac{1}{2}},j}^k\approx {\displaystyle \frac{H_{i+1,j}^k{u}_{i+1,j}^k+H_{i,j}^k{u}_{i,j}^k}{2}};$$

$${\left(Hu\right)}_{i-{\displaystyle \frac{1}{2}},j}^k=H_{i-{\displaystyle \frac{1}{2}},j}^k{u}_{i-{\displaystyle \frac{1}{2}},j}^k\approx {\displaystyle \frac{H_{i-1,j}^k{u}_{i-1,j}^k+H_{i,j}^k{u}_{i,j}^k}{2}};$$

(3)

$${\left(Hv\right)}_{i,j+{\displaystyle \frac{1}{2}}}^k=H_{i,j+{\displaystyle \frac{1}{2}}}^k{v}_{i,j+{\displaystyle \frac{1}{2}}}^k\approx {\displaystyle \frac{H_{i,j+1}^k{v}_{i,j+1}^k+H_{i,j}^k{H}_{i,j}^k}{2}};$$

$${\left(Hv\right)}_{i,j-{\displaystyle \frac{1}{2}}}^k=H_{i,j-{\displaystyle \frac{1}{2}}}^k{v}_{i,j-{\displaystyle \frac{1}{2}}}^k\approx {\displaystyle \frac{H_{i,j-1}^k{v}_{i,j-1}^k+H_{i,j}^k{H}_{i,j}^k}{2}}.$$

Let us substitute expressions (3) into (2) for the first equation of the shallow water equation system (1). As a result, we obtain:

$$H_{i,j}^{k+1}=H_{i,j}^k-\mathsf{\Delta }t\left(\begin{array}{c}{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(\begin{array}{c}{\displaystyle \frac{H_{i+1,j}^k{u}_{i+1,j}^k+H_{i,j}^k{u}_{i,j}^k}{2}}-{\displaystyle \frac{H_{i-1,j}^k{u}_{i-1,j}^k+H_{i,j}^k{u}_{i,j}^k}{2}}\end{array}\right)+\\ +{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(\begin{array}{c}{\displaystyle \frac{H_{i,j+1}^k{v}_{i,j+1}^k+H_{i,j}^k{H}_{i,j}^k}{2}}-{\displaystyle \frac{H_{i,j-1}^k{v}_{i,j-1}^k+H_{i,j}^k{H}_{i,j}^k}{2}}\end{array}\right)\end{array}\right).$$

(4)

The final finite difference scheme for the discrete version of the first equation of the shallow water equation system of the form (1):

$$H_{i,j}^{k+1}=H_{i,j}^k-{\displaystyle \frac{\mathsf{\Delta }t}{2}}\left(\begin{array}{c}{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(\begin{array}{c}{H}_{i+1,j}^k{u}_{i+1,j}^k-H_{i-1,j}^k{u}_{i-1,j}^k\end{array}\right)+\\ +{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(H_{i,j+1}^k{v}_{i,j+1}^k-H_{i,j-1}^k{v}_{i,j-1}^k\right)\end{array}\right).$$

(5)

In expression (5), the central element (node (i, j, k+1)) is expressed in terms of nodes of the form (*, *, k).

Similarly, the discrete analogue of the second equation of the shallow water equation system, presented in conservative form (1), is formulated. We write it as follows:

$${\displaystyle \frac{H_{i,j}^{k+1}{u}_{i,j}^{k+1}-H_{i,j}^k{u}_{i,j}^k}{\mathsf{\Delta }t}}=$$

$$=-{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(\begin{array}{c}{H}_{i+{\displaystyle \frac{1}{2}},j}^k{\left(u_{i+{\displaystyle \frac{1}{2}},j}^k\right)}^2+{\displaystyle \frac{g}{2}}{\left(H_{i+{\displaystyle \frac{1}{2}},j}^k\right)}^2-\\ -H_{i-{\displaystyle \frac{1}{2}},j}^k{\left(u_{i-{\displaystyle \frac{1}{2}},j}^k\right)}^2-{\displaystyle \frac{g}{2}}{\left(H_{i-{\displaystyle \frac{1}{2}},j}^k\right)}^2\end{array}\right)-{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(\begin{array}{c}{H}_{i,j+{\displaystyle \frac{1}{2}}}^k{u}_{i,j+{\displaystyle \frac{1}{2}}}^k{v}_{i,j+{\displaystyle \frac{1}{2}}}^k-\\ -H_{i,j-{\displaystyle \frac{1}{2}}}^k{u}_{i,j-{\displaystyle \frac{1}{2}}}^k{v}_{i,j-{\displaystyle \frac{1}{2}}}^k\end{array}\right).$$

(6)

In expression (6), the expressions for the central difference derivatives are taken into account:

${\left.{\displaystyle \frac{\partial }{\partial x}}\left(H{u}^2+{\displaystyle \frac{g}{2}}{H}^2\right)\right|}_{i,j}^k\approx {\displaystyle \frac{{\left.\left(H{u}^2+\frac{g}{2}{H}^2\right)\right|}_{i+\frac{1}{2},j}^k-{\left.\left(H{u}^2+\frac{g}{2}{H}^2\right)\right|}_{i-\frac{1}{2},j}^k}{\mathsf{\Delta }x}};$${\left.{\displaystyle \frac{\partial \left(Huv\right)}{\partial y}}\right|}_{i,j}^k\approx {\displaystyle \frac{{\left.\left(Huv\right)\right|}_{i,j+\frac{1}{2}}^k-{\left.\left(Huv\right)\right|}_{i,j-\frac{1}{2}}^k}{\mathsf{\Delta }y}};$

$${\left.\left(H{u}^2+{\displaystyle \frac{g}{2}}{H}^2\right)\right|}_{i+{\displaystyle \frac{1}{2}},j}^k=H_{i+{\displaystyle \frac{1}{2}},j}^k{\left(u_{i+{\displaystyle \frac{1}{2}},j}^k\right)}^2+{\displaystyle \frac{g}{2}}{\left(H_{i+{\displaystyle \frac{1}{2}},j}^k\right)}^2;$$

(7)

${\left.\left(H{u}^2+{\displaystyle \frac{g}{2}}{H}^2\right)\right|}_{i-{\displaystyle \frac{1}{2}},j}^k=H_{i-{\displaystyle \frac{1}{2}},j}^k{\left(u_{i-{\displaystyle \frac{1}{2}},j}^k\right)}^2+{\displaystyle \frac{g}{2}}{\left(H_{i-{\displaystyle \frac{1}{2}},j}^k\right)}^2;$${\left.\left(Huv\right)\right|}_{i,j+{\displaystyle \frac{1}{2}}}^k=H_{i,j+{\displaystyle \frac{1}{2}}}^k{u}_{i,j+{\displaystyle \frac{1}{2}}}^k{v}_{i,j+{\displaystyle \frac{1}{2}}}^k;$${\left.\left(Huv\right)\right|}_{i,j-{\displaystyle \frac{1}{2}}}^k=H_{i,j-{\displaystyle \frac{1}{2}}}^k{u}_{i,j-{\displaystyle \frac{1}{2}}}^k{v}_{i,j-{\displaystyle \frac{1}{2}}}^k.$For expression (7), we finally obtain the expression for $u_{i,j}^{k+1}$:

$u_{i,j}^{k+1}={\displaystyle \frac{H_{i,j}^k}{H_{i,j}^{k+1}}}{u}_{i,j}^k-{\displaystyle \frac{1}{H_{i,j}^{k+1}}}\cdot {\displaystyle \frac{\mathsf{\Delta }t}{\mathsf{\Delta }x}}\left(\begin{array}{c}\left(\begin{array}{c}{\displaystyle \frac{H_{i+1,j}^k{\left(u_{i+1,j}^k\right)}^2+H_{i,j}^k{\left(u_{i,j}^k\right)}^2}{2}}+\\ +{\displaystyle \frac{g\left({\left(H_{i+1,j}^k\right)}^2+{\left(H_{i,j}^k\right)}^2\right)}{4}}\end{array}\right)-\\ -\left(\begin{array}{c}{\displaystyle \frac{H_{i-1,j}^k{\left(u_{i-1,j}^k\right)}^2+H_{i,j}^k{\left(u_{i,j}^k\right)}^2}{2}}+\\ +{\displaystyle \frac{g\left({\left(H_{i-1,j}^k\right)}^2+{\left(H_{i,j}^k\right)}^2\right)}{4}}\end{array}\right)\end{array}\right)-$

$$-{\displaystyle \frac{\mathsf{\Delta }t}{\mathsf{\Delta }y}}\left(\begin{array}{c}{\displaystyle \frac{H_{i,j+1}^k{u}_{i,j+1}^k{v}_{i,j+1}^k+H_{i,j}^k{u}_{i,j}^k{v}_{i,j}^k}{H_{i,j}^{k+1}}}-\\ -{\displaystyle \frac{H_{i,j-1}^k{u}_{i,j-1}^k{v}_{i,j-1}^k+H_{i,j}^k{u}_{i,j}^k{v}_{i,j}^k}{H_{i,j}^{k+1}}}\end{array}\right).$$

(8)

In (8), just as in (4), the following was taken into account:

$H_{i+{\displaystyle \frac{1}{2}},j}^k{\left(u_{i+{\displaystyle \frac{1}{2}},j}^k\right)}^2\approx {\displaystyle \frac{H_{i+1,j}^k{\left(u_{i+1,j}^k\right)}^2+H_{i,j}^k{\left(u_{i,j}^k\right)}^2}{2}};$$H_{i-{\displaystyle \frac{1}{2}},j}^k{\left(u_{i-{\displaystyle \frac{1}{2}},j}^k\right)}^2\approx {\displaystyle \frac{H_{i-1,j}^k{\left(u_{i-1,j}^k\right)}^2+H_{i,j}^k{\left(u_{i,j}^k\right)}^2}{2}};$${\left(Huv\right)}_{i,j+{\displaystyle \frac{1}{2}}}^k=H_{i,j+{\displaystyle \frac{1}{2}}}^k{u}_{i,j+{\displaystyle \frac{1}{2}}}^k{v}_{i,j+{\displaystyle \frac{1}{2}}}^k\approx {\displaystyle \frac{H_{i,j+1}^k{u}_{i,j+1}^k{v}_{i,j+1}^k+H_{i,j}^k{u}_{i,j}^k{v}_{i,j}^k}{2}};$${\left(Huv\right)}_{i,j-{\displaystyle \frac{1}{2}}}^k=H_{i,j-{\displaystyle \frac{1}{2}}}^k{u}_{i,j-{\displaystyle \frac{1}{2}}}^k{v}_{i,j-{\displaystyle \frac{1}{2}}}^k\approx {\displaystyle \frac{H_{i,j-1}^k{u}_{i,j-1}^k{v}_{i,j-1}^k+H_{i,j}^k{u}_{i,j}^k{v}_{i,j}^k}{2}}.$Thus, we finally have:

$u_{i,j}^{k+1}={\displaystyle \frac{H_{i,j}^k}{H_{i,j}^{k+1}}}{u}_{i,j}^k-{\displaystyle \frac{1}{H_{i,j}^{k+1}}}\cdot {\displaystyle \frac{\mathsf{\Delta }t}{2\mathsf{\Delta }x}}\left(\begin{array}{c}\left(H_{i+1,j}^k{\left(u_{i+1,j}^k\right)}^2+{\displaystyle \frac{g{\left(H_{i+1,j}^k\right)}^2}{2}}\right)-\\ -\left({\displaystyle \frac{H_{i-1,j}^k{\left(u_{i-1,j}^k\right)}^2}{2}}+{\displaystyle \frac{g{\left(H_{i-1,j}^k\right)}^2}{2}}\right)\end{array}\right)-$

$$-{\displaystyle \frac{1}{H_{i,j}^{k+1}}}\cdot {\displaystyle \frac{\mathsf{\Delta }t}{2\mathsf{\Delta }y}}\left(H_{i,j+1}^k{u}_{i,j+1}^k{v}_{i,j+1}^k+H_{i,j-1}^k{u}_{i,j-1}^k{v}_{i,j-1}^k\right).$$

(9)

The final step is to formulate a discrete analysis for the third equation of the shallow water equation system (1). We shall use central difference schemes to approximate the derivatives with respect to the spatial coordinates$x$ and$y$ . In this case, we obtain:

${\displaystyle \frac{H_{i,j}^{k+1}{v}_{i,j}^{k+1}-H_{i,j}^k{v}_{i,j}^k}{\mathsf{\Delta }t}}=-{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(\begin{array}{c}{H}_{i+{\displaystyle \frac{1}{2}},j}^k{u}_{i+{\displaystyle \frac{1}{2}},j}^k{v}_{i+{\displaystyle \frac{1}{2}},j}^k-\\ -H_{i-{\displaystyle \frac{1}{2}},j}^k{u}_{i-{\displaystyle \frac{1}{2}},j}^k{v}_{i-{\displaystyle \frac{1}{2}},j}^k\end{array}\right)-$

$$-{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(\begin{array}{c}{H}_{i,j+{\displaystyle \frac{1}{2}}}^k{\left(v_{i,j+{\displaystyle \frac{1}{2}}}^k\right)}^2+{\displaystyle \frac{g}{2}}{\left(H_{i,j+{\displaystyle \frac{1}{2}}}^k\right)}^2-\\ -H_{i,j-{\displaystyle \frac{1}{2}}}^k{\left(v_{i,j-{\displaystyle \frac{1}{2}}}^k\right)}^2-{\displaystyle \frac{g}{2}}{\left(H_{i,j-{\displaystyle \frac{1}{2}}}^k\right)}^2\end{array}\right).$$

(10)

Multiplying both sides of (10) by $\mathsf{\Delta }t$, and expanding the terms with indices of the form (*, *, k+1), we obtain:

$H_{i,j}^{k+1}{v}_{i,j}^{k+1}=H_{i,j}^k{v}_{i,j}^k-{\displaystyle \frac{\mathsf{\Delta }t}{\mathsf{\Delta }x}}\left(\begin{array}{c}{H}_{i+{\displaystyle \frac{1}{2}},j}^k{u}_{i+{\displaystyle \frac{1}{2}},j}^k{v}_{i+{\displaystyle \frac{1}{2}},j}^k-H_{i-{\displaystyle \frac{1}{2}},j}^k{u}_{i-{\displaystyle \frac{1}{2}},j}^k{v}_{i-{\displaystyle \frac{1}{2}},j}^k\end{array}\right)-$$-{\displaystyle \frac{\mathsf{\Delta }t}{\mathsf{\Delta }y}}\left(\begin{array}{c}{H}_{i,j+{\displaystyle \frac{1}{2}}}^k{\left(v_{i,j+{\displaystyle \frac{1}{2}}}^k\right)}^2+{\displaystyle \frac{g}{2}}{\left(H_{i,j+{\displaystyle \frac{1}{2}}}^k\right)}^2-H_{i,j-{\displaystyle \frac{1}{2}}}^k{\left(v_{i,j-{\displaystyle \frac{1}{2}}}^k\right)}^2-{\displaystyle \frac{g}{2}}{\left(H_{i,j-{\displaystyle \frac{1}{2}}}^k\right)}^2\end{array}\right).$Now we obtain the final expression for $v_{i,i}^{k+1}$. It takes the following form:

$v_{i,j}^{k+1}={\displaystyle \frac{H_{i,j}^k}{H_{i,j}^{k+1}}}{v}_{i,j}^k-{\displaystyle \frac{1}{H_{i,j}^{k+1}}}\cdot {\displaystyle \frac{\mathsf{\Delta }t}{2\mathsf{\Delta }x}}\left(\begin{array}{c}{H}_{i+1,j}^k{u}_{i+1,j}^k{v}_{i+1,j}^k-\\ -H_{i-1,j}^k{u}_{i-1,j}^k{v}_{i-1,j}^k\end{array}\right)-$

$$-{\displaystyle \frac{1}{H_{i,j}^{k+1}}}\cdot {\displaystyle \frac{\mathsf{\Delta }t}{2\mathsf{\Delta }y}}\left(\begin{array}{c}{H}_{i,j+1}^k{\left(v_{i,j+1}^k\right)}^2+{\displaystyle \frac{g\left({\left(H_{i,j+1}^k\right)}^2\right)}{2}}\\ -H_{i,j-1}^k{\left(v_{i,j-1}^k\right)}^2-{\displaystyle \frac{g\left({\left(H_{i,j-1}^k\right)}^2\right)}{2}}\end{array}\right).$$

(11)

Finally, we have three difference equations for the explicit difference scheme:

$H_{i,j}^{k+1}=H_{i,j}^k-{\displaystyle \frac{\mathsf{\Delta }t}{2}}\left(\begin{array}{c}{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(\begin{array}{c}{H}_{i+1,j}^k{u}_{i+1,j}^k-H_{i-1,j}^k{u}_{i-1,j}^k\end{array}\right)+\\ +{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(H_{i,j+1}^k{v}_{i,j+1}^k-H_{i,j-1}^k{v}_{i,j-1}^k\right)\end{array}\right);$$u_{i,j}^{k+1}={\displaystyle \frac{H_{i,j}^k}{H_{i,j}^{k+1}}}{u}_{i,j}^k-{\displaystyle \frac{1}{H_{i,j}^{k+1}}}\cdot {\displaystyle \frac{\mathsf{\Delta }t}{2\mathsf{\Delta }x}}\left(\begin{array}{c}\left(H_{i+1,j}^k{\left(u_{i+1,j}^k\right)}^2+{\displaystyle \frac{g{\left(H_{i+1,j}^k\right)}^2}{2}}\right)-\\ -\left({\displaystyle \frac{H_{i-1,j}^k{\left(u_{i-1,j}^k\right)}^2}{2}}+{\displaystyle \frac{g{\left(H_{i-1,j}^k\right)}^2}{2}}\right)\end{array}\right)-$

$$-{\displaystyle \frac{1}{H_{i,j}^{k+1}}}\cdot {\displaystyle \frac{\mathsf{\Delta }t}{2\mathsf{\Delta }y}}\left(H_{i,j+1}^k{u}_{i,j+1}^k{v}_{i,j+1}^k+H_{i,j-1}^k{u}_{i,j-1}^k{v}_{i,j-1}^k\right);$$

(12)

$v_{i,j}^{k+1}={\displaystyle \frac{H_{i,j}^k}{H_{i,j}^{k+1}}}{v}_{i,j}^k-{\displaystyle \frac{1}{H_{i,j}^{k+1}}}\cdot {\displaystyle \frac{\mathsf{\Delta }t}{2\mathsf{\Delta }x}}\left(\begin{array}{c}{H}_{i+1,j}^k{u}_{i+1,j}^k{v}_{i+1,j}^k-\\ -H_{i-1,j}^k{u}_{i-1,j}^k{v}_{i-1,j}^k\end{array}\right)-$$-{\displaystyle \frac{1}{H_{i,j}^{k+1}}}\cdot {\displaystyle \frac{\mathsf{\Delta }t}{2\mathsf{\Delta }y}}\left(\begin{array}{c}{H}_{i,j+1}^k{\left(v_{i,j+1}^k\right)}^2+{\displaystyle \frac{g\left({\left(H_{i,j+1}^k\right)}^2\right)}{2}}\\ -H_{i,j-1}^k{\left(v_{i,j-1}^k\right)}^2-{\displaystyle \frac{g\left({\left(H_{i,j-1}^k\right)}^2\right)}{2}}\end{array}\right).$We now proceed to derive an implicit difference scheme for (1). It will be based on the fundamental idea of the Crank-Nicolson scheme, which has proven to be quite effective in modelling complex systems and processes described by hyperbolic and parabolic equations with respect to time.

For the first equation, we obtain the following expression based on the Crank-Nicholson scheme:

${\displaystyle \frac{H_{i,j}^{k+1}-H_{i,j}^k}{\mathsf{\Delta }t}}={\displaystyle \frac{1}{2}}\left(\begin{array}{c}-{\displaystyle \frac{{\left(Hu\right)}_{i+\frac{1}{2},j}^k-{\left(Hu\right)}_{i-\frac{1}{2},j}^k}{\mathsf{\Delta }x}}-{\displaystyle \frac{{\left(Hv\right)}_{i,j+\frac{1}{2}}^k-{\left(Hv\right)}_{i,j-\frac{1}{2}}^k}{\mathsf{\Delta }y}}-\\ -{\displaystyle \frac{{\left(Hu\right)}_{i+\frac{1}{2},j}^{k+1}-{\left(Hu\right)}_{i-\frac{1}{2},j}^{k+1}}{\mathsf{\Delta }x}}-{\displaystyle \frac{{\left(Hv\right)}_{i,j+\frac{1}{2}}^{k+1}-{\left(Hv\right)}_{i,j-\frac{1}{2}}^{k+1}}{\mathsf{\Delta }y}}\end{array}\right).$Therefore, the discrete form of the first equation in (1) is

${\displaystyle \frac{H_{i,j}^{k+1}-H_{i,j}^k}{\mathsf{\Delta }t}}={\displaystyle \frac{1}{2}}\left(\begin{array}{c}-{\displaystyle \frac{H_{i+\frac{1}{2},j}^k{u}_{i+\frac{1}{2},j}^k-H_{i-\frac{1}{2},j}^k{u}_{i-\frac{1}{2},j}^k}{\mathsf{\Delta }x}}-{\displaystyle \frac{H_{i,j+\frac{1}{2}}^k{v}_{i,j+\frac{1}{2}}^k-H_{i,j-\frac{1}{2}}^k{v}_{i,j-\frac{1}{2}}^k}{\mathsf{\Delta }y}}-\\ -{\displaystyle \frac{H_{i+\frac{1}{2},j}^{k+1}{u}_{i+\frac{1}{2},j}^{k+1}-H_{i-\frac{1}{2},j}^{k+1}{u}_{i-\frac{1}{2},j}^{k+1}}{\mathsf{\Delta }x}}+{\displaystyle \frac{H_{i,j+\frac{1}{2}}^{k+1}{v}_{i,j+\frac{1}{2}}^{k+1}-H_{i,j-\frac{1}{2}}^{k+1}{v}_{i,j-\frac{1}{2}}^{k+1}}{\mathsf{\Delta }y}}\end{array}\right).$Taking into account expressions of the form (3), the first equation (1) can be rewritten as:

$${\displaystyle \frac{H_{i,j}^{k+1}-H_{i,j}^k}{\mathsf{\Delta }t}}={\displaystyle \frac{1}{4}}\left(\begin{array}{c}-{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(H_{i+1,j}^k{u}_{i+1,j}^k-H_{i-1,j}^k{u}_{i-1,j}^k\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(H_{i,j+1}^k{v}_{i,j+1}^k-H_{i,j-1}^k{v}_{i,j-1}^k\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(H_{i+1,j}^{k+1}{u}_{i+1,j}^{k+1}-H_{i-1,j}^{k+1}{u}_{i-1,j}^{k+1}\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(H_{i,j+1}^{k+1}{v}_{i,j+1}^{k+1}-H_{i,j-1}^{k+1}{v}_{i,j-1}^{k+1}\right)\end{array}\right).$$

(13)

In expression (13), the terms with superscripts (k+1) have been deliberately left without using a formula of the form (3). This was done deliberately, as it is simpler from a mathematical and computational point of view to linearise these terms first and then approximate the elements with non-integer indices, rather than the reverse – approximating the elements using formulas of the form (3) and then linearising them. In the latter case, there are two main drawbacks: the accuracy of the approximation will be lower, and the number of terms requiring linearisation after expanding the brackets will be significantly greater.

Similarly, a discrete analogue of the second equation can be constructed for a system of equations of type (1). It takes the form:

${\displaystyle \frac{H_{i,j}^{k+1}{u}_{i,j}^{k+1}-H_{i,j}^k{u}_{i,j}^k}{\mathsf{\Delta }t}}=$$={\displaystyle \frac{1}{2}}\left(\begin{array}{c}-{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(\begin{array}{c}{H}_{i+{\displaystyle \frac{1}{2}},j}^k{\left(u_{i+{\displaystyle \frac{1}{2}},j}^k\right)}^2+{\displaystyle \frac{g}{2}}{\left(H_{i+{\displaystyle \frac{1}{2}},j}^k\right)}^2-\\ -H_{i-{\displaystyle \frac{1}{2}},j}^k{\left(u_{i-{\displaystyle \frac{1}{2}},j}^k\right)}^2-{\displaystyle \frac{g}{2}}{\left(H_{i-{\displaystyle \frac{1}{2}},j}^k\right)}^2\end{array}\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(\begin{array}{c}{H}_{i,j+{\displaystyle \frac{1}{2}}}^k{u}_{i,j+{\displaystyle \frac{1}{2}}}^k{v}_{i,j+{\displaystyle \frac{1}{2}}}^k-\\ -H_{i,j-{\displaystyle \frac{1}{2}}}^k{u}_{i,j-{\displaystyle \frac{1}{2}}}^k{v}_{i,j-{\displaystyle \frac{1}{2}}}^k\end{array}\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(\begin{array}{c}{H}_{i+{\displaystyle \frac{1}{2}},j}^{k+1}{\left(u_{i+{\displaystyle \frac{1}{2}},j}^{k+1}\right)}^2+{\displaystyle \frac{g}{2}}{\left(H_{i+{\displaystyle \frac{1}{2}},j}^{k+1}\right)}^2-\\ -H_{i-{\displaystyle \frac{1}{2}},j}^{k+1}{\left(u_{i-{\displaystyle \frac{1}{2}},j}^{k+1}\right)}^2-{\displaystyle \frac{g}{2}}{\left(H_{i-{\displaystyle \frac{1}{2}},j}^{k+1}\right)}^2\end{array}\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(\begin{array}{c}{H}_{i,j+{\displaystyle \frac{1}{2}}}^{k+1}{u}_{i,j+{\displaystyle \frac{1}{2}}}^{k+1}{v}_{i,j+{\displaystyle \frac{1}{2}}}^{k+1}-\\ -H_{i,j-{\displaystyle \frac{1}{2}}}^{k+1}{u}_{i,j-{\displaystyle \frac{1}{2}}}^{k+1}{v}_{i,j-{\displaystyle \frac{1}{2}}}^{k+1}\end{array}\right)\end{array}\right).$If we substitute the formula of type (3) into the latter expression, we obtain:

${\displaystyle \frac{H_{i,j}^{k+1}{u}_{i,j}^{k+1}-H_{i,j}^k{u}_{i,j}^k}{\mathsf{\Delta }t}}=$

$$={\displaystyle \frac{1}{4}}\left(\begin{array}{c}-{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(\begin{array}{c}\left(H_{i+1,j}^k{\left(u_{i+1,j}^k\right)}^2-H_{i-1,j}^k{\left(u_{i-1,j}^k\right)}^2\right)+\\ +{\displaystyle \frac{g\left({\left(H_{i+1,j}^k\right)}^2-{\left(H_{i-1,j}^k\right)}^2\right)}{2}}\end{array}\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(\begin{array}{c}{H}_{i,j+1}^k{u}_{i,j+1}^k{v}_{i,j+1}^k-\\ -H_{i,j-1}^k{u}_{i,j-1}^k{v}_{i,j-1}^k\end{array}\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(\begin{array}{c}\left(H_{i+1,j}^{k+1}{\left(u_{i+1,j}^{k+1}\right)}^2-H_{i-1,j}^{k+1}{\left(u_{i-1,j}^{k+1}\right)}^2\right)+\\ +{\displaystyle \frac{g\left({\left(H_{i+1,j}^{k+1}\right)}^2-{\left(H_{i-1,j}^{k+1}\right)}^2\right)}{2}}\end{array}\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(\begin{array}{c}{H}_{i,j+1}^{k+1}{u}_{i,j+1}^{k+1}{v}_{i,j+1}^{k+1}-\\ -H_{i,j-1}^{k+1}{u}_{i,j-1}^{k+1}{v}_{i,j-1}^{k+1}\end{array}\right)\end{array}\right).$$

(14)

The same principle is used to discretise the third equation of system (1). We have:

${\displaystyle \frac{H_{i,j}^{k+1}{v}_{i,j}^{k+1}-H_{i,j}^k{v}_{i,j}^k}{\mathsf{\Delta }t}}=$$={\displaystyle \frac{1}{2}}\left(\begin{array}{c}-{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(\begin{array}{c}{H}_{i+{\displaystyle \frac{1}{2}},j}^k{u}_{i+{\displaystyle \frac{1}{2}},j}^k{v}_{i+{\displaystyle \frac{1}{2}},j}^k-\\ -H_{i-{\displaystyle \frac{1}{2}},j}^k{u}_{i-{\displaystyle \frac{1}{2}},j}^k{v}_{i-{\displaystyle \frac{1}{2}},j}^k\end{array}\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(\begin{array}{c}{H}_{i,j+{\displaystyle \frac{1}{2}}}^k{\left(v_{i,j+{\displaystyle \frac{1}{2}}}^k\right)}^2+{\displaystyle \frac{g}{2}}{\left(H_{i,j+{\displaystyle \frac{1}{2}}}^k\right)}^2-\\ -H_{i,j-{\displaystyle \frac{1}{2}}}^k{\left(v_{i,j-{\displaystyle \frac{1}{2}}}^k\right)}^2-{\displaystyle \frac{g}{2}}{\left(H_{i,j-{\displaystyle \frac{1}{2}}}^k\right)}^2\end{array}\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(\begin{array}{c}{H}_{i+{\displaystyle \frac{1}{2}},j}^{k+1}{u}_{i+{\displaystyle \frac{1}{2}},j}^{k+1}{v}_{i+{\displaystyle \frac{1}{2}},j}^{k+1}-\\ -H_{i-{\displaystyle \frac{1}{2}},j}^{k+1}{u}_{i-{\displaystyle \frac{1}{2}},j}^{k+1}{v}_{i-{\displaystyle \frac{1}{2}},j}^{k+1}\end{array}\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(\begin{array}{c}{H}_{i,j+{\displaystyle \frac{1}{2}}}^{k+1}{\left(v_{i,j+{\displaystyle \frac{1}{2}}}^{k+1}\right)}^2+{\displaystyle \frac{g}{2}}{\left(H_{i,j+{\displaystyle \frac{1}{2}}}^{k+1}\right)}^2-\\ -H_{i,j-{\displaystyle \frac{1}{2}}}^{k+1}{\left(v_{i,j-{\displaystyle \frac{1}{2}}}^{k+1}\right)}^2-{\displaystyle \frac{g}{2}}{\left(H_{i,j-{\displaystyle \frac{1}{2}}}^{k+1}\right)}^2\end{array}\right)\end{array}\right).$The final difference equation for the third shallow water equation after the transformations will be:

${\displaystyle \frac{H_{i,j}^{k+1}{v}_{i,j}^{k+1}-H_{i,j}^k{v}_{i,j}^k}{\mathsf{\Delta }t}}=$

$$={\displaystyle \frac{1}{4}}\left(\begin{array}{c}-{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(\begin{array}{c}{H}_{i+1,j}^k{u}_{i+1,j}^k{v}_{i+1,j}^k-\\ -H_{i-1,j}^k{u}_{i-1,j}^k{v}_{i-1,j}^k\end{array}\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(\begin{array}{c}{H}_{i,j+1}^k{\left(v_{i,j+1}^k\right)}^2-H_{i,j-1}^k{\left(v_{i,j-1}^k\right)}^2+\\ +{\displaystyle \frac{g\left({\left(H_{i,j+1}^k\right)}^2-{\left(H_{i,j-1}^k\right)}^2\right)}{2}}-\end{array}\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }x}}\left(\begin{array}{c}{H}_{i+1,j}^{k+1}{u}_{i+1,j}^{k+1}{v}_{i+1,j}^{k+1}-\\ -H_{i-1,j}^{k+1}{u}_{i-1,j}^{k+1}{v}_{i-1,j}^{k+1}\end{array}\right)-\\ -{\displaystyle \frac{1}{\mathsf{\Delta }y}}\left(\begin{array}{c}{H}_{i,j+1}^{k+1}{\left(v_{i,j+1}^{k+1}\right)}^2-H_{i,j-1}^{k+1}{\left(v_{i,j-1}^{k+1}\right)}^2+\\ +{\displaystyle \frac{g\left({\left(H_{i,j+1}^{k+1}\right)}^2-{\left(H_{i,j-1}^{k+1}\right)}^2\right)}{2}}-\end{array}\right)\end{array}\right).$$

(15)

In expression (15), just as in expression (13), the terms with superscripts (k+1) have been deliberately left without using a formula of type (3).

The next step is to linearise the non-linear terms in the circuit. We have the following non-linear functions of the form:$F\left(H\right)=H^2,G_1\left(H,u\right)=Hu,G_2\left(H,v\right)=Hv,G_3\left(H,u\right)=H{u}^2,G_4\left(H,v\right)=H{v}^2,R\left(H,u,v\right)=Huv.$ First, we shall implement the non-linear function with one unknown$F\left(H\right)=H^2.$ Based on Newton’s linearisation, we obtain:

${\left[H^2\right]}^{\left(n+1\right)}\approx {\left[H^2\right]}^{\left(n\right)}+{\left[{\left(H^2\right)}^{\text{'}}\right]}^{\left(n\right)}\left({\left[H\right]}^{\left(n+1\right)}-{\left[H\right]}^{\left(n\right)}\right);$${\left[H^2\right]}^{\left(n+1\right)}\approx {\left[H^2\right]}^{\left(n\right)}+{\left[2H\right]}^{\left(n\right)}\left({\left[H\right]}^{\left(n+1\right)}-{\left[H\right]}^{\left(n\right)}\right).$Therefore:

$${\left[H^2\right]}^{\left(n+1\right)}\approx 2{\left[H\right]}^{\left(n\right)}{\left[H\right]}^{\left(n+1\right)}-{\left[H^2\right]}^{\left(n\right)}.$$

(16)

For the node$\left(i+1,j\right)$ at time step k+1, we have:

${\left[{\left(H_{i+1,j}^{k+1}\right)}^2\right]}^{\left(n+1\right)}\approx 2{\left[H_{i+1,j}^{k+1}\right]}^{\left(n\right)}{\left[H_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}-{\left[{\left(H_{i+1,j}^{k+1}\right)}^2\right]}^{\left(n\right)}.$For the node$\left(i-1,j\right)$ at time step k+1, we have:

${\left[{\left(H_{i-1,j}^{k+1}\right)}^2\right]}^{\left(n+1\right)}\approx 2{\left[H_{i-1,j}^{k+1}\right]}^{\left(n\right)}{\left[H_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}-{\left[{\left(H_{i-1,j}^{k+1}\right)}^2\right]}^{\left(n\right)}.$For the node$\left(i,j+1\right)$ at time step k+1, we have:

${\left[{\left(H_{i,j+1}^{k+1}\right)}^2\right]}^{\left(n+1\right)}\approx 2{\left[H_{i,j+1}^{k+1}\right]}^{\left(n\right)}{\left[H_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}-{\left[{\left(H_{i,j+1}^{k+1}\right)}^2\right]}^{\left(n\right)}.$For the node$\left(i,j-1\right)$ at time step k+1, we have:

${\left[{\left(H_{i,j-1}^{k+1}\right)}^2\right]}^{\left(n+1\right)}\approx 2{\left[H_{i,j-1}^{k+1}\right]}^{\left(n\right)}{\left[H_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}-{\left[{\left(H_{i,j-1}^{k+1}\right)}^2\right]}^{\left(n\right)}.$Next, we perform linearisation using Newton’s method for the non-linear function $G_1\left(H,u\right)=Hu$. We obtain:

${\left[Hu\right]}^{\left(n+1\right)}\approx {\left[Hu\right]}^{\left(n\right)}+{\left[{\left(Hu\right)}_H^{\text{'}}\right]}^{\left(n\right)}\left({\left[H\right]}^{\left(n+1\right)}-{\left[H\right]}^{\left(n\right)}\right)+$$+{\left[{\left(Hu\right)}_u^{\text{'}}\right]}^{\left(n\right)}\left({\left[u\right]}^{\left(n+1\right)}-{\left[u\right]}^{\left(n\right)}\right);$${\left[Hu\right]}^{\left(n+1\right)}\approx {\left[Hu\right]}^{\left(n\right)}+{\left[u\right]}^{\left(n\right)}\left({\left[H\right]}^{\left(n+1\right)}-{\left[H\right]}^{\left(n\right)}\right)+{\left[H\right]}^{\left(n\right)}\left({\left[u\right]}^{\left(n+1\right)}-{\left[u\right]}^{\left(n\right)}\right).$Therefore:

$${\left[Hu\right]}^{\left(n+1\right)}\approx {\left[H\right]}^{\left(n+1\right)}{\left[u\right]}^{\left(n\right)}+{\left[u\right]}^{\left(n+1\right)}{\left[H\right]}^{\left(n\right)}-{\left[H\right]}^{\left(n\right)}{\left[u\right]}^{\left(n\right)}.$$

(17)

Using (17) for the node$\left(i+1,j\right)$ at time step k+1, we have:

${\left[H_{i+1,j}^{k+1}{u}_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}\approx {\left[H_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[u_{i+1,j}^{k+1}\right]}^{\left(n\right)}+$$+{\left[u_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i+1,j}^{k+1}\right]}^{\left(n\right)}-{\left[H_{i+1,j}^{k+1}\right]}^{\left(n\right)}{\left[u_{i+1,j}^{k+1}\right]}^{\left(n\right)}.$Using (17) for the node$\left(i-1,j\right)$ at time step k+1, we have:

${\left[H_{i-1,j}^{k+1}{u}_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}\approx {\left[H_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[u_{i-1,j}^{k+1}\right]}^{\left(n\right)}+$$+{\left[u_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i-1,j}^{k+1}\right]}^{\left(n\right)}-{\left[H_{i-1,j}^{k+1}\right]}^{\left(n\right)}{\left[u_{i-1,j}^{k+1}\right]}^{\left(n\right)}.$Using (17) for the node$\left(i,j+1\right)$ at time step k+1, we have:

${\left[H_{i,j+1}^{k+1}{u}_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}\approx {\left[H_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}{\left[u_{i,j+1}^{k+1}\right]}^{\left(n\right)}+$$+{\left[u_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i,j+1}^{k+1}\right]}^{\left(n\right)}-{\left[H_{i,j+1}^{k+1}\right]}^{\left(n\right)}{\left[u_{i,j+1}^{k+1}\right]}^{\left(n\right)}.$Using (17) for the node$\left(i,j-1\right)$ at time step k+1, we have:

${\left[H_{i,j-1}^{k+1}{u}_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}\approx {\left[H_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}{\left[u_{i,j-1}^{k+1}\right]}^{\left(n\right)}+$$+{\left[u_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i,j-1}^{k+1}\right]}^{\left(n\right)}-{\left[H_{i,j-1}^{k+1}\right]}^{\left(n\right)}{\left[u_{i,j-1}^{k+1}\right]}^{\left(n\right)}.$Newton’s linearisation of the function$G_2\left(H,v\right)=Hv$ takes the following form:

$${\left[Hv\right]}^{\left(n+1\right)}\approx {\left[H\right]}^{\left(n+1\right)}{\left[v\right]}^{\left(n\right)}+{\left[v\right]}^{\left(n+1\right)}{\left[H\right]}^{\left(n\right)}-{\left[H\right]}^{\left(n\right)}{\left[v\right]}^{\left(n\right)}.$$

(18)

Then using (18) for the node$\left(i+1,j\right)$ at time step k+1:

${\left[H_{i+1,j}^{k+1}{v}_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}\approx {\left[H_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[v_{i+1,j}^{k+1}\right]}^{\left(n\right)}+$$+{\left[v_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i+1,j}^{k+1}\right]}^{\left(n\right)}-{\left[H_{i+1,j}^{k+1}\right]}^{\left(n\right)}{\left[v_{i+1,j}^{k+1}\right]}^{\left(n\right)};$for the node $\left(i-1,j\right)$ at time step k+1:

${\left[H_{i-1,j}^{k+1}{v}_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}\approx {\left[H_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[v_{i-1,j}^{k+1}\right]}^{\left(n\right)}+$$+{\left[v_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i-1,j}^{k+1}\right]}^{\left(n\right)}-{\left[H_{i-1,j}^{k+1}\right]}^{\left(n\right)}{\left[v_{i-1,j}^{k+1}\right]}^{\left(n\right)};$for the node$\left(i,j+1\right)$ at time step k+1:

${\left[H_{i,j+1}^{k+1}{v}_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}\approx {\left[H_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}{\left[v_{i,j+1}^{k+1}\right]}^{\left(n\right)}+$$+{\left[v_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i,j+1}^{k+1}\right]}^{\left(n\right)}-{\left[H_{i,j+1}^{k+1}\right]}^{\left(n\right)}{\left[v_{i,j+1}^{k+1}\right]}^{\left(n\right)};$for the node$\left(i,j-1\right)$ at time step k+1:

${\left[H_{i,j-1}^{k+1}{v}_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}\approx {\left[H_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}{\left[v_{i,j-1}^{k+1}\right]}^{\left(n\right)}+$$+{\left[v_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i,j-1}^{k+1}\right]}^{\left(n\right)}-{\left[H_{i,j-1}^{k+1}\right]}^{\left(n\right)}{\left[v_{i,j-1}^{k+1}\right]}^{\left(n\right)}.$Now let us perform linearisation using Newton’s method for the non-linear function $G_3\left(H,u\right)=H{u}^2$. We obtain:

${\left[H{u}^2\right]}^{\left(n+1\right)}\approx {\left[H{u}^2\right]}^{\left(n\right)}+{\left[{\left(H{u}^2\right)}_H^{\text{'}}\right]}^{\left(n\right)}\left({\left[H\right]}^{\left(n+1\right)}-{\left[H\right]}^{\left(n\right)}\right)+$$+{\left[{\left(H{u}^2\right)}_u^{\text{'}}\right]}^{\left(n\right)}\left({\left[u\right]}^{\left(n+1\right)}-{\left[u\right]}^{\left(n\right)}\right),$${\left[H{u}^2\right]}^{\left(n+1\right)}\approx {\left[H{u}^2\right]}^{\left(n\right)}+{\left[u^2\right]}^{\left(n\right)}\left({\left[H\right]}^{\left(n+1\right)}-{\left[H\right]}^{\left(n\right)}\right)+$$+2{\left[H\right]}^{\left(n\right)}{\left[u\right]}^{\left(n\right)}\left({\left[u\right]}^{\left(n+1\right)}-{\left[u\right]}^{\left(n\right)}\right).$Therefore:

$${\left[H{u}^2\right]}^{\left(n+1\right)}\approx {\left[H\right]}^{\left(n+1\right)}{\left[u^2\right]}^{\left(n\right)}+2{\left[u\right]}^{\left(n+1\right)}{\left[H\right]}^{\left(n\right)}{\left[u\right]}^{\left(n\right)}-2{\left[H\right]}^{\left(n\right)}{\left[u^2\right]}^{\left(n\right)}.$$

(19)

Then using (19) for the node $\left(i+1,j\right)$ at time step k+1:

${\left[H_{i+1,j}^{k+1}{\left(u_{i+1,j}^{k+1}\right)}^2\right]}^{\left(n+1\right)}\approx {\left[H_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[{\left(u_{i+1,j}^{k+1}\right)}^2\right]}^{\left(n\right)}+$$+2{\left[u_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i+1,j}^{k+1}\right]}^{\left(n\right)}{\left[u_{i+1,j}^{k+1}\right]}^{\left(n\right)}-2{\left[H_{i+1,j}^{k+1}\right]}^{\left(n\right)}{\left[{\left(u_{i+1,j}^{k+1}\right)}^2\right]}^{\left(n\right)};$for the node$\left(i-1,j\right)$ at time step k+1:

${\left[H_{i-1,j}^{k+1}{\left(u_{i-1,j}^{k+1}\right)}^2\right]}^{\left(n+1\right)}\approx {\left[H_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[{\left(u_{i-1,j}^{k+1}\right)}^2\right]}^{\left(n\right)}+$$+2{\left[u_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i-1,j}^{k+1}\right]}^{\left(n\right)}{\left[u_{i-1,j}^{k+1}\right]}^{\left(n\right)}-2{\left[H_{i-1,j}^{k+1}\right]}^{\left(n\right)}{\left[{\left(u_{i-1,j}^{k+1}\right)}^2\right]}^{\left(n\right)};$for the node$\left(i,j+1\right)$ at time step k+1:

${\left[H_{i,j+1}^{k+1}{\left(u_{i,j+1}^{k+1}\right)}^2\right]}^{\left(n+1\right)}\approx {\left[H_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}{\left[{\left(u_{i,j+1}^{k+1}\right)}^2\right]}^{\left(n\right)}+$$+2{\left[u_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i,j+1}^{k+1}\right]}^{\left(n\right)}{\left[u_{i,j+1}^{k+1}\right]}^{\left(n\right)}-2{\left[H_{i,j+1}^{k+1}\right]}^{\left(n\right)}{\left[{\left(u_{i,j+1}^{k+1}\right)}^2\right]}^{\left(n\right)};$for the node$\left(i,j-1\right)$ at time step k+1:

${\left[H_{i,j-1}^{k+1}{\left(u_{i,j-1}^{k+1}\right)}^2\right]}^{\left(n+1\right)}\approx {\left[H_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}{\left[{\left(u_{i,j-1}^{k+1}\right)}^2\right]}^{\left(n\right)}+$$+2{\left[u_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i,j-1}^{k+1}\right]}^{\left(n\right)}{\left[u_{i,j-1}^{k+1}\right]}^{\left(n\right)}-2{\left[H_{i,j-1}^{k+1}\right]}^{\left(n\right)}{\left[{\left(u_{i,j-1}^{k+1}\right)}^2\right]}^{\left(n\right)}.$Let us perform linearisation using Newton’s method for the non-linear function $G_3\left(H,u\right)=H{v}^2$. We obtain:

$${\left[H{v}^2\right]}^{\left(n+1\right)}\approx {\left[H\right]}^{\left(n+1\right)}{\left[v^2\right]}^{\left(n\right)}+2{\left[v\right]}^{\left(n+1\right)}{\left[H\right]}^{\left(n\right)}{\left[v\right]}^{\left(n\right)}-2{\left[H\right]}^{\left(n\right)}{\left[v^2\right]}^{\left(n\right)}.$$

(20)

Then using (20) for the node $\left(i+1,j\right)$ at time step k+1:

${\left[H_{i+1,j}^{k+1}{\left(v_{i+1,j}^{k+1}\right)}^2\right]}^{\left(n+1\right)}\approx {\left[H_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[{\left(v_{i+1,j}^{k+1}\right)}^2\right]}^{\left(n\right)}+$$+2{\left[v_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i+1,j}^{k+1}\right]}^{\left(n\right)}{\left[v_{i+1,j}^{k+1}\right]}^{\left(n\right)}-2{\left[H_{i+1,j}^{k+1}\right]}^{\left(n\right)}{\left[{\left(v_{i+1,j}^{k+1}\right)}^2\right]}^{\left(n\right)};$for the node$\left(i-1,j\right)$ at time step k+1:

${\left[H_{i-1,j}^{k+1}{\left(v_{i-1,j}^{k+1}\right)}^2\right]}^{\left(n+1\right)}\approx {\left[H_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[{\left(v_{i-1,j}^{k+1}\right)}^2\right]}^{\left(n\right)}+$$+2{\left[v_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i-1,j}^{k+1}\right]}^{\left(n\right)}{\left[v_{i-1,j}^{k+1}\right]}^{\left(n\right)}-2{\left[H_{i-1,j}^{k+1}\right]}^{\left(n\right)}{\left[{\left(v_{i-1,j}^{k+1}\right)}^2\right]}^{\left(n\right)};$for the node $\left(i,j+1\right)$ at time step k+1:

${\left[H_{i,j+1}^{k+1}{\left(v_{i,j+1}^{k+1}\right)}^2\right]}^{\left(n+1\right)}\approx {\left[H_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}{\left[{\left(v_{i,j+1}^{k+1}\right)}^2\right]}^{\left(n\right)}+$$+2{\left[v_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i,j+1}^{k+1}\right]}^{\left(n\right)}{\left[v_{i,j+1}^{k+1}\right]}^{\left(n\right)}-2{\left[H_{i,j+1}^{k+1}\right]}^{\left(n\right)}{\left[{\left(v_{i,j+1}^{k+1}\right)}^2\right]}^{\left(n\right)};$for the node $\left(i,j-1\right)$ at time step k+1:

${\left[H_{i,j-1}^{k+1}{\left(v_{i,j-1}^{k+1}\right)}^2\right]}^{\left(n+1\right)}\approx {\left[H_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}{\left[{\left(v_{i,j-1}^{k+1}\right)}^2\right]}^{\left(n\right)}+$$+2{\left[v_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i,j-1}^{k+1}\right]}^{\left(n\right)}{\left[v_{i,j-1}^{k+1}\right]}^{\left(n\right)}-2{\left[H_{i,j-1}^{k+1}\right]}^{\left(n\right)}{\left[{\left(v_{i,j-1}^{k+1}\right)}^2\right]}^{\left(n\right)}.$Finally, we perform linearisation using Newton’s method for the non-linear function $G_4\left(H,u,v\right)=Huv$. We obtain:

${\left[Huv\right]}^{\left(n+1\right)}\approx {\left[Huv\right]}^{\left(n\right)}+{\left[{\left(Huv\right)}_H^{\text{'}}\right]}^{\left(n\right)}\left({\left[H\right]}^{\left(n+1\right)}-{\left[H\right]}^{\left(n\right)}\right)+$$+{\left[{\left(Huv\right)}_u^{\text{'}}\right]}^{\left(n\right)}\left({\left[u\right]}^{\left(n+1\right)}-{\left[u\right]}^{\left(n\right)}\right)+{\left[{\left(Huv\right)}_v^{\text{'}}\right]}^{\left(n\right)}\left({\left[v\right]}^{\left(n+1\right)}-{\left[v\right]}^{\left(n\right)}\right),$${\left[Huv\right]}^{\left(n+1\right)}\approx {\left[Huv\right]}^{\left(n\right)}+{\left[uv\right]}^{\left(n\right)}\left({\left[H\right]}^{\left(n+1\right)}-{\left[H\right]}^{\left(n\right)}\right)+$$+{\left[Hv\right]}^{\left(n\right)}\left({\left[u\right]}^{\left(n+1\right)}-{\left[u\right]}^{\left(n\right)}\right)+{\left[Hu\right]}^{\left(n\right)}\left({\left[v\right]}^{\left(n+1\right)}-{\left[v\right]}^{\left(n\right)}\right).$Therefore:

$${\left[Huv\right]}^{\left(n+1\right)}\approx {\left[H\right]}^{\left(n+1\right)}{\left[u\right]}^{\left(n\right)}{\left[v\right]}^{\left(n\right)}+{\left[u\right]}^{\left(n+1\right)}{\left[H\right]}^{\left(n\right)}{\left[v\right]}^{\left(n\right)}+$$

(21)

$+{\left[v\right]}^{\left(n+1\right)}{\left[H\right]}^{\left(n\right)}{\left[u\right]}^{\left(n\right)}-2{\left[H\right]}^{\left(n\right)}{\left[u\right]}^{\left(n\right)}{\left[v\right]}^{\left(n\right)}.$Then using (21) for the node$\left(i+1,j\right)$ at time step k+1:

${\left[H_{i+1,j}^{k+1}{u}_{i+1,j}^{k+1}{v}_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}\approx {\left[H_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[u_{i+1,j}^{k+1}\right]}^{\left(n\right)}{\left[v_{i+1,j}^{k+1}\right]}^{\left(n\right)}+$$+{\left[u_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i+1,j}^{k+1}\right]}^{\left(n\right)}{\left[v_{i+1,j}^{k+1}\right]}^{\left(n\right)}+{\left[v_{i+1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i+1,j}^{k+1}\right]}^{\left(n\right)}{\left[u_{i+1,j}^{k+1}\right]}^{\left(n\right)}-$$-2{\left[H_{i+1,j}^{k+1}\right]}^{\left(n\right)}{\left[u_{i+1,j}^{k+1}\right]}^{\left(n\right)}{\left[v_{i+1,j}^{k+1}\right]}^{\left(n\right)};$for the node$\left(i-1,j\right)$ at time step k+1:

${\left[H_{i-1,j}^{k+1}{u}_{i-1,j}^{k+1}{v}_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}\approx {\left[H_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[u_{i-1,j}^{k+1}\right]}^{\left(n\right)}{\left[v_{i-1,j}^{k+1}\right]}^{\left(n\right)}+$$+{\left[u_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i-1,j}^{k+1}\right]}^{\left(n\right)}{\left[v_{i-1,j}^{k+1}\right]}^{\left(n\right)}+{\left[v_{i-1,j}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i-1,j}^{k+1}\right]}^{\left(n\right)}{\left[u_{i-1,j}^{k+1}\right]}^{\left(n\right)}-$$-2{\left[H_{i-1,j}^{k+1}\right]}^{\left(n\right)}{\left[u_{i-1,j}^{k+1}\right]}^{\left(n\right)}{\left[v_{i-1,j}^{k+1}\right]}^{\left(n\right)};$for the node$\left(i,j+1\right)$ at time step k+1:

${\left[H_{i,j+1}^{k+1}{u}_{i,j+1}^{k+1}{v}_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}\approx {\left[H_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}{\left[u_{i,j+1}^{k+1}\right]}^{\left(n\right)}{\left[v_{i,j+1}^{k+1}\right]}^{\left(n\right)}+$$+{\left[u_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i,j+1}^{k+1}\right]}^{\left(n\right)}{\left[v_{i,j+1}^{k+1}\right]}^{\left(n\right)}+{\left[v_{i,j+1}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i,j+1}^{k+1}\right]}^{\left(n\right)}{\left[u_{i,j+1}^{k+1}\right]}^{\left(n\right)}-$$-2{\left[H_{i,j+1}^{k+1}\right]}^{\left(n\right)}{\left[u_{i,j+1}^{k+1}\right]}^{\left(n\right)}{\left[v_{i,j+1}^{k+1}\right]}^{\left(n\right)};$for the node$\left(i,j-1\right)$ at time step k+1:

${\left[H_{i,j-1}^{k+1}{u}_{i,j-1}^{k+1}{v}_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}\approx {\left[H_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}{\left[u_{i,j-1}^{k+1}\right]}^{\left(n\right)}{\left[v_{i,j-1}^{k+1}\right]}^{\left(n\right)}+$$+{\left[u_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i,j-1}^{k+1}\right]}^{\left(n\right)}{\left[v_{i,j-1}^{k+1}\right]}^{\left(n\right)}+{\left[v_{i,j-1}^{k+1}\right]}^{\left(n+1\right)}{\left[H_{i,j-1}^{k+1}\right]}^{\left(n\right)}{\left[u_{i,j-1}^{k+1}\right]}^{\left(n\right)}-$$-2{\left[H_{i,j-1}^{k+1}\right]}^{\left(n\right)}{\left[u_{i,j-1}^{k+1}\right]}^{\left(n\right)}{\left[v_{i,j-1}^{k+1}\right]}^{\left(n\right)}.$Finally, we obtain an iterative implicit scheme for (1), consisting of discrete analogues of equations (13)–(15) and an iterative linearisation procedure (16)–(21).

4. Developed Simulation Software

The developed software module is intended for numerical modelling of two-dimensional hydrodynamic processes in shallow water bodies, based on the nonlinear shallow-water equations. The module is implemented in the MATLAB environment as a single executable file, which provides a full cycle of calculations - from initialization of model parameters to saving and visualization of results. It is based on a fully implicit finite-difference scheme on a uniform rectangular grid, which guarantees unconditional stability of the calculations and allows the use of significant time-integration steps, regardless of Courant restrictions. Spatial discretization is performed with second-order accuracy, and integration over time is performed using the implicit Euler scheme of the first order.

To overcome the nonlinearity of the original system at each time step, Newton’s iterative method is used, which provides quadratic convergence in the vicinity of the solution. In the process of linearization, the partial derivatives of all nonlinear flow terms are analytically calculated, thereby avoiding errors from numerical differentiation. At each iteration, a system of linear algebraic equations with a sparse strip matrix is formed, which is diagonally dominant due to the locality of finite-difference patterns.м The solution of this system is performed by the direct method using the built-in MATLAB tools for working with sparse matrices, which ensures high reliability and accuracy of calculations. The maximum absolute error criterion across all grid nodes controls the convergence of the iterations.

Flexible settings characterise the software module: the user can specify the dimensions of the computational domain, grid spacing, time integration step, relaxation parameters, and initial conditions as a localised perturbation of the free surface. Initial velocities are taken equal to zero, which corresponds to a state of rest until the moment of perturbation. No-flow conditions are applied at the domain boundaries to simulate the reflection of waves from the impermeable walls of a closed reservoir. The modelling results are presented as two-dimensional and three-dimensional colour maps of the free-surface height distribution, as well as an automatically generated animation that allows you to track the evolution of the wave field over time. The data is stored in a MAT file for further analysis.

A feature of the module is its orientation towards application in information and measurement technologies for monitoring water bodies. It can be integrated with additional blocks that implement the impurity transport equation, adaptive mesh thickening, parallel calculations and data assimilation. The use of an implicit scheme in combination with an efficient sparse systems solver allows you to perform calculations on large-dimensional grids with acceptable time costs, which opens up prospects for using the module in real-time systems, in particular for solving inverse problems of identifying pollution sources by jointly using direct hydrodynamic calculations and optimization methods.

Source [25] presents a project containing the program code for a nonlinear numerical process modelling scheme based on the shallow water equations.

5. Numerical Experiments

Let us consider a problem with the following mathematical formulation:

$$\left\{\begin{array}{c}{\displaystyle \frac{\partial H}{\partial t}}+{\displaystyle \frac{\partial \left(Hu\right)}{\partial x}}+{\displaystyle \frac{\partial \left(Hv\right)}{\partial y}}=0;\\ {\displaystyle \frac{\partial \left(Hu\right)}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left(H{u}^2+{\displaystyle \frac{g}{2}}{H}^2\right)+{\displaystyle \frac{\partial \left(Huv\right)}{\partial y}}=0;\\ {\displaystyle \frac{\partial \left(Hv\right)}{\partial t}}+{\displaystyle \frac{\partial \left(Huv\right)}{\partial x}}+{\displaystyle \frac{\partial }{\partial y}}\left(H{v}^2+{\displaystyle \frac{g}{2}}{H}^2\right)=0.\end{array}\right.$$

(22)

Computational domain: $\mathsf{\Omega }:\left(x;y\right)\in {\left[0;100\right]}^2,t\in \left[0;25\right].$Initial conditions for the model:

$H\left(x,y,0\right)=2-min\left({\displaystyle \frac{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2}{100}},1\right),x_0=27,y_0=23;$

$$u\left(x,y,0\right)=0,v\left(x,y,0\right)=0.$$

(23)

The boundary conditions in the model are as follows:

$\left.{\displaystyle \frac{\partial H}{\partial \overline{n}}}\right|=0,\forall t>0,$

$$Hu\left(0,y,t\right)=0,Hu\left(100,y,t\right)=0;$$

(24)

$Hu\left(0,y,t\right)=0,Hu\left(100,y,t\right)=0;$$Hv\left(x,0,t\right)=0,Hv\left(x,100,t\right)=0.$Let us recall the notation in the model (22)–(24):

• $H\left(x,y,t\right)$ – total water depth;
• $u\left(x,y,t\right)$ – velocity along$x$;
• $v\left(x,y,t\right)$ – velocity along the$y$;
• g – acceleration due to gravity (gravitational acceleration).

To find the solution to problem (22)–(24) using an explicit scheme, a step size of $\mathsf{\Delta }t\le 0.01$ (the Courant–Friedrichs–Lewy condition). For the implicit difference scheme, the step size $\mathsf{\Delta }t=0.1,$ $\mathsf{\Delta }t=0.5,$ $\mathsf{\Delta }t=0.9$ was chosen. The results of modelling the problem using the implicit difference scheme are shown in Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7.

The developed application software as a module easily integrates with real data to detect flooding and other complex hydrodynamic processes.

6. Conclusions

This paper presents a theoretical generalisation and a practical implementation of numerical modelling of hydrodynamic processes based on the two-dimensional shallow-water equations. In the context of the modern paradigm of integrating mathematical models with geographic information systems, this study confirms that the use of shallow-water equations in combination with spatial data (digital elevation models, bed surface parameters) enables the creation of a powerful toolkit for flood forecasting, risk assessment, and water resource management.

Computational experiments on a test problem with an initial local rise in water level (simulating a breach or wave pulse) confirmed the effectiveness of the developed implicit scheme. The resulting graphical dependencies (Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7) illustrate the correct propagation of excitation waves across the domain, accounting for boundary reflections, consistent with the physical meaning of the problem and the specified boundary conditions.

The modified implicit difference scheme proposed in this work is an effective tool for constructing hydrodynamic models intended for further use in a geoinformation environment. Its application enables improved stability and efficiency in calculations when modelling complex water systems, an important step in developing decision-support systems for emergency prevention and the rational use of water resources. Prospects for further research lie in implementing this scheme on grids and directly integrating the software code with existing GIS platforms.