Backward Integration of Nonlinear Shallow Water Model: Part I: Solitary Rossby Waves

Abstract

The inviscid, nonlinear shallow water model developed by Sun was applied to study the inverse of equatorial Rossby solitons, which can be represented by the Korteweg–De Vries equation (KdV equation). The model was integrated forward in time, then the results were used as initial conditions for backward integration by just changing time step from positive to negative. The detailed structure, secondary circulation, and propagating speed of waves from both integrations are in good agreement with analytic solutions. The total mass, energy, and enstrophy are also well conserved. The procedure is much simpler and the results are more accurate than other backward integrations of 2D nonlinear models, which require significant modification of the model and can be contaminated by unwanted diffusion in forward–backward integrations or time-consuming iterative methods. This paper is also different from the numerical method for solving the inverse of the KdV equation.

1. Introduction

Here we applied Sun’s model [1] to study the reverse of fluid dynamics. The flow of Newtonian fluid at a low Reynolds number can be reversible; when the fluid between concentric cylinders is sheared by boundary motion that is subsequently reversed, then fluid elements return to their starting positions [2], as demonstrated by Taylor [3] in the experiment of Taylor–Couette flow [4]. It is well known that the numerical prediction of atmospheric and oceanic motions involves solving Navier–Stokes equations [5], a set of nonlinear partial differential equations based on conservations of mass, energy, and momentum, as well as equations of state [6]. Typically, equations of fluid dynamics are integrated forward in time. They require well-defined, accurate initial and boundary conditions, which cannot be fully provided by observations. While the inverse problem is important to numerical weather prediction, it is also important in mathematics, physics, geoscience, engineering, medicine, environment, and other areas. However, Navier–Stokes equations include viscosity/diffusion and nonlinear terms. Systems governed by time-reversible equations of motion often have irreversible behavior [2]. Traditionally, diffusion equation, a parabolic partial differential equation, cannot be integrated backward in time for a positive thermal conductivity [7], unless it has been converted to a different form as discussed in Sun and Sun [8]. In an inviscid fluid, Eckhardt and Hawcoat [9] and Fang et al. [10] suggested that if velocity satisfies u(−t) = −u(t), then the advected parcels reverse their path and return to their initial positions. The quadratic nonlinear terms in the equations break down reversibility [9]. Furthermore, nonlinear terms easily create numerical instability, especially when equations are integrated backward. Hence, scientists frequently apply diffusion/smoothing to filter instability. For example, it is suggested that the inverse of an 1D Burgers’ equation, ${\displaystyle \frac{\partial u}{\partial t}}=-u{\displaystyle \frac{\partial u}{\partial x}}+\nu {\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}$, can be integrated backward by changing the sign of viscosity $\nu$ [11]. Consequently, the backward integration results in a diluted, widespread source instead of a well-defined, concentrated source, because of the property of the heat equation [7]. Besides, nonlinear terms easily create nonlinear instability. The calculated initial value is smoothed twice during both forward and backward integrations by viscosity, which can be quite different from the actual initial condition. It can also introduce unexpected numerical instability [12]. In this article, we will show that we can integrate the inviscid nonlinear model forward and backward by simply changing the time step Δt from a positive to a negative. We also prove that the system returns to its original location if the change of velocity satisfies, du(−t) = −du(t).

In meteorology and oceanography, observational data may come from stations, atmospheric/oceanic soundings, ships, buoys, radar, satellite, GPS, airplanes, etc., at regular and irregular temporal intervals and spatial distributions. Furthermore, those data may be incomplete, especially over the ocean and remote areas. Hence, analyzing the data from different formats, locations, and observation periods; filling the gaps; and providing consistent initial and boundary conditions for numerical forecasting models are challenging inverse problems. Data assimilation methods using an optimal blend of observed values and a background where the first guess is given by the previous model forecast have become popular. The weights of the blend are determined by the error covariance matrix [13]. The three-dimensional variation assimilation (3D Var) and four-dimensional variation assimilation (4D Var) (in space and time [14]), based on variational analysis and optimal control theory [15], are frequently used to present the state of the atmosphere as the initial conditions for numerical weather prediction models [16,17]. Variational data assimilation takes models as constraints and calculates the gradient of cost functional with respect to initial conditions, boundary conditions, and parameters by the adjoint method, so that the solutions to the models fit the observations as accurately as possibly. This approach includes several modifications. Pu et al. [18] and Kalnay et al. [11] applied a quasi-inverse linear method to study forecast errors in the initial conditions for the National Centers for Environmental Prediction’s (NCEP) global spectral model. Their inverse is approximated by running the tangent linear model (TLM) of the nonlinear forecast model with a negative time step but reversing the sign of friction and diffusion terms to avoid the computational instability that would be associated with these terms if they were run backward, as discussed in the 1D Burgers’ equation previously. They also found backward integration is more unstable than forward integration, and the results can be contaminated by double diffusion during forward and backward integrations. Nabi et al. [19] developed the One-Shot method to solve the coupled system of direct, adjoint, and design equations in one procedure to cut down on computing of direct-adjoint-looping (DAL); Simple Backward Integration (SBI) was applied by Fang et al. [10] to study the energy cascade of turbulence. O’Connor et al. [20] used Simple Backward Integration (SBI) and the Iterative Quasi-Reversible Method (QRM) to simulate the inverse of the Korteweg–de Vries–Burgers (KdVB) equation.

VAR is used to obtain initial, boundary conditions and parameters of models statically or dynamically. Generally, VAR is an ill-posed, inverse problem [17]. To overcome the difficulty of ill-posedness, an additional term is usually added to the cost functional as a stabilized functional, for example, O’Connor et al. [20] added two more terms to the equation in Simple Backward Integration (SBI) to study the inverse of the KdVB.

Other approaches have also been applied to solve the inverse problems. Keller [21] used a pair of stream functions downstream to predict the time-independent upstream flow. Kaus and Podadchikov [22] applied a combined spectral/finite-difference scheme to study the Rayleigh–Taylor instability at the interface separating two fluids. They were able to reproduce the initial 2D pattern from a deformed 3D structure when the amplitude of perturbations remained small, and the system is still controlled by the linear terms of equations. The authors of [8] formulated the diffusion of a heat equation as the exchange of substances carried by subgrid ${\rho }_m^j$ (the m^th specific in the j^th cell with a fixed volume V_j) due to the subgrid-turbulent velocities $u_{i,j}^{\prime }$ moving across the surfaces A_i,j between the j^th and i^th cells becomes:

$${\displaystyle \frac{\partial {\rho }_m^j}{\partial t}}={\displaystyle \frac{1}{V_j}}\left\{-{\rho }_m^j{\displaystyle \sum _{i surround of j}{A}_{i,j}{u}_{i,j}^{\prime }}+{\displaystyle \sum _{i surround of j}{\rho }_m^i{A}_{i,j}{u}_{i,j}^{\prime }}\right\}$$

(1)

Each cell contained a different substance with high density initially at t = 0. Subgrid turbulence transferred and mixed those substances among the cells. After 4.6 days of forward integration, each cell was filled with the mixed low-density substances, which was used as the initial condition for backward integration by changing the time step from positive to negative in the same model. After 4.6 days of backward integration, we precisely reproduced the initial distribution of specifics, i.e., each cell contains a specific substance with high density as its initial value, according to Equation (1). Because the subgrid turbulence can simulate the counter gradient flow [23], different substances were extracted separately from those mixed substances in each cell, then assembled and put back to their originated cell. This avoids double counts of diffusion used in the traditional forward–backward integrations [11] The subgrid turbulence pollution model also invalidates the popular Lagrangian backward-trajectory method used to assess the property of the source from downstream, an important inverse problem in meteorology and pollution [23].

Sun and Sun [24] integrated Sun’s 2011 model forward in time to reproduce the westward propagating solitary Rossby wave. Here, we will apply the same model to study the inverse of Rossby solitons, which can be represented by the Korteweg–de Vries (KdV) equation [25,26,27,28,29],

$${\displaystyle \frac{\partial \varphi }{\partial t}}+{\displaystyle \frac{{\partial }^3\varphi }{\partial x^3}}+6\varphi {\displaystyle \frac{\partial \varphi }{\partial x}}=0$$

(2)

The solution of the KdV is

$$\varphi (x,t)=-{\displaystyle \frac{1}{2}}c{\mathrm{sech}}^2\left[{\displaystyle \frac{\sqrt{c}}{2}}(x-ct-a)\right]$$

(3)

where c is phase speed and $a$ is an arbitrary constant. The soliton propagates at a fixed phase speed, without a change in shape, through a delicate balance between linear wave dynamics and weak nonlinearity. The inverse of the KdV has been studied by Muccino and Bennett [30] by minimizing a quadratic cost functional. O’Connor et al. [20] applied the SBI by adding extra terms and the QRM to solve the inverse of a 1D KdV with viscosity.

2. Basic Equations and Numerical Model

The flux form of the shallow water equations is

$${\displaystyle \frac{\partial h}{\partial t}}=-{\displaystyle \frac{\partial \psi }{\partial x}}-{\displaystyle \frac{\partial \phi }{\partial y}}$$

(4)

$${\displaystyle \frac{\partial \psi }{\partial t}}+{\displaystyle \frac{\partial \left(u\psi \right)}{\partial x}}+{\displaystyle \frac{\partial \left(v\psi \right)}{\partial y}}-f\varphi =-gh{\displaystyle \frac{\partial (h+h_s)}{\partial x}}$$

(5)

$${\displaystyle \frac{\partial \phi }{\partial t}}+{\displaystyle \frac{\partial \left(u\phi \right)}{\partial x}}+{\displaystyle \frac{\partial \left(v\phi \right)}{\partial y}}+f\psi =-gh{\displaystyle \frac{\partial (h+h_s)}{\partial y}}$$

(6)

where h is the water depth and h_s is the terrain surface height. Thus, the free surface height can be determined using h_f = h + h_s, ψ = hu, and ϕ = hv, where u and v are x and y component velocities, respectively. Equations are solved using Sun’s scheme [1], which not only has neutral stability but also does not require physical or artificial smoothing/diffusion to control numerical instability. It has been applied to the shockwave in dam-breaking [1], vortices merging [31], etc., as discussed in [32].

3. Backward Integration

The model can be applied to both forward and backward integrations. The procedure of backward integration is the same as forward integration. We calculate $h_{p,q}^n$ from a given data set at the (n + 1)^th time step as follows:

$${\displaystyle \frac{h_{p,q}^n-h_{p,q}^{n+1}}{\left(-\Delta t\right)}}=-{\left[{\displaystyle \frac{\partial \left(\psi +\left(-\frac{\Delta t}{2}\right)\frac{\partial \psi }{\partial t}\right)}{\partial x}}\right]}_{p,q}^{n+1}-{\left[{\displaystyle \frac{\partial \left(\phi +\left(-\frac{\Delta t}{2}\right)\frac{\partial \phi }{\partial t}\right)}{\partial y}}\right]}_{p,q}^{n+1}$$

(7)

Then, we use $h_{p,q}^{n+1/2}=0.5(h_{p,q}^n+h_{p,q}^{n+1})$, ${\psi }_{p*,q}^{n+1}$ and ${\phi }_{p,q*}^{n+1}$ to calculate ${\psi }_{p*,q}^{*}$ and ${\phi }_{p,q*}^{*}$, i.e.,

$$\begin{array}{l}{\displaystyle \frac{{\psi }_{p*,q}^{*}-{\psi }_{p*,q}^{n+1}}{\left(-\Delta t\right)}}=-{\left[{\displaystyle \frac{\partial {\left(u\psi \right)}^{n+1/2}}{\partial x}}\right]}_{p*,q}-{\left[{\displaystyle \frac{\partial {\left(v\psi \right)}^{n+1/2}}{\partial y}}\right]}_{p*,q}-g{\left[h{\displaystyle \frac{\partial \left(h+h_s\right)}{\partial x}}\right]}_{p*,q}^{n+1/2}\\ =-{\displaystyle \frac{\partial }{\partial x}}{\left[u\psi +\left(-{\displaystyle \frac{\Delta t}{2}}\right){\displaystyle \frac{\partial \left(u\psi \right)}{\partial t}}\right]}_{p*,q}^{n+1}-{\displaystyle \frac{\partial }{\partial y}}{\left[v\psi +\left(-{\displaystyle \frac{\Delta t}{2}}\right){\displaystyle \frac{\partial \left(v\psi \right)}{\partial t}}\right]}_{p*,q}^{n+1}-g{\left[h^{n+1/2}{\displaystyle \frac{\partial \left(h^{n+1/2}+h_s\right)}{\partial x}}\right]}_{p*,q}\end{array}$$

(8)

$$\begin{array}{l}{\displaystyle \frac{{\phi }_{p,q*}^{*}-{\phi }_{p,q*}^{n+1}}{\left(-\Delta t\right)}}=-{\left[{\displaystyle \frac{\partial {\left(u\phi \right)}^{n+1/2}}{\partial x}}\right]}_{p,q*}-{\left[{\displaystyle \frac{\partial {\left(v\phi \right)}^{n+1/2}}{\partial y}}\right]}_{p,q*}-g{\left[h{\displaystyle \frac{\partial \left(h+h_s\right)}{\partial y}}\right]}_{p,q*}^{n+1/2}\\ =-{\displaystyle \frac{\partial }{\partial x}}{\left[u\phi -{\displaystyle \frac{\Delta t}{2}}{\displaystyle \frac{\partial \left(u\phi \right)}{\partial t}}\right]}_{p,q*}^{n+1}-{\displaystyle \frac{\partial }{\partial y}}{\left[v\phi -{\displaystyle \frac{\Delta t}{2}}{\displaystyle \frac{\partial \left(v\phi \right)}{\partial t}}\right]}_{p,q*}^{n+1}-g{\left[h^{n+1/2}{\displaystyle \frac{\partial \left(h^{n+1/2}+h_s\right)}{\partial y}}\right]}_{p,q*}\end{array}$$

(9)

Finally, we solve the Coriolis force implicitly:

$${\displaystyle \frac{{\psi }^n-{\psi }^{*}}{-\Delta t}}={\displaystyle \frac{f\left({\phi }^{*}+{\phi }^n\right)}{2}}$$

(10)

and

$${\displaystyle \frac{{\phi }^n-{\phi }^{*}}{-\Delta t}}=-{\displaystyle \frac{f\left({\psi }^{*}+{\psi }^n\right)}{2}}$$

(11)

4. Initial Conditions and Numerical Results

4.1. Initial Condition of Solitary Rossby Wave in Forward Integration

After defining s = x − ct, the nonlinear shallow water wave equations on the equatorial beta plane can be transformed into a frame of reference moving with the linear wave. The zero-order of flow variables (u, v, h) for the mode-1 wave are

$$u_0(s,y,t)=\zeta (s,t){\displaystyle \frac{\left(-9+6y^2\right)}{4}}{e}^{-0.5y^2}$$

(12)

$$v_0(s,y,t)=2y{\displaystyle \frac{\partial \zeta (s,t)}{\partial s}}{e}^{-0.5y^2}$$

(13)

$$h_0(s,y,t)=\zeta (s,t){\displaystyle \frac{\left(3+6y^2\right)}{4}}{e}^{-0.5y^2}$$

(14)

and the variable ϛ(s,t) satisfies

$${\displaystyle \frac{\partial \varsigma }{\partial t}}-1.5366\varsigma {\displaystyle \frac{\partial \varsigma }{\partial s}}-0.098765{\displaystyle \frac{{\partial }^3\varsigma }{\partial s^3}}=0$$

(15)

Like Equations (2) and (3), Equation (15) is the KdV equation with the exact solution [28]:

$$\zeta (s,t)=0.772B^2{\sec \mathrm{h}}^2(B(s+0.395B^{2}t))$$

(16)

Boyd [29] provided the sum of zeroth- and first-order solutions, i.e., $\phi ={\phi }_o+{\phi }_1$. The first-order solution is

$$u_1(s,y,t)=c_1\zeta (s,t){\displaystyle \frac{9\left(3+2y^2\right)}{16}}{e}^{-0.5y^2}+{\zeta }^2{U}_1(y)$$

(17)

$$v_1(s,y,t)=\zeta (s,t){\displaystyle \frac{\partial \zeta (s,t)}{\partial s}}{V}_1$$

(18)

and

$$h_1(s,y,t)=c_1\zeta (s,t){\displaystyle \frac{9\left(-5+2y^2\right)}{16}}{e}^{-0.5y^2}+{\zeta }^2(s,t)H_1(y)$$

(19)

where U₁, V₁, and H₁ are given by the infinite Hermite series as shown in Table 1 in [29]. We also impose

$${\displaystyle \frac{\max |v_1/(\partial \zeta /\partial s)|}{\max |v_0/(\partial \zeta /\partial s)|}}=0.277B^2$$

(20)

$${\displaystyle \frac{\max |u_1/\zeta |}{\max |u_0/\zeta |}}=0.334B^2$$

(21)

and

$${\displaystyle \frac{\max |h_1/\zeta |}{\max |h_0/\zeta |}}=0.924B^2$$

(22)

The initial condition is given by (12)–(22) at t = 0. The periodic condition is applied to the east–west boundaries, and the open boundary at the north and south. The speed of the soliton is C_B= −1/3 − 0.395 × B². The domain consists of 480 × 240 grids with the spatial intervals Δx = Δy = 0.1. The east–west length of the domain is 48. It takes 120 time units for the soliton with B = 0.395 and C_B= 0.395 to propagate westward and complete a cycle; and 100 time units for the wave with B = 0.6 and C_B= 0.476 to complete a cycle and return to the original location. We simulated height with the initial condition containing both zeroth and the first orders. More detailed forward integrations from the nonlinear model with different numerical schemes or viscosity and from the linearized equations can be found in [24].

4.2. Result with B = 0.395 and C_B = 0.395

Figure 1a shows the simulated height and location of the soliton at the beginning, x = 0 and t = 0 (black), and at t = 30 (red), 60 (green), 90 (blue) and 120 (purple) from forward integration. The soliton propagates westward and returns to its initial location because of a periodic boundary. The front of the simulated soliton at t = 120 matches well with the analytic solution shown in black, which represents waves at both t = 0 and t = 120. But the tail of the simulated wave slightly lags behind the analytic solution due to the slower phase speed of short waves. Simulations can be improved by using more accurate advection schemes [33,34], as shown in [24].

Using forward integration at t = 120 (purple in Figure 1a or black in Figure 1b) as the initial condition, the Rossby wave moves eastward as time decreases, initially at t = 120 (black), t = 90 (red), t = 60 (green), t = 30 (blue), and t = 0 (purple), as shown in Figure 1b, from backward integration. The heights and circulations after 120 time units from backward integration at t = 0 (Figure 2b) are quite comparable to its initial conditions at t = 120 in forward integration (Figure 2a). The soliton does return to its original position, but the velocity does not follow u(−t) = −u(t) as suggested by Eckhardt and Hascoët [9] or Fang et al. [10]. We conclude that the system is reversable if the change in velocity and height satisfies d ξ(−t) = −d ξ(t), which is consistent with Equations (7)–(11). The results also prove that the model without modification can be applied to both forward and backward integrations.

The total mass (Ma), energy (Eng), potential vorticity (PV), and enstrophy (He) [31] are

$$Ma(t)={\displaystyle \sum _{j=1}^{Jm}{\displaystyle \sum _{i=1}^{\mathrm{Im}}{h}_{i,j}(t)\Delta x\Delta y}}$$

(23)

$$Eng(t)={\displaystyle \frac{1}{2}}{\displaystyle \sum _{j=1}^{Jm}{\displaystyle \sum _{i=1}^{\mathrm{Im}}\left[h_{i,j}(t)\left(u_{i,j}^2+v_{i,j}^2\right)+g{\left(h_{i,j}+h{s}_{i,j}\right)}^2\right]\Delta x\Delta y}}$$

(24)

$$PV(t)={\displaystyle \sum _{j=1}^{Jm}{\displaystyle \sum _{i=1}^{\mathrm{Im}}\Pi \Delta x\Delta y}}={\displaystyle \sum _{j=1}^{Jm}{\displaystyle \sum _{i=1}^{\mathrm{Im}}\left\{\left[\left(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}\right)+f\right]/h\right\}\Delta x\Delta y}}$$

(25)

and

$$He(t)={\displaystyle \frac{1}{2}}{\displaystyle \sum _{j=1}^{Jm}{\displaystyle \sum _{i=1}^{\mathrm{Im}}{h}_{i,j}{q}_{i,j}^2\Delta x\Delta y}}={\displaystyle \frac{1}{2}}{\displaystyle \sum _{j=1}^{Jm}{\displaystyle \sum _{i=1}^{\mathrm{Im}}\frac{{\left({\varsigma }_{i,j}+f_{i,j}\right)}^2}{h_{i,j}}\Delta x\Delta y}}$$

(26)

Figure 3a shows the time sequences of the ratio of the variable with respect to its initial value, Ma(t)/Ma(t = 0) by green, Eng(t)/Eng(t = 0) by red, PV(t)/PV(t = 0) by blue, and He(t)/He(t = 0) by black according to (23)–(26) from forward integration. Figure 3b shows the same, except it is from backward integration. They remain almost constant with respect to time in both integrations, although an open boundary condition is applied to the northern and southern boundaries, because the wave is trapped near the equator.

4.3. Results with B = 0.6 and C_B = 0.476

The soliton propagates westward from forward integration as shown in Figure 4a: red designates the location at t = 20, green at t = 40, blue at t = 60, light blue at t = 80, and purple at t = 100. Then, the result at t = 100 is used as an initial condition for backward integration. Red in Figure 4b shows the locations at t = 80, green at t = 60, blue at t = 40, light blue at t = 20, and purple at t = 0, as the soliton moves backward (eastward) as time decreases in backward integration. The propagating speed of the soliton is faster than in the previous case. The change in heights between t = 0 (black) and t = 100 (purple) during forward integration in Figure 4a and between t = 100 (black) and t = 0 (purple) in backward integration (Figure 4b) is larger than in the previous case with B = 0.395.

The KdV equation was derived based on asymptotic expansion of small perturbation i.e., $\epsilon =O(\varsigma /H(\mathrm{mean} \mathrm{height}))<<1$. According to (16) with B = 0.6, we obtain $\epsilon =0.278$, which is consistent with the perturbation of water height shown in Figure 5a,b, but does not satisfy $\epsilon <<1$. Hence, it needs higher-order terms to represent the KdV solution, as discussed in [29,35,36,37] and in shallow water modeling [32].

The secondary circulations with respect to the moving soliton, (u-C_B, v), are like the simulations of Sun and Sun [24] and discussed in [27,28]. It is also noted that the longer the integration, the wider the spreading between long waves and short waves in advection equations [34]. Hence, distortion becomes visible after 100 time units from backward integration at t = 0 (purple) in Figure 4b. The error of phase speed can be reduced by using a more accurate advection scheme [34], the scheme combining the QUICK [38] and fourth-order schemes as shown in the simulations of Sun and Sun [24], or semi-Lagrange schemes [39,40]. The time sequences of the ratios, Ma(t)/Ma(t = 0) (green), Eng(t)/Eng(t = 0) (red), PV(t)/PV(t = 0) (blue), and He(t)/He(t = 0) (black) remain near unity in both forward and backward integrations (Figure 6a,b). The wave is trapped around the equator and the slight decrease in the total energy comes from the truncational error of the fourth-order scheme and the dispersion of phase speed among different waves.

5. Discussion

In this paper, the analytic solutions of solitary Rossby waves are used to verify the backward integrations because the shape and propagating speed of the waves do not change with time. We will present the inverse of vortices merging [31] and nonlinear shear instability [41] in Part II, in which the nonlinear terms become dominated. In the future, we will combine this inviscid nonlinear model and the subgrid turbulence diffusion models [8,23] to study the inverse of fluid dynamic systems with diffusion/viscosity. The results should be close to realistic initial values, because the subgrid turbulence diffusion model can overcome the difficulty of counter-gradient flow. It also avoids extra diffusion by changing the sign of viscosity in the traditional forward–backward integrations or the time-consuming iterative data assimilations discussed previously.

6. Summary

The shallow water model developed by Sun [1] was applied to simulate Rossby solitons in the equator. The model was integrated forward in time for 120 time units (with B = 0.395) or 100 time units (with B = 0.6), then the simulations were used as initial conditions for backward integrations. The detailed structure, secondary circulation, and propagating speed of solitons from both integrations are in good agreement with the analytic solutions of KdV. The total mass, energy, and enstrophy are also well conserved, although a slight decrease in total kinetic energy may be due to truncational errors of the fourth-order scheme and the dispersion of phase speed among different waves.

This inviscid, nonlinear model, which has neutral stability and does not need any viscosity/smoothing to control instability, can be applied to both forward and backward integrations by just changing the time step from positive to negative. Hence, not only is the procedure much simpler, but the results are also more accurate than the traditional inverse methods which can be contaminated by unwanted diffusion in forward–backward integrations, or inaccurate guesses of initial values in time-consuming iterative data assimilations.