# Optimization of Open Boundary Conditions in a 3D Internal Tidal Model with the Adjoint Method around Hawaii

Authors: Anzhou Cao; Haibo Chen; Jicai Zhang; Xianqing Lv
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

Based on the theory of inverse problem, the optimization of open boundary conditions (OBCs) in a 3D internal tidal model is investigated with the adjoint method. Fourier coefficients ofM2 internal tide on four open boundaries, which are regarded as OBCs, are inverted simultaneously. During the optimization, the steepest descent method is used to minimize cost function. The reasonability and feasibility of the model are tested by twin experiments (TEs). In TE1, OBCs on four open boundaries are successfully inverted by using independent point (IP) strategy, suggesting that IP strategy is useful in parameter estimation. Results of TE2 indicate that the model is effective even by assimilating inaccurate “observations.” Based on conclusions of TEs, the M2 internal tide around Hawaii is simulated by assimilating T/P data in practical experiment.  The simulated cochart shows good agreement with that obtained from the Oregon State University tidal model and T/P observations. Careful inspection shows that the major difference between simulated results and OSU model results is short-scale fluctuations superposed on coamplitude lines, which can be treated as the surface manifestation modulated by the internal tide. The computed surface manifestation along T/P tracks is comparable to the estimation in previous work.


---

## Body

## 1. Introduction

Internal tide is an important intermediate step of the tides-to-turbulence cascade [1] and plays an important role in the dissipation of surface tides and mixing in the deep ocean so as to have strong influence on the global thermohaline circulation [2]. The Hawaiian Ridge is one of the most suitable areas to investigate the internal tides. Rudnick et al. [1], Ray and Mitchum [3, 4], Munk [5], Martin et al. [6], and Zhao et al. [7–9] investigated internal tide along the Hawaiian Ridge by using data of T/P satellite altimetry, Hawaii Ocean Mixing Experiment (HOME), and many other observations such as acoustics Doppler current profilers (ADCPs). Kang et al. [10], Simmons et al. [11], Carter et al. [12], Powell et al. [13], Rainville et al. [14], and Zaron and Egbert [15] used numerical models to simulate internal tides and study the mechanism of internal tide. Through their work, the generation, propagation, vertical structure, energy budget, and other characteristics of internal tide around Hawaii have been described preliminarily.Among all data assimilation methods, the 4D variational data assimilation is one of the most effective approaches and widely used in meteorological and oceanographic models [16–32]. Based on the theory of inverse problem, it is an advanced data assimilation method which involves the adjoint method and has the advantage of assimilating various observations disturbed in time and space [22–37]. Lardner [22], Myers and Baptista [23], Zhang and Lu [24], Guo et al. [25], and Cao et al. [26] studied the inverse problem of OBCs in barotropic tidal models. Chen et al.’s [27] constructed a 3D isopycnic-coordinate internal tidal model with adjoint method and studied inverse problem of OBCs at 14 kinds of bottom topographies. As for the inverse problem of other parameters, Lu and Zhang [28] and Zhang et al. [29] inverted spatially varying bottom friction coefficient in the 2D tidal model and Yu and O’Brien [30] estimated wind drag coefficient and eddy viscosity in the Ekman layer model.The work of this paper is an extension of Chen et al. [27], in which Fourier coefficients of M2 internal tide on four open boundaries are inverted simultaneously. During the optimization, steepest descent method with special step factor is used to minimize cost function. This paper is organized as follows. The 3D internal tidal model with adjoint method, minimization algorithm, independent point (IP) strategy and model settings are introduced in Section 2. Identical twin experiments (TEs) are carried out to test the model in Section 3. Practical experiment (PE) and its analysis are given in Section 4. Finally, the summary and conclusions are provided in Section 5.


## 2. Numerical Model and Settings

### 2.1. Forward Model

The governing equations of internal tide in spherical coordinate are described as follows:(1)∂qk∂t+1a∂(qkuk)∂λ+1a∂(qkvkcosφ)∂φ=0,∂uk∂t+L(uk)-fvk=-1a∂Mk∂λ-ga∂ζ-∂λ+τ(k-1)x-τkxq1+AhΔuk,∂vk∂t+L(vk)+fuk=-1r∂Mk∂φ-gr∂ζ-∂φ+τ(k-1)y-τkyq1+AhΔvk,
in which,
(2)L(uk)=uka∂uk∂λ+vkr∂uk∂φ-ukvktanφr,L(vk)=uka∂vk∂λ+vkr∂vk∂φ+uk2tanφr,Δuk=1a2∂2uk∂λ2+1ra∂∂φ(∂uk∂φcosφ),Δvk=1a2∂2vk∂λ2+1ra∂∂φ(∂vk∂φcosφ),M1=g(qρ-D),(3)Mk=g(qρ-D)+g∑m=1k-1[(1ρk-1ρm)qm],
where k=1,2,…,N is vertical layer index, t is time, λ and φ are longitude and latitude, qk, uk and vk are water mass and horizontal velocity of the kth layer, respectively, f is the Coriolis parameter, ζ- is the tidal potential, Ah is the horizontal viscosity coefficient, r is the radius of the earth, a=rcosφ, g is the gravity, D is the total depth at rest, q=∑k=1Nqk is the total mass of the water column, ρ=q/∑k=1N(qk/ρk) is the mean density, τ0x and τ0y are wind stress components (in the paper, they are equal to zero), τNx and τNy are bottom friction components (τNx=κρNuNuN2+vN2, τNy=κρNvNuN2+vN2 in which κ is the nondimensional bottom friction coefficient), and τkx and τky (k=1,…,N-1) are interfacial friction between the kth layer and the (k+1)th layer, and they are defined as
(4)τkx=μk(uk-uk+1),τky=μk(vk-vk+1),μk=Avρk+ρk+1hk+hk+1,
where Av is the vertical viscosity coefficient; hk is the undisturbed thickness of the kth layer.Equations (1) and its vertical integration, which are called the internal and external mode, respectively, compose the forward model and are used to simulate the internal tide. Miao et al. [31] used this model to investigate the characteristic of generation and propagation of M2 and K1 internal tides in the South China Sea.The Flather condition [38] can yield good results [27, 31] when applied to the open boundaries in this model. In this paper, an adaptation of the Flather condition is installed in the external mode of the forward model, which is given as follows:
(5)ζ=ζ~±(1-f2ω2)HgU,ζ=ζ~±(1-f2ω2)HgV,
where ζ~ is external data beyond the model boundary representing the clamped surface elevation relating to the boundary barotropic tidal force, H is the time-varying total water depth, and U and V are barotropic velocity. The signs in (5) are positive for eastern and northern boundaries and negative for western and southern boundaries. Equations (5) can be determined once the external data ζ~ is described. The M2 tidal force ζ~ at the nth time step along the open boundary can be described as
(6)ζ~ii,jjn=aii,jjcos(ωnΔt)+bii,jjsin(ωnΔt),
where a and b are Fourier coefficients of OBCs, ω is the frequency of M2 tide, and Δt is time step. So the inverse problem of OBCs is turned into the inverse problem of Fourier coefficients a and b.For the internal mode, the relaxation conditions are applied. A relaxation region of several grids close to the boundary is defined. In this region, tidal values including the interface elevation and fluid velocity of both barotropic and baroclinic tides are first computed and then relaxed to final values displayed as follow(7)x(i,j)=xBt(i,j)+[1-γ(i,j)]·xBc(i,j),
where x represents the tidal state ηk, uk, or vk, xBt and xBc represent the barotropic and baroclinic tides, respectively, and γ(i,j)∈[0,1] is the relaxation factor in grid (i,j). In our model, γ is increased from 0 to 1 as the grid is getting close to the open boundary.In addition, the initial conditions are that the initial surface elevation, interfacial elevation,x- and y-velocity at all layers are zeros.


### 2.2. Adjoint Model

The adjoint model is similar to that in Chen et al. [27], which can be used to optimize OBCs and/or other spatially varying parameters and displayed as follows.The cost function is defined by(8)J(q,u,v;p)=12∫Σ(qρ-D-ζ^)2dσ,
where Σ is the model domain of both time and space, ζ^ is the observations, and p represents the control variables and is defined by the OBCs, which are denoted by a series of Fourier coefficients in this work.The Lagrangian function is defined by(9)L(q,u,v,qa,ua,va;p)=J(q,u,v;p)+∫Σ[qa1(+1a∂(q1v1cosφ)∂φ∂q1∂t+1a∂(q1u1)∂λ=J(q,u,v;p)(q,u,v)+∫Σ+1a∂(q1v1cosφ)∂φ)]dσ+∫Σ[ua1(+ga∂ζ-∂λ-τ0x-τ1xq1-AhΔu1∂u1∂t+L(u1)-fv1+1a∂M1∂λ+∫Σ[ua1+ga∂ζ-∂λ-τ0x-τ1xq1-AhΔu1)]dσ+∫Σ[va1(+gr∂ζ-∂φ-τ0y-τ1yq1-AhΔv1∂v1∂t+L(v1)+fu1+1r∂M1∂φ+[(1r∂M1∂φ+gr∂ζ-∂φ-τ0y-τ1yq1-AhΔv1)]dσ+⋯+∫Σ[qaN(∂qN∂t+1a∂(qNuN)∂λ+1a∂(qNvNcosφ)∂φ)]dσ+∫Σ[uaN(+ga∂ζ-∂λ-τ(N-1)x-τNxq1-AhΔuN∂uN∂t+L(uN)-fvN+1a∂MN∂λ+∫Σ[(∂MN∂λ+ga∂ζ-∂λ-τ(N-1)x-τNxq1-AhΔuN)]dσ+∫Σ[vaN(+gr∂ζ-∂φ-τ(N-1)y-τNyq1-AhΔvN∂vN∂t+L(vN)+fuN+1r∂MN∂φ+∫Σ[(1r∂MN∂φ+gr∂ζ-∂φ-τ(N-1)y-τNyq1-AhΔvN)]dσ,
where qak, uak and vak are adjoint variables of qk, uk, and vk, respectively.According to the typical theory of Lagrangian multiplier method, we have the following first-order derivates of Lagrangian function with respect to all the variables and parameters:(10)∂L∂qak=0,∂L∂uak=0,∂L∂vak=0,(11)∂L∂qk=0,∂L∂uk=0,∂L∂vk=0,(12)∂L∂p=0.Equations (10) are (1) actually. The adjoint equations can be derived from (11). From (12), the gradients of the cost function with respect to control variables can be obtained and displayed as follows:
(13)∂J∂aii,jj+∑nTii,jjncos(ωnΔt)=0,∂J∂bii,jj+∑nTii,jjnsin(ωnΔt)=0.
Details on computing Tii,jjn are given in [27].


### 2.3. Minimization Algorithm

Considering that the aim is to make the simulated results closer to the observations, there is the optimization problem as follows:(14)min:J(x),
where J is the cost function; x represents Fourier coefficients a and b. Suppose that the dimension of x is N and the gradient of J is ∇J(xik) with respect to the ith parameter (1≤i≤N) and the kth iteration step. According to the typical steepest descent method,
(15)xik+1=xik+αk·pik,
where pik is the negative of ∇J(xik); αk is the step factor at the kth iteration step. The initial parameter xi0 should be prescribed (in the following experiments, xi0=0). In order to minimize J, the step factor αk used in this paper is performed as follows:
(16)αk=α0∥∇J(xi1)∥2,
where α0 is a constant which can be obtained by prior experience; ∥∇J(xi1)∥2 is the L2 norm of gradient at the first iteration step. Then formula (6) is changed into
(17)xik+1=xik-α0·∇J(xik)∥∇J(xi1)∥2orxik+1=xik-α0·∥∇J(xik)∥2∥∇J(xi1)∥2·∇J(xik)∥∇J(xik)∥2,
where ∇J(xik)/∥∇J(xik)∥2 is the normalized gradient at the kth iteration step and ∥∇J(xik)∥2/∥∇J(xi1)∥2 is the ratio between L2 norm of gradient at the kth iteration step and that at the first iteration step, which could be a decreasing value as iteration step is increased.


### 2.4. Independent Point Strategy

Independent point strategy is considered to be a useful tool in the estimation of OBCs and spatially varying parameters. On one hand, adopting IP strategy could decrease the the number of variables so as to increase the well posedness of inverse problem. It is generally accepted that number of observations should be larger than that of variables and if number of observations is almost equal to or smaller than that of variables, the inverse problem will become ill-posed and the inverse results may be bad. On the other hand, adopting IP strategy could ensure the continuity of variables, which accords with the physical properties of variables. Lu and Zhang [28] studied the inversion problem of spatially varying bottom frication coefficients by using adjoint method and IP stragety. Results indicated that adopting IP strategy in which IPs are distributed according to depth could lead to a better simulation. Based on the theory of inverse problem, Zhang and Lu [24] simulated 3D tidal current in the Bohai and Yellow Sea by assimilating the satellite altimetry. They found that too large number of variables would lead to inaccurate inversion results, while too small number of IPs could not figure out the real condition of ocean. Guo et al. [25] and Cao et al. [26] investigated different IP strategies in a 2D tidal adjoint model. They tested the sensitivity of the model to the IP distribution, interpolation radius, and relation between prescribed OBCs and IP strategies.IP strategy can be described as follows. It is assumed that several IPs exist along the open boundary (corresponding to the inversion of OBCs) or in the domain area (corresponding to the inversion of spatially varying parameters). Values at the IPs are optimized with adjoint method and those at other points are determined by linearly interpolating the values at IPs.


### 2.5. Settings

The area of interest is around Hawaii, from 196°E to 206°E and from 17°N to 25°N. Figure1 shows the bathymetry and T/P tracks in this area. In this model, horizontal resolution is 10′×10′ and time step is one sixtieth of the M2 tidal period. According to area-averaged potential density and buoyancy frequency calculated from temperature and salinity data of World Ocean Atlas 2005 (WOA05), six uneven vertical layers are set varying from 100 m at the top to more than 4000 m at the bottom (shown in Figure 2). The initial density is set to be horizontally homogeneous. The non-dimensional bottom friction coefficient κ, horizontal eddy viscosity Ah, and vertical eddy viscosity Av are 0.003, 1000 m2/s, and 0.005 m2/s, respectively.Figure 1
Bathymetry and T/P tracks around Hawaii (“o” represents grid point on the open boundary).Area-averaged potential density and buoyancy frequency.
(a)
(b)


## 2.1. Forward Model

The governing equations of internal tide in spherical coordinate are described as follows:(1)∂qk∂t+1a∂(qkuk)∂λ+1a∂(qkvkcosφ)∂φ=0,∂uk∂t+L(uk)-fvk=-1a∂Mk∂λ-ga∂ζ-∂λ+τ(k-1)x-τkxq1+AhΔuk,∂vk∂t+L(vk)+fuk=-1r∂Mk∂φ-gr∂ζ-∂φ+τ(k-1)y-τkyq1+AhΔvk,
in which,
(2)L(uk)=uka∂uk∂λ+vkr∂uk∂φ-ukvktanφr,L(vk)=uka∂vk∂λ+vkr∂vk∂φ+uk2tanφr,Δuk=1a2∂2uk∂λ2+1ra∂∂φ(∂uk∂φcosφ),Δvk=1a2∂2vk∂λ2+1ra∂∂φ(∂vk∂φcosφ),M1=g(qρ-D),(3)Mk=g(qρ-D)+g∑m=1k-1[(1ρk-1ρm)qm],
where k=1,2,…,N is vertical layer index, t is time, λ and φ are longitude and latitude, qk, uk and vk are water mass and horizontal velocity of the kth layer, respectively, f is the Coriolis parameter, ζ- is the tidal potential, Ah is the horizontal viscosity coefficient, r is the radius of the earth, a=rcosφ, g is the gravity, D is the total depth at rest, q=∑k=1Nqk is the total mass of the water column, ρ=q/∑k=1N(qk/ρk) is the mean density, τ0x and τ0y are wind stress components (in the paper, they are equal to zero), τNx and τNy are bottom friction components (τNx=κρNuNuN2+vN2, τNy=κρNvNuN2+vN2 in which κ is the nondimensional bottom friction coefficient), and τkx and τky (k=1,…,N-1) are interfacial friction between the kth layer and the (k+1)th layer, and they are defined as
(4)τkx=μk(uk-uk+1),τky=μk(vk-vk+1),μk=Avρk+ρk+1hk+hk+1,
where Av is the vertical viscosity coefficient; hk is the undisturbed thickness of the kth layer.Equations (1) and its vertical integration, which are called the internal and external mode, respectively, compose the forward model and are used to simulate the internal tide. Miao et al. [31] used this model to investigate the characteristic of generation and propagation of M2 and K1 internal tides in the South China Sea.The Flather condition [38] can yield good results [27, 31] when applied to the open boundaries in this model. In this paper, an adaptation of the Flather condition is installed in the external mode of the forward model, which is given as follows:
(5)ζ=ζ~±(1-f2ω2)HgU,ζ=ζ~±(1-f2ω2)HgV,
where ζ~ is external data beyond the model boundary representing the clamped surface elevation relating to the boundary barotropic tidal force, H is the time-varying total water depth, and U and V are barotropic velocity. The signs in (5) are positive for eastern and northern boundaries and negative for western and southern boundaries. Equations (5) can be determined once the external data ζ~ is described. The M2 tidal force ζ~ at the nth time step along the open boundary can be described as
(6)ζ~ii,jjn=aii,jjcos(ωnΔt)+bii,jjsin(ωnΔt),
where a and b are Fourier coefficients of OBCs, ω is the frequency of M2 tide, and Δt is time step. So the inverse problem of OBCs is turned into the inverse problem of Fourier coefficients a and b.For the internal mode, the relaxation conditions are applied. A relaxation region of several grids close to the boundary is defined. In this region, tidal values including the interface elevation and fluid velocity of both barotropic and baroclinic tides are first computed and then relaxed to final values displayed as follow(7)x(i,j)=xBt(i,j)+[1-γ(i,j)]·xBc(i,j),
where x represents the tidal state ηk, uk, or vk, xBt and xBc represent the barotropic and baroclinic tides, respectively, and γ(i,j)∈[0,1] is the relaxation factor in grid (i,j). In our model, γ is increased from 0 to 1 as the grid is getting close to the open boundary.In addition, the initial conditions are that the initial surface elevation, interfacial elevation,x- and y-velocity at all layers are zeros.


## 2.2. Adjoint Model

The adjoint model is similar to that in Chen et al. [27], which can be used to optimize OBCs and/or other spatially varying parameters and displayed as follows.The cost function is defined by(8)J(q,u,v;p)=12∫Σ(qρ-D-ζ^)2dσ,
where Σ is the model domain of both time and space, ζ^ is the observations, and p represents the control variables and is defined by the OBCs, which are denoted by a series of Fourier coefficients in this work.The Lagrangian function is defined by(9)L(q,u,v,qa,ua,va;p)=J(q,u,v;p)+∫Σ[qa1(+1a∂(q1v1cosφ)∂φ∂q1∂t+1a∂(q1u1)∂λ=J(q,u,v;p)(q,u,v)+∫Σ+1a∂(q1v1cosφ)∂φ)]dσ+∫Σ[ua1(+ga∂ζ-∂λ-τ0x-τ1xq1-AhΔu1∂u1∂t+L(u1)-fv1+1a∂M1∂λ+∫Σ[ua1+ga∂ζ-∂λ-τ0x-τ1xq1-AhΔu1)]dσ+∫Σ[va1(+gr∂ζ-∂φ-τ0y-τ1yq1-AhΔv1∂v1∂t+L(v1)+fu1+1r∂M1∂φ+[(1r∂M1∂φ+gr∂ζ-∂φ-τ0y-τ1yq1-AhΔv1)]dσ+⋯+∫Σ[qaN(∂qN∂t+1a∂(qNuN)∂λ+1a∂(qNvNcosφ)∂φ)]dσ+∫Σ[uaN(+ga∂ζ-∂λ-τ(N-1)x-τNxq1-AhΔuN∂uN∂t+L(uN)-fvN+1a∂MN∂λ+∫Σ[(∂MN∂λ+ga∂ζ-∂λ-τ(N-1)x-τNxq1-AhΔuN)]dσ+∫Σ[vaN(+gr∂ζ-∂φ-τ(N-1)y-τNyq1-AhΔvN∂vN∂t+L(vN)+fuN+1r∂MN∂φ+∫Σ[(1r∂MN∂φ+gr∂ζ-∂φ-τ(N-1)y-τNyq1-AhΔvN)]dσ,
where qak, uak and vak are adjoint variables of qk, uk, and vk, respectively.According to the typical theory of Lagrangian multiplier method, we have the following first-order derivates of Lagrangian function with respect to all the variables and parameters:(10)∂L∂qak=0,∂L∂uak=0,∂L∂vak=0,(11)∂L∂qk=0,∂L∂uk=0,∂L∂vk=0,(12)∂L∂p=0.Equations (10) are (1) actually. The adjoint equations can be derived from (11). From (12), the gradients of the cost function with respect to control variables can be obtained and displayed as follows:
(13)∂J∂aii,jj+∑nTii,jjncos(ωnΔt)=0,∂J∂bii,jj+∑nTii,jjnsin(ωnΔt)=0.
Details on computing Tii,jjn are given in [27].


## 2.3. Minimization Algorithm

Considering that the aim is to make the simulated results closer to the observations, there is the optimization problem as follows:(14)min:J(x),
where J is the cost function; x represents Fourier coefficients a and b. Suppose that the dimension of x is N and the gradient of J is ∇J(xik) with respect to the ith parameter (1≤i≤N) and the kth iteration step. According to the typical steepest descent method,
(15)xik+1=xik+αk·pik,
where pik is the negative of ∇J(xik); αk is the step factor at the kth iteration step. The initial parameter xi0 should be prescribed (in the following experiments, xi0=0). In order to minimize J, the step factor αk used in this paper is performed as follows:
(16)αk=α0∥∇J(xi1)∥2,
where α0 is a constant which can be obtained by prior experience; ∥∇J(xi1)∥2 is the L2 norm of gradient at the first iteration step. Then formula (6) is changed into
(17)xik+1=xik-α0·∇J(xik)∥∇J(xi1)∥2orxik+1=xik-α0·∥∇J(xik)∥2∥∇J(xi1)∥2·∇J(xik)∥∇J(xik)∥2,
where ∇J(xik)/∥∇J(xik)∥2 is the normalized gradient at the kth iteration step and ∥∇J(xik)∥2/∥∇J(xi1)∥2 is the ratio between L2 norm of gradient at the kth iteration step and that at the first iteration step, which could be a decreasing value as iteration step is increased.


## 2.4. Independent Point Strategy

Independent point strategy is considered to be a useful tool in the estimation of OBCs and spatially varying parameters. On one hand, adopting IP strategy could decrease the the number of variables so as to increase the well posedness of inverse problem. It is generally accepted that number of observations should be larger than that of variables and if number of observations is almost equal to or smaller than that of variables, the inverse problem will become ill-posed and the inverse results may be bad. On the other hand, adopting IP strategy could ensure the continuity of variables, which accords with the physical properties of variables. Lu and Zhang [28] studied the inversion problem of spatially varying bottom frication coefficients by using adjoint method and IP stragety. Results indicated that adopting IP strategy in which IPs are distributed according to depth could lead to a better simulation. Based on the theory of inverse problem, Zhang and Lu [24] simulated 3D tidal current in the Bohai and Yellow Sea by assimilating the satellite altimetry. They found that too large number of variables would lead to inaccurate inversion results, while too small number of IPs could not figure out the real condition of ocean. Guo et al. [25] and Cao et al. [26] investigated different IP strategies in a 2D tidal adjoint model. They tested the sensitivity of the model to the IP distribution, interpolation radius, and relation between prescribed OBCs and IP strategies.IP strategy can be described as follows. It is assumed that several IPs exist along the open boundary (corresponding to the inversion of OBCs) or in the domain area (corresponding to the inversion of spatially varying parameters). Values at the IPs are optimized with adjoint method and those at other points are determined by linearly interpolating the values at IPs.


## 2.5. Settings

The area of interest is around Hawaii, from 196°E to 206°E and from 17°N to 25°N. Figure1 shows the bathymetry and T/P tracks in this area. In this model, horizontal resolution is 10′×10′ and time step is one sixtieth of the M2 tidal period. According to area-averaged potential density and buoyancy frequency calculated from temperature and salinity data of World Ocean Atlas 2005 (WOA05), six uneven vertical layers are set varying from 100 m at the top to more than 4000 m at the bottom (shown in Figure 2). The initial density is set to be horizontally homogeneous. The non-dimensional bottom friction coefficient κ, horizontal eddy viscosity Ah, and vertical eddy viscosity Av are 0.003, 1000 m2/s, and 0.005 m2/s, respectively.Figure 1
Bathymetry and T/P tracks around Hawaii (“o” represents grid point on the open boundary).Area-averaged potential density and buoyancy frequency.
(a)
(b)


## 3. Twin Experiments

TEs are designed to test the reasonability and feasibility of the model. The results and conclusions of TEs are useful for PE. The process of TE is designed as follow. (a) Run the forward model with prescribed OBCs. The simulated sea surface elevation recorded at the grid points on T/P satellite tracks is taken as the “observations.” Here only the positions of T/P tracks are used and the “observations” are not the real observations but the simulated results of forward model based on the prescribed OBCs. (b) Initial values of Fourier coefficients (taken as 0 here) are assigned to the OBCs to run the forward model. The values of state variables are obtained. (c) The difference of sea surface elevation between simulated results and “observations” serves as the external force of the adjoint model. The adjoint variables are obtained through the backward integration of the adjoint equations. (d) After obtaining the values of state variables and adjoint variables, the gradient of Fourier coefficients can be computed and the Fourier coefficients can be updated. Repeat (b), (c), and (d); the Fourier coefficients will be optimized and the difference of sea surface elevation between simulated results and “observations” will decrease. Once some convergence criterion is met, the process of inversion terminates.The iteration will be terminated once a convergence criterion is met. The criterion could be that the last two values of the cost function are sufficiently close, the magnitude of the gradient is sufficiently small, the discrepancy between the updated and old parameters is sufficiently small, or a combination of these [32]. In our work, the criterion of TEs is that the number of iteration steps is equal to 30 exactly.


### 3.1. Inversion of Fourier Coefficients on Four Open Boundaries

In this section, we study the inverse problem of Fourier coefficients on four open boundaries simultaneously, which means more variables are estimated compared with the work of Chen et al. [27] where only Fourier coefficient a on west open boundary is inverted. In addition, the feasibility and efficiency of IP strategy are tested in this section.Figure3 is the snapshot of prescribed OBCs and the inverse values. In TE1a, IP strategy in which IPs are distributed evenly every four grid points on all four open boundaries is employed, while in TE1b no IP strategy is used. As a whole, the shapes of inverted OBCs are similar to those of prescribed ones in both TE1a and TE1b. But it cannot be denied that inverted OBCs are disordered with fluctuations in TE1b where no IP strategy is used, implying the ill-posedness of inverse problem in TE1b and revealing feasibility and efficiency of IP strategy. Figure 4 shows the cost function and gradient norm which are normalized by their initial values at the first iteration step. As iteration step increased, the values of cost function and gradient norm are decreased obviously and at the end of iteration could reach 1% and 10% of their initial values, respectively, implying the optimizing ability of the model. In addition, as we use the step factor mentioned in Section 2, the cost function decreases efficiently without fluctuations. Tables 1 and 2 report the mean absolute errors (MAEs) and relative coefficients (RCs) between inverted OBCs and prescribed ones in TE1a and TE1b, respectively. All MAEs and RCs are less than 3 cm and greater than 0.7, respectively, for both Fourier coefficients a and b and on every open boundary, suggesting the reasonability and feasibility of the model. In addition, the inverse results in TE1a are better than those in TE1b, revealing the feasibility and efficiency of IP strategy. Combining all these results, it can be concluded that the 3D internal tidal model involving adjoint method is suitable to invert OBCs and the IP strategy is useful in the estimation of OBCs.Table 1
MAEs between inverted OBCs and prescribed ones in TE1 (unit: cm).
TE1a
TE1b
a
b
a
b
West
1.5
1.7
1.7
1.8
South
1.1
1.2
1.3
1.4
East
0.8
2.6
2.0
2.7
North
1.9
1.2
2.0
2.6Table 2
RCs between inverted OBCs and prescribed ones in TE1.
TE1a
TE1b
a
b
a
b
West
0.94
0.99
0.80
0.97
South
0.99
0.79
0.96
0.72
East
0.98
0.96
0.86
0.88
North
0.97
0.86
0.96
0.70Prescribed OBCs and inverted ones in TE1.
(a)
(b)Cost function and gradient norm in TE1.
(a)
(b)


### 3.2. Inversion by Assimilating “Inaccurate Observations”

Considering that practice in situ observations may contain noises, artificial random errors are added to “observations.” In this section, based on the “observations” in TE1a, five TEs are carried out and the maximum percentages of error are 10%, 20%, 30%, 40%, and 50% in TE2a, TE2b, Te2c, TE2d and TE2e, respectively. Then these “observations” are assimilated to invert OBCs. Results are displayed in Figures5, 6, and 7.(a) Prescribed OBCs and inverted ones on the west open boundary in TE2. (b) Prescribed OBCs and inverted ones on the south open boundary in TE2. (c) Prescribed OBCs and inverted ones on the east open boundary in TE2. (d) Prescribed OBCs and inverted ones on the north open boundary in TE2.
(a)
(b)
(c)
(d)Cost function and gradient norm in TE2.
(a)
(b)MAEs and RCs in TE2 (aw, bw, as, bs, ae, be, an, and bn represent Fourier coefficients on west, south, east, and north open boundaries, respectively).
(a)
(b)Figure5 illustrates the inverse results. All the shapes of inverse OBCs are consistent with those of the prescribed ones. It is found that with increased artificial errors, the inverted OBCs become further away from the prescribed values, implying that larger errors will result in worse inversion of OBCs. The cost function and its gradient norm are plotted in Figure 6. After assimilation, the cost function is about 1% of its initial value in TE1a and about 3% in TE2e. Combined with MAEs and RCs between inverted OBCs and prescribed ones in Figure 7, where all MAEs are less than 3.5 cm and RCs are larger than 0.6, it is reasonable to believe that the prescribed OBCs have been inverted successfully though assimilating the “observations” is not accurate. But it cannot be denied that as artificial errors increased, the inversion results become worse with larger MAEs and smaller RCs. In addition, in these TEs, the inversion results of a are better than those of b.


## 3.1. Inversion of Fourier Coefficients on Four Open Boundaries

In this section, we study the inverse problem of Fourier coefficients on four open boundaries simultaneously, which means more variables are estimated compared with the work of Chen et al. [27] where only Fourier coefficient a on west open boundary is inverted. In addition, the feasibility and efficiency of IP strategy are tested in this section.Figure3 is the snapshot of prescribed OBCs and the inverse values. In TE1a, IP strategy in which IPs are distributed evenly every four grid points on all four open boundaries is employed, while in TE1b no IP strategy is used. As a whole, the shapes of inverted OBCs are similar to those of prescribed ones in both TE1a and TE1b. But it cannot be denied that inverted OBCs are disordered with fluctuations in TE1b where no IP strategy is used, implying the ill-posedness of inverse problem in TE1b and revealing feasibility and efficiency of IP strategy. Figure 4 shows the cost function and gradient norm which are normalized by their initial values at the first iteration step. As iteration step increased, the values of cost function and gradient norm are decreased obviously and at the end of iteration could reach 1% and 10% of their initial values, respectively, implying the optimizing ability of the model. In addition, as we use the step factor mentioned in Section 2, the cost function decreases efficiently without fluctuations. Tables 1 and 2 report the mean absolute errors (MAEs) and relative coefficients (RCs) between inverted OBCs and prescribed ones in TE1a and TE1b, respectively. All MAEs and RCs are less than 3 cm and greater than 0.7, respectively, for both Fourier coefficients a and b and on every open boundary, suggesting the reasonability and feasibility of the model. In addition, the inverse results in TE1a are better than those in TE1b, revealing the feasibility and efficiency of IP strategy. Combining all these results, it can be concluded that the 3D internal tidal model involving adjoint method is suitable to invert OBCs and the IP strategy is useful in the estimation of OBCs.Table 1
MAEs between inverted OBCs and prescribed ones in TE1 (unit: cm).
TE1a
TE1b
a
b
a
b
West
1.5
1.7
1.7
1.8
South
1.1
1.2
1.3
1.4
East
0.8
2.6
2.0
2.7
North
1.9
1.2
2.0
2.6Table 2
RCs between inverted OBCs and prescribed ones in TE1.
TE1a
TE1b
a
b
a
b
West
0.94
0.99
0.80
0.97
South
0.99
0.79
0.96
0.72
East
0.98
0.96
0.86
0.88
North
0.97
0.86
0.96
0.70Prescribed OBCs and inverted ones in TE1.
(a)
(b)Cost function and gradient norm in TE1.
(a)
(b)


## 3.2. Inversion by Assimilating “Inaccurate Observations”

Considering that practice in situ observations may contain noises, artificial random errors are added to “observations.” In this section, based on the “observations” in TE1a, five TEs are carried out and the maximum percentages of error are 10%, 20%, 30%, 40%, and 50% in TE2a, TE2b, Te2c, TE2d and TE2e, respectively. Then these “observations” are assimilated to invert OBCs. Results are displayed in Figures5, 6, and 7.(a) Prescribed OBCs and inverted ones on the west open boundary in TE2. (b) Prescribed OBCs and inverted ones on the south open boundary in TE2. (c) Prescribed OBCs and inverted ones on the east open boundary in TE2. (d) Prescribed OBCs and inverted ones on the north open boundary in TE2.
(a)
(b)
(c)
(d)Cost function and gradient norm in TE2.
(a)
(b)MAEs and RCs in TE2 (aw, bw, as, bs, ae, be, an, and bn represent Fourier coefficients on west, south, east, and north open boundaries, respectively).
(a)
(b)Figure5 illustrates the inverse results. All the shapes of inverse OBCs are consistent with those of the prescribed ones. It is found that with increased artificial errors, the inverted OBCs become further away from the prescribed values, implying that larger errors will result in worse inversion of OBCs. The cost function and its gradient norm are plotted in Figure 6. After assimilation, the cost function is about 1% of its initial value in TE1a and about 3% in TE2e. Combined with MAEs and RCs between inverted OBCs and prescribed ones in Figure 7, where all MAEs are less than 3.5 cm and RCs are larger than 0.6, it is reasonable to believe that the prescribed OBCs have been inverted successfully though assimilating the “observations” is not accurate. But it cannot be denied that as artificial errors increased, the inversion results become worse with larger MAEs and smaller RCs. In addition, in these TEs, the inversion results of a are better than those of b.


## 4. Practical Experiments

The 3D internal tidal model is used to simulate theM2 internal tide around Hawaii by assimilating T/P data with OBCs optimized. IP strategy is used, which is the same as that used in TE1a. Both locations and values of observations are used in PE, which is different from TEs. In PE, the iterative assimilation terminates once the number of iteration steps reaches 50 exactly.After assimilation, the cost function is about 3% of its initial value. The simulated surface cochart is illustrated in Figure8. Considering the co-phase lines are disordered due to the influence of internal tide, only coamplitude lines are contoured. On the whole, the simulated cochart shows good agreement with that obtained from Oregon State University (OSU) tidal model in which only barotropic tide is simulated. The sea surface amplitude is about 25 cm at the southeast of Hawaii and smaller than 15 cm at the northwest. Careful inspection shows that the major difference between simulated results and OSU model results is in the short-scale fluctuations superposed on co-amplitude lines. These fluctuations can be considered as the surface manifestation modulated by the internal tide, from which the propagation of internal tide could be found. The M2 internal tide generated along the Hawaiian Ridge mainly propagates towards to northward and southeastward.Simulated co-chart (unit: cm, (a)) and that from OSU tidal model (b).
(a)
(b)Figure9 illustrates the difference of surface amplitude along tracks 006, 045, 082, 147, 184, and 223. Simulated results show good agreement with observations. Table 2 shows the corresponding MAEs. All MAEs of amplitude and phase are less than 2 cm and 8° (Table 3), respectively. Relative errors of amplitude are about 0.1, smaller than those in Ray and Mitchum [3], suggesting the advantage of this 3D internal tidal model involving adjoint method in simulating internal tides.Table 3
MAEs of amplitude and phase between observations and simulated results (values in brackets are mean relative errors of amplitude).
006
045
082
147
182
223
Amplitude (cm)
1.0 (0.07)
0.8 (0.04)
1.1 (0.06)
1.4 (0.11)
1.8 (0.14)
1.6 (0.10)
Phase (°)
4.3
3.0
4.7
6.0
7.7
5.9Observations and simulated results along tracks 006, 045, 082, 147, 184, and 223.
(a)
(b)Compared with the barotropic results from OSU tidal model, we can calculate the surface manifestation modulated by the internal tide, the maximum of which are 4.0, 4.2, 5.6, 3.6, 9.1, and 6.7 cm, corresponding to tracks 006, 045, 082, 147, 184, and 223, respectively. Except the manifestation of 9.1 cm along track 184, the values along other tracks are consistent with the estimation of Ray and Mitchum [3, 4].


## 5. Conclusion and Summary

In this paper, the optimization of OBCs in a 3D internal tidal model is investigated with the adjoint method. Fourier coefficients ofM2 internal tide on four open boundaries are inverted simultaneously, which means more variables are estimated in the inverse problem. Special step factor is used in the steepest descent method in order to minimize cost function efficiently without fluctuations. The reasonability and feasibility of the model and IP strategy are tested by TE1. In TE2, artificial random errors, of which the percentage is from 10% to 50%, are added to observations. After assimilation, the cost function reaches less than 3% of its initial values MAEs and RCs between inverted OBCs and prescribed ones are less than 3.5 cm and larger than 0.6, respectively, indicating that the prescribed OBCs are inverted successfully. But it cannot be denied that as artificial errors increased, the inversion results become worse gradually.In PE, the OBCs are optimized and theM2 internal tide around Hawaii is simulated by assimilating T/P data. The simulated co-chart shows good agreement with that from Oregon State University (OSU) tidal model and T/P observations on the whole. Careful inspection shows that the major difference between simulated results and OSU model results is in the short-scale fluctuations superposed on co-amplitude lines, which can be treated as the surface manifestation modulated by the internal tide. The maximum manifestation along tracks 006, 045, 082, 147, 184, and 223 are 4.0, 4.2, 5.6, 3.6, 9.1, and 6.7 cm, receptively. Except the manifestation of 9.1 cm along track 184, the values along other tracks are consistent with the estimation of Ray and Mitchum [3, 4].



---

*Source: 950926-2013-03-31.xml*