**6. Iterative Algorithm**

As it was shown in the Sections 3 and 5, the problem to find Φ(*i*) = (*U*- (*i*), *ξ*(*i*))*<sup>T</sup>*, *i* = 1, 2 and the additional boundary functions *v*, *ds* is uniquely and densely solvable. Therefore the functions <sup>Φ</sup>(*i*)(*α*), *<sup>v</sup>*(*α*), *ds*(*α*) satisfying Equations (9) and (28)–(30) could be taken as an approximation to the solution of the original problem [26]. An approximation to <sup>Φ</sup>(*i*)(*α*), *<sup>v</sup>*(*α*), *ds*(*α*) could be found with the following iterative algorithm.

1. Let *uk* = (*dks*, *v<sup>k</sup>*)*<sup>T</sup>* be found. Solve the problems in each subdomain Ω*i*, *i* = 1, 2:

$$\begin{cases} \frac{\mathbf{U}^{(i),k}}{\Delta t} + \begin{bmatrix} 0 & -\ell \\ \ell & 0 \end{bmatrix} \mathbf{J}^{(i),k} + R\_f \mathbf{J}^{(i),k} - \mathbf{g} \cdot \mathbf{grad} \,\xi^{(i),k} = \tilde{\boldsymbol{f}}' & \text{in } \Omega\_{i\prime} \\\ \frac{\xi^{(i),k}}{\Delta t} - \mathbf{div} \left( H \Omega^{(i),k} \right) = f\_3, & \text{in } \Omega\_{i\prime} \\\ H \boldsymbol{\Omega}^{(1),k} \cdot \vec{\boldsymbol{n}}\_1 + m\_{op} \sqrt{\mathfrak{g}} \boldsymbol{H}^{(1),k} = m\_{op} \sqrt{\mathfrak{g}} \boldsymbol{H}^k\_s + m\_{in} \sqrt{\mathfrak{g}} \boldsymbol{H}^k, & \text{on } \partial \Omega\_{1\prime} \\\ H \boldsymbol{\Omega}^{(2),k} \cdot \vec{\boldsymbol{n}}\_2 = -m\_{in} \sqrt{\mathfrak{g}} \boldsymbol{H}^k, & \text{on } \partial \Omega\_2. \end{cases} \tag{31}$$

2. Solve the adjoint problems in Ω*i*, *i* = 1, 2:

$$\begin{cases} \frac{\tilde{\mathbf{U}}^{\ast,(i),k}}{\Delta t} - \begin{bmatrix} 0 & -\ell \\ \ell & 0 \end{bmatrix} \dot{\Omega}^{\ast,(i),k} + R\_{f} \dot{\Omega}^{\ast,(i),k} + \mathcal{g} \cdot \mathbf{grad} \, \xi^{\ast,(i),k} = 0, \quad \text{in } \Omega\_{i\prime} \\\frac{\tilde{\mathbf{g}}^{\ast,(i),k}}{\Delta t} + \text{div} \left( H \dot{\Omega}^{\ast,(i),k} \right) = 0, \quad \text{in } \Omega\_{i\prime} \\\ -H \dot{\Omega}^{\ast,(1),k} \cdot \vec{\boldsymbol{n}}\_{1} + m\_{\partial P} \sqrt{\tilde{g} H} \xi^{\ast,(1),k} = m\_{\partial P} \sqrt{\tilde{g} H} (\xi^{(1),k} - \underline{\boldsymbol{\varepsilon}}\_{\text{obs}}) + m\_{\text{in}} \sqrt{\tilde{g} H} (\xi^{(1),k} - \xi^{(2),k}), \quad \text{on } \partial \Omega\_{1\prime}, \Gamma\_{1} \\\ -H \dot{\Omega}^{\ast,(2),k} \cdot \vec{\boldsymbol{n}}\_{2} = m\_{\text{in}} \sqrt{\tilde{g} H} (\xi^{(1),k} - \xi^{(2),k}), \quad \text{on } \partial \Omega\_{2}. \end{cases} \tag{32}$$

3. Find new *uk*+<sup>1</sup> = (*dk*+<sup>1</sup> *s* , *vk*+<sup>1</sup>)*<sup>T</sup>* by:

$$d\_s^{k+1} = d\_s^k - \tau\_k (ad\_s^k + \xi^{\*,(1),k}), \quad \text{ on } \Gamma\_{op\prime} \tag{33}$$

$$
\sigma^{k+1} = \upsilon^k - \tau\_k (a\upsilon^k + (\xi^{\*,(1),k} + \xi^{\*,(2),k})), \quad \text{on } \Gamma\_{in}.\tag{34}
$$

Using assertions of [26] we obtain the validity of the following theorem:

**Theorem 1.** *(I) If ξ*(0) *obs – exact data, ξobs – measured data of observations (possibly containing errors), ξobs* − *ξ*(0) *obs<sup>L</sup>*2(<sup>Γ</sup>*op*) ≤ *δobs, δobs* > 0*, ξ*(1)(*α*) *and ξ*(2)(*α*) *are the solutions of the optimality systems of Equations* (9) *and* (28)*–*(30)*, i* = 1, 2*, then the following assessment is valid:*

$$\left( \left\| \boldsymbol{\xi}^{(1)}(\boldsymbol{a}) - \boldsymbol{\xi}\_{\mathrm{obs}} \right\|\_{L^{\mathbb{M}}\_{2}(\Gamma\_{\mathrm{op}})}^{2} + \left\| \boldsymbol{\xi}^{(1)}(\boldsymbol{a}) - \boldsymbol{\xi}^{(2)}(\boldsymbol{a}) \right\|\_{L^{\mathbb{M}}\_{2}(\Gamma\_{\mathrm{in}})}^{2} \right)^{1/2} \leq c\_{\mathsf{a}} \sqrt{\mathsf{a}} + c\_{\mathsf{b}} \delta\_{\mathrm{obs}} \tag{35}$$

*where cα* = *const* > 0*, cδ* = *const* > 0*.*

*(II) If ϕ* ∈ *<sup>R</sup>*(*A*)*, then the inverse problem of Equation* (22) *has a unique normal solution u*0 = (*ds*,0, *<sup>v</sup>*0)*T. In that case for enough small τk* = *τ* > 0 *vector function <sup>u</sup><sup>k</sup>*(*α*)=(*dks*(*α*), *v<sup>k</sup>*(*α*))*<sup>T</sup> tends to u*0 *in* **H***c, when k* → <sup>∞</sup>*, α* → +0*.*

When the stopping criterion of iterative algorithm is satisfied, the functions *dks*, *<sup>v</sup>k*, *ξ*(*i*),*k*, *U*- (*i*),*k* are taken as an approximate solution of the considered Equations (5), (8) and (9) in Ω*i*, *i* = 1, 2. Suitable stopping criterion should be chosen depending on the iterative parameters of the algorithm, size of the modelled region, numerical method, and a given accuracy. More often it is chosen as a limitation of iterations or of residual value. Equations (31)–(34) converge for the small enough parameter *τk* = *τ* = *const* > 0. However, *τk* may be chosen as follows [26]:

$$\tau\_{k} = \frac{I\_{a}(\boldsymbol{u}\_{a}^{k})}{||\boldsymbol{f}\_{a}^{\prime}(\boldsymbol{u}\_{a}^{k})||^{2}} = \frac{1}{2} \frac{\int\_{\boldsymbol{\tau}\_{op}} \sqrt{\boldsymbol{g}^{\prime}(\boldsymbol{\xi}^{(1)} - \boldsymbol{\xi}\_{obs})^{2}} \, d\boldsymbol{\Gamma} + \int\_{\boldsymbol{\tau}\_{in}} \sqrt{\boldsymbol{g}^{\prime}\boldsymbol{H}} (\boldsymbol{\xi}^{(1)} - \boldsymbol{\xi}^{(2)})^{2} \, d\boldsymbol{\Gamma}}{\int\_{\boldsymbol{\tau}\_{op}} \sqrt{\boldsymbol{g}^{\prime}\boldsymbol{H}} (\boldsymbol{\xi}^{\*,(1)})^{2} \, d\boldsymbol{\Gamma} + \int\_{\boldsymbol{\tau}\_{in}} \sqrt{\boldsymbol{g}^{\prime}\boldsymbol{H}} (\boldsymbol{\xi}^{\*,(1)} + \boldsymbol{\xi}^{\*,(2)})^{2} \, d\boldsymbol{\Gamma}} \,. \tag{36}$$

This choice of the parameter *τk* could be helpful to reduce the number of iterations required.

#### **7. Numerical Experiments and Discussion**

The described approach to domain decomposition and variational DA is applied to the system of shallow water systems of Equation (4). Here Ω is a domain on the plane, and (*<sup>x</sup>*, *y*) are Cartesian coordinates.

To set initial conditions, a preliminary calculation without the domain decomposition and variational DA methods is carried out. In this case the domain is represented by the rectangular [−*L*, *L*] × [0, *L*] with *L* = 100 m. The value of the physical parameters are *g* = 9.81 m/s2, *l* = 0, *Rf* = 0 and *H* = −0.7*x*/*L* + 1 (m). The forces ˜ - *f* and (-*<sup>F</sup>*)3 are equal to 0. The boundary conditions are *H*(*U*- ·*n*) = 0. For this preliminary problem the conditions at *t* = 0 are given by (*p* is a subscript for the preliminary results):

$$\begin{aligned} u\_p &= 0, \ v\_p = 0, \\ \xi\_p &= 10 \exp\left(-\frac{(x-25)^2 + (y-30)^2}{100}\right). \end{aligned}$$

To discretize in time, an implicit scheme is used. The time step is 0.5 s. The finite difference method with the first order space discretization is applied. The spatial grid is uniform, the spatial grid steps are 1 m. The results of simulation after 25 s are taken as initial conditions to the next series of numerical experiments.

To set the observational data on <sup>Γ</sup>*op*, the preliminary experiment is continued on the next 5 s. At each time step after the first 25 s the observational data are chosen as *ξobs* = *ξ p*(0, *y*) · (1 + 0.1*a* − 0.1*b*), where *a* and *b* are random numbers in the range [0, <sup>1</sup>). The numbers are received by pseudo-random generators for the uniform distribution law. Thus, the obtained observational data are artificially noisy (noise level is 0.1).

So, we ge<sup>t</sup> the initial conditions and the observational data for the next experiments with previously described approach of application of DA and DDMs.

For the numerical experiments the domain Ω is the square [0, *L*] × [0, *<sup>L</sup>*]. It is decomposed into the two subdomains Ω1 = [0, *L*/2] × [0, *L*] and Ω2 = [*L*/2, *L*] × [0, *L*] without overlap. We denote the "inner" boundary by Γ*in* = *L*/2 × [0, *L*] and the "outer" boundary by <sup>Γ</sup>*op* = 0 × [0, *<sup>L</sup>*]. Here *L* = 100 m. The value of the physical parameters are *g* = 9.81 m/s2, *l* = 0, *Rf* = 0. The depth of the modelling domain is *H* = −0.7*x*/*L* + 1 (m). The forces ˜ - *f* and (-*<sup>F</sup>*)3 are equal to 0.

To discretize in time, an implicit scheme is chosen. The time step is 0.5 s. The finite difference method with the first order space discretization is applied. The spatial grid is uniform, with grid steps 1 m. The described algorithm of domain decomposition and variational DA is implemented at each time step.

Table 1 gives the residual value depending on the noise level and the stopping criterion (the number of iterations less than 10 or less than 50). The residual value here is defined by

$$\text{Res} = \left( \int\_{\Gamma\_{\mathcal{P}}} \sqrt{\text{g}H} (\mathfrak{J}^{(1)} - \mathfrak{J}\_{\text{obs}})^2 \, d\Gamma + \int\_{\Gamma\_{\text{in}}} \sqrt{\text{g}H} (\mathfrak{J}^{(1)} - \mathfrak{J}^{(2)})^2 \, d\Gamma \right)^{1/2}$$

From Table 1 one can see that the noise level has a significant impact on the residual value at the 50th iteration and has almost no affect at the 10th iteration.

**Table 1.** The residual value (*Res*) depending on the stopping criterion and noise level


Further in this section the stopping criterion is the number of iterations less than 50. Figure 2 shows the results of modelling for the preliminary problem and the results of modelling using domain decomposition and variational DA. The "inner" boundary is presented as a white line. As seen in Figure 2a,b, the results differ from each other. It is connected with the application of variational DA on the open boundary, because the observational data *ξobs* used as closure condition Equation (8) have some additional noise and the results reproduce it (as can be seen in Figure 3).

**Figure 2.** Sea level function (m) at *t* = 5 s after the preliminary simulation: (**a**) Part of the preliminary calculation. (**b**) Simulation with domain decomposition and data assimilation (DA) methods.

**Figure 3.** Sea level function on <sup>Γ</sup>*op* and observational data (1—observational data; 2—results of simulation; 3—results of the preliminary calculation).

The iterative algorithm converges with the chosen parameter *τk* from Equation (36). Figure 4 shows the results of sea level functions at the first, the ninth and the last iterations. For clarity, the observational data are also presented.

**Figure 4.** Sea level function on <sup>Γ</sup>*op* at different iterations and observational data (1—observational data; 2—the first iteration; 3—the last iteration; 4—the ninth iteration).

In order to analyse the impact on the simulation results of the DDM, the comparison of sea level function results at the first and the last iterations on Γ*in* is presented in Figure 5. We could note that *ξ*(1) almost coincides with *ξ*(2) on the "inner" boundary at the last iteration.

**Figure 5.** Sea level function on Γ*in* at different iterations (1—*ξ*(1) at the first iteration; 2—*ξ*(2) at the first iteration; 3—*ξ*(1) at the last iteration; 4—*ξ*(2) at the last iteration).

In addition, we study the dependence on the regularization parameter *α*. According to the theory [27], an increase in the regularization parameter *α* leads to a decrease in the influence of additional noise on the solution. Table 2 provides the error

$$\|\|\boldsymbol{\xi}^{(1)} - \boldsymbol{\xi}\_{\mathcal{P}}\|\|\_{L^{\mathcal{N}}\_{\boldsymbol{\xi}}\left(\Gamma\_{\mathcal{P}\boldsymbol{\rho}}\right)} = \left(\int\_{\Gamma\_{\mathcal{P}}} \sqrt{\boldsymbol{g}\boldsymbol{H}} \left(\boldsymbol{\xi} - \boldsymbol{\xi}\_{\mathcal{P}}\right)^{2} d\boldsymbol{\Gamma}\right)^{1/2}$$

depending on the regularization parameter *α* and the stopping criterion (the number of iterations less than 10 or less than 50). As expected, the error increases at the beginning, because the accuracy of the method is not sufficient to reproduce the additional noise (see Figure 4 and Table 1). When the stopping criterion is the number of iterations less than 50, the regularization parameter smooth the result.

**Table 2.** The error *ξ*(1) − *ξ <sup>p</sup>LW*2 (<sup>Γ</sup>*op* ) depending on the stopping criterion and regularization parameter *α*.


Thus, application of DDMs in variational DA problem is considered. The time required to solve DA problem with using DDM slightly increases (by 5%) in contrast to one required to solve only DA procedure. The developing of a parallel algorithm based on the described approach may decrease the computation time.
