Other schemes can be used, instead.
if opt == 1: #model-space approach
Bi = np.linalg.inv(B)
Ri = np.linalg.inv(R)
A = Bi + (H.T)@Ri@H
b = Bi@ub + (H.T)@Ri@w
ua = np.linalg.solve(A,b) #solve a linear system
elif opt == 2: #model-space incremental approach
Bi = np.linalg.inv(B)
Ri = np.linalg.inv(R)
A = Bi + (H.T)@Ri@H
b = (H.T)@Ri@(w-H@ub)
ua = ub + np.linalg.solve(A,b) #solve a linear system
elif opt == 3: #observation-space incremental approach
A = R + H@B@(H.T)
b = (w-H@ub)
ua = ub + B@(H.T)@np.linalg.solve(A,b) #solve a linear system
return ua
```
#### *3.2. Nonlinear Case*

On the other hand, if *h*(**u**) is a nonlinear function, Equation (5) implies the solution of a system of nonlinear equations. Unlike linear systems, few algorithms are available to directly solve nonlinear systems and their convergence and stability are usually questionable. Alternatively, we can use Taylor series to expand *h*(**u**) around an initial estimate of **u***a*, denoted as **u***c*, where **u***<sup>a</sup>* = **u***<sup>c</sup>* + Δ**u**. The first-order approximation of *h*(**u***a*) can be written as

$$h(\mathbf{u}\_a) \approx h(\mathbf{u}\_c) + \mathbf{D}\_h(\mathbf{u}\_c)\Delta\mathbf{u}\_\prime \tag{9}$$

and Equation (5) can be approximated as

$$\mathbf{D}\_h^T(\mathbf{u}\_\mathfrak{c})\mathbf{R}^{-1}(\mathbf{w} - h(\mathbf{u}\_\mathfrak{c}) - \mathbf{D}\_h(\mathbf{u}\_\mathfrak{c})\Delta\mathbf{u}) = \mathbf{B}^{-1}(\mathbf{u}\_\mathfrak{c} + \Delta\mathbf{u} - \mathbf{u}\_\mathfrak{b}).\tag{10}$$

Thus, the correction to the initial guess of **u***a* can be computed by solving the following system of linear equations

$$\left(\mathbf{B}^{-1} + \mathbf{D}\_{\hbar}^{T}(\mathbf{u}\_{\varepsilon})\mathbf{R}^{-1}\mathbf{D}\_{\hbar}(\mathbf{u}\_{\varepsilon})\right)\Delta\mathbf{u} = \left(\mathbf{B}^{-1}(\mathbf{u}\_{\delta} - \mathbf{u}\_{\varepsilon}) + \mathbf{D}\_{\hbar}(\mathbf{u}\_{\varepsilon})^{T}\mathbf{R}^{-1}(\mathbf{w} - h(\mathbf{u}\_{\varepsilon}))\right),\tag{11}$$

and a new guess of **u***<sup>a</sup>* is estimated as **u***<sup>c</sup>* + Δ**u**, which is then plugged back into Equation (11) and the computations are repeated until convergence is reached. Python implementation of the 3DVAR in nonlinear observation operator is presented in Listing 2 Although we only present the first order approximation of *h*(**u**), higher order expansions can be utilized for increased accuracy [27].

**Listing 2.** Implementation of the 3DVAR for with a nonlinear observation operator, using first-order approximation. ✞ ☎

```
import numpy as np
def NonLin3dvar(ub,w,ObsOp,JObsOp,R,B):