Adjoint Data Assimilation of the Flow on the Southern Flank of Georges Bank: March–June 1999

Abstract

An adjoint data assimilation method was incorporated into the ECOM-si coastal ocean circulation model and applied to assimilate the flow field on the southern flank of Georges Bank from March to June 1999. The model was driven by tidal forcing consisting of ten tidal constituents at the open boundary, observed winds, and surface heat fluxes. Numerical experiments were conducted following a strategy to understand the critical issues affecting the efficiency and accuracy of the assimilated flow field. The adjoint data assimilation method significantly improved the computational accuracy of subtidal currents, especially for the along-isobath velocity. The integration window length and iteration number were two parameters affecting the assimilation convergence rate toward the observations. In such a nonlinear dynamical system, using a window length close to the M₂ tidal period could make the adjoint model difficult to converge, no matter how many iterations were made. Reducing the time interval for the window length speeded up the convergence rate, but it was paid out with the sacrifice of statistical confidence. This assimilation experiment used a 6-h window length, which led to a faster convergence rate in the first ten iterations. The best-assimilated fields that satisfied the error criteria were obtained as the iteration number increased.

1. Introduction

One of the primary goals of developing a numerical model is to simulate and predict the temporospatial variability of temperature, salinity, and currents in the ocean. The four-dimensional (4-D) data assimilation is aimed at integrating observational data into a model for improving the simulation systematically. A data assimilation system, such as nudging, optimal interpolation, an adjoint model, and Kalman filters, consists of three components: an observational dataset, a hydrodynamical model, and an assimilation scheme. Since the data have measurement uncertainty errors and the model dynamics and parameters are imperfect, a well-constructed melding algorithm should enable to improve the match between observed and simulated variables within an error range of observations and modeling [1] Therefore, once a hydrodynamical model is developed, an appropriate assimilation algorithm becomes more critical in ensuring the success of its application to the ocean.

Optimal interpolation (IO), ensemble Kalman filtering (EnKF), and variance or adjoint methods are the most popular algorithms used in data assimilation. The IO and EnKF are based on the estimation theory [2]. They involve a sequence of estimation processes based on the minimization of the expected estimation errors in terms of the statistics for both the model and observations. The IO builds a linear combination of the observed and model-predicted variables by calculating their error covariances [3]. It is the most straightforward assimilation algorithm requiring a priori statistical assumption about the model noises and measurement uncertainties in observational data. The EnKF is the most sophisticated method, but ensemble simulations are computationally intensive to run [4,5,6,7,8]. The variance or adjoint method is developed following the control theory. In this algorithm, a cost function is defined by the difference between model-predicted and measured quantities. An optimal assimilation can be achieved by minimizing the cost function equations through a least-square approach [9,10,11,12,13]. This method is conceptually simple, less computationally intensive, and can be easily applied to linear and nonlinear models to improve model parameters, including initial and boundary conditions.

Various data assimilation algorithms have received much attention in the U.S. GLOBEC Northwest Atlantic/Georges Bank (GB) program. A real-time hindcast/forecast data assimilation finite-element circulation model (QUODDY-CASCO) was developed for the GB/Gulf of Maine (GOM) region [14,15,16,17,18]. Assimilations were conducted in hindcast simulation experiments during three GLOBEC cruises on GB, including 14–25 April, 4 May–8 June, and 10–31 May 1999. These assimilation experiments improved the model-simulated flow fields for the study of the ecosystem in the GB/GOM region.

A seven-current meter mooring array was deployed on the southern flank of GB in March 1999, with an across-bank layout centered at the 60 m isobath (the position of the tidally mixed front) (Figure 1). Instruments were in the water for about six months and recovered on 1 September 1999. All tripods were located between 55 and 70 isobaths, with a horizontal separation spacing of 2.0 to 3.0 km. Two settings were designed to cover spring and summer. The spring one started on 20 March and ended on 20 June and the summer one was from 20 June to 1 September. The spring was a period of no or intermittent weak vertical stratification, and the summer was a period of permanent vertical stratification except for the tidal-mixed region on the top of the bank. The scallop fishing in Closed Area II on GB was unanticipatedly opened in summer, 1999. The array was moved ~8 nm westward over that period. Except for tripod #4, all moorings in the array were recovered from the summer setting. The current sensors were acoustic doppler current meters (ADCPs), operating at 300 KHz, located on bottom-mounted tripods. The velocity was recorded with a 1.0-m vertical bin throughout the water column every 15 min. High-quality data were received in a depth range of 6–7 m from the sea surface and 3–4 m from the bottom. The measurement uncertainty was within 0.7 cm/s.

We incorporated an adjoint data assimilation method into the ECOM-si model. The ECOM-si model was the semi-implicit version of the Princeton Ocean Model (POM). This model was initially developed by [19] and modified by [20,21] to apply for the GB/GOM. The adjoint ECOM-si code was developed using the tangent linear and adjoint model compiler (TAMC) developed by [22]. The numerical experiments were conducted to assimilate the flow field over March–June 1999, with an objective to understand critically essential issues affecting the efficiency, accuracy, and convergence of this algorithm for the case with the inclusion of nonlinear advections in the inverse model.

This paper is organized as follows. Section 2 describes the adjoint assimilation method and the designs of numerical experiments. Section 3 presents forward model simulation results of tidal and subtidal currents. Section 4 shows assimilated total and subtidal currents. Section 5 discusses the issues raised from assimilation experiments. Finally, Section 6 summarizes the primary findings from this study.

2. The Adjoint Data Assimilation Model

ECOM-si used in this study is a three-dimensional, free-surface, σ-coordinated vertical transformation coastal ocean circulation model [23]. In ECOM-si, the Mellor and Yamada level 2.5 turbulent closure scheme was used for vertical viscosity parameterization [24]. Unlike the time-splitting solver in POM, ECOM-si is integrated using a semi-implicit scheme, which makes the Fortran code simpler by removing explicit time smoothing [19].

The adjoint ECOM-si model was conducted using a variational method. The primitive equations ECOM-si can be written in the vector form given as

$$\frac{ \partial \mathit{x}}{\partial t}=\mathit{F}\left(\mathit{x},\mathit{c}\right),$$

(1)

where x is a matrix array consisting of dependent variables such as the free-surface water level ($\zeta$), water temperature ($\theta$), water salinity (s) and the three-dimensional (3D) velocity components (u, v, and w). F is a nonlinear operator representing the ocean-governing equations, including the advective, Coriolis forcing, barotropic and baroclinic pressure gradient forcings, and diffusive terms in the momentum equations, as well as advective and diffusive terms for temperature and salinity equations. c is a matrix array containing the model parameters such as surface and bottom drag coefficients, light attenuation lengths, and open boundary and initial conditions. The model uses the Cartesian coordinates, in which x, y, and z are the cross-bank (northward: positive), along-bank (eastward: positive), and vertical axes of the orthogonal coordinate). The cost function in the adjoint model is defined as

$$\mathit{J}\left(\mathit{x},\mathit{c}\right)={{\displaystyle \int }}_{{\mathit{t}}_1}^{{\mathit{t}}_2}{{\displaystyle \int }}_{\Omega }^{}[\frac{\mathit{K}}{2}{\left(\mathit{x}-{\mathit{x}}_{\mathit{o}}\right)}^2+\frac{{\mathit{K}}_{\mathit{c}}}{2}{\left(\mathit{c}-{\mathit{c}}_{\mathit{o}}\right)}^2]d\mathsf{\Omega }dt,$$

(2)

where ${\mathit{x}}_{\mathit{o}}$ and ${\mathit{c}}_{\mathit{o}}$ are the observed vectors for dependent variables and model parameters, respectively; and K and K_c are validity coefficients, which were chosen as K = 0.5 for u and v, and K_c was taken as 1 in our numerical experiments. $\Omega$ is the numerical computational domain. t represents time, and (t₁,t₂) is the time assimilation window. The cost function measures a “distance” (error) between observed and simulated quantities. A Lagrange function was defined as

$$\mathit{L}(\mathit{x},\mathit{\lambda },\mathit{c})={\displaystyle \underset{{\mathit{t}}_1}{\overset{{\mathit{t}}_2}{\int }}{\displaystyle \underset{\Omega }{\int }\{\frac{\mathit{K}}{2}}}{(\mathit{x}-{\mathit{x}}_{\mathit{o}})}^2+\frac{{\mathit{K}}_{\mathit{c}}}{2}{(\mathit{c}-{\mathit{c}}_{\mathit{o}})}^2+{\mathit{\lambda }}^{\mathit{T}}[\frac{\partial \mathit{x}}{\partial \mathit{t}}-\mathit{F}(\mathit{x},\mathit{c}\left)\right]\}d\Omega dt,$$

(3)

where $\mathit{\lambda }$ is a Lagrange multiplier’s matrix array for x (also denoted as adjoint variables), and the superscript “T” is the sign for the transposed matrix. The variational method used in the adjoint assimilation model was designed to determine x with minimum L with respect to $\mathit{\lambda }$, x, and c in terms of least-square fitting, i.e.,

$$\delta \mathit{L}(\mathit{\lambda },\mathit{x},\mathit{c})=\frac{\partial \mathit{L}}{\partial \mathit{\lambda }}\partial \mathit{\lambda }+\frac{\partial \mathit{L}}{\partial \mathit{x}}\delta \mathit{x}+\frac{\partial \mathit{L}}{\partial \mathit{c}}\delta \mathit{c}=0.$$

(4)

This condition was equivalent to solving the Lagrange–Euler equations that satisfied the constraints in the form of

$$\frac{\partial \mathit{L}}{\partial \mathit{\lambda }}=0, \frac{\partial \mathit{L}}{\partial \mathit{x}}=0, \mathrm{and} \frac{\partial \mathit{L}}{\partial \mathit{c}}=0.$$

(5)

Three equations can be derived from these three conditions of restraint. They are the Lagrange–Euler equations from the original (forward) model equations of

$$\frac{\partial x}{\partial t}=F(t,x,c),$$

(6)

adjoint equations of

$$\left\{\begin{array}{c}\frac{\partial \mathit{\lambda }}{\partial t}=-\mathit{\lambda }\frac{\partial \mathit{F}(t,\mathit{x},\mathit{c})}{\partial \mathit{x}}+\mathit{K}(\mathit{x}-{\mathit{x}}_o)\hfill \\ {\left.\mathit{\lambda }\right|}_{\Sigma }=0\hfill \end{array}\right.,$$

(7)

where Σ is the boundary of the computational domain. The gradient of the cost function can be derived by $\partial \mathit{L}/\partial \mathit{c}=0$. That is,

$$\frac{\partial \mathit{F}}{\partial \mathit{c}}=\frac{{\mathit{K}}_{\mathit{c}}}{\mathit{\lambda }}(\mathit{c}-{\mathit{c}}_o).$$

(8)

Equation (8) can be rewritten as an iterative form as

$${\mathit{c}}^{n+1}={\mathit{c}}^n+{\left(\frac{\lambda }{{\mathit{K}}_{\mathit{c}}}\frac{\partial \mathit{F}}{\partial \mathit{c}}\right)}^n.$$

(9)

where ${\mathit{c}}^{\mathit{n}}$ is the previous nth iteration value of c, which is also equivalent to ${\mathit{c}}_{\mathit{o}}$. n indicates the iteration number.

The forward model of ECOM-si was converted to the adjoint model using the tangent linear adjoint model compiler (OpenAD) developed by [22]. The OpenAD is a numerical computer program that could automatically convert the differentiation equations to the adjoint code (https://www.mcs.anl.gov/OpenAD/, accessed on 24 November 2023). We used this software to automatically convert all subroutines in ECOM-si into adjoint codes. Then, the adjoint model was built by assembling all these subroutines based on the mathematical approach shown in Equations (6)–(9). A detailed description of OpenAD can be found at the website listed above.

A schematic of the adjoint assimilation approach is shown in Figure 2. First, the forward model of Equation (6) was run with measured wind forcing and surface heat flux in a tidal environment. Second, the model-simulated u and v within each integration window were output and compared with observed u and v at all available measurement positions. Then, the model-predicted and observed u and v differences were input into the adjoint model of Equation (7). Third, Equation (7) was run backward within the assimilation window, which provided a new set of $\lambda$. A new set of c (initial condition of u and v and the bottom drag coefficient) at n + 1 was then determined by solving Equation (9) with the input of $\lambda$ from the adjoint model. A criterion based on the difference ($\epsilon$) between ${\mathit{c}}^{\mathit{n}+1}$ and ${\mathit{c}}^{\mathit{n}}$ was used to judge if the model provided an optimal flow field of u and v. If not, the simulation returned to the forward model with a new set of c as its initial condition and all three procedures were repeated. The iteration was conducted until the criterion was satisfied. This approach assumed that the inaccurate simulation results were due to the improper model setup in parameterization [25]. In this study, the model parameters only considered the bottom stress coefficient since it was specified initially from an empirical formulation. However, in this case, a test showed that it required at least 40 iterations to meet the criterion of $\epsilon$ < 2.5 cm/s. For an assimilation window of 6 h, this iteration number produced the most rapid convergence of the assimilated currents to the observations. A detailed discussion on the accuracy of the assimilation regarding iteration number and assimilation window is given in Section 3.

Adjoint assimilation modeling was conducted for a two-dimensional (2-D) numerical experiment on the assumption that the change of all variables in the along-isobath direction was negligible. The experiments used the 2-D version of ECOM-si developed by [26]. This model was configured for a cross-bank section through the moored ADCP array, where the water depth varied from 100 to 35 m on the southern flank (Figure 1). [26] examined the dynamics attributed to tidal fronts over GB. They found that the tidal mixing front formation and along-frontal current were predominated by 2-D dynamical features, where the along-frontal barotropic and baroclinic pressure gradient could be ignored. The 2-D model, assuming no along-bank variations, was robust enough to simulate the stratified frontal currents observed on the south–north transect across the central bank (Figure 1). The tripod layout was on the cross-bank transect used in the 2-D experiments. This 2-D simplification made it easy to incorporate the observed data from the current meter array. It also helped avoid internal wave reflection at the northern open boundary.

The model used the σ-coordinate transformation defined by $\sigma =\left(z-\zeta \right)/\left(H+\zeta \right)$ where σ varied from 0 at the surface to −1 at the bottom, and H was the mean sea level depth. A total of 50 uniform σ-layers were specified in the vertical, with a resolution of 2 m at the 100 m isobath and 0.7 m at the 33.4 m isobath. The cross-bank domain was configured with a nonuniform horizontal grid. The horizontal resolution (∆x) was 471 m near and across the bank area and increased linearly to 10 km over 33 grid points outside the domain of interest. Twenty additional grid points with a resolution of 10 km were added at the northern end to extend the northern domain 200 km longer. A gravity wave radiation condition was specified for the surface elevation and currents at the northern open boundary. The internal model time step was 110.4 s, 405 time steps over an M₂ (12.42 h.) tidal period.

Under a stratified condition, a free-surface model can generate surface and internal gravity waves. The wavelengths of internal waves were much shorter than surface gravity waves. Since the radiation condition could not transfer all wave energies from the boundary, reducing the horizontal resolution near the northern open boundary could help filter short-wavelength interval waves [26].

The model was driven by the tidal forcing at the southern open boundary and calculated wind stress plus heat flux at the sea surface. The tidal forcing consisted of ten tidal elevation constituents, including M₂ (12.42 h), S₂ (12 h), N₂ (12.66 h), L₂ (12.19 h), K₂ (11.97 h), MU₂ (12.87 h), K₁ (23.93 h), P₁ (24.07 h), O₁ (25.82 h), and MSF (14.77 days). The amplitudes and phases of these tidal constituents were specified by the output from the global 0.5° × 0.5° tidal assimilation model developed by [27]. These values were adjusted to produce an optimal fit to tidal currents observed from the ADCP measurements.

At first, the model ran with tidal forcing only for three months after it was fully ramped up from the initial. After tidal currents in the model reached an equilibrium state, the calculated wind stress and surface heat flux were added at the correct tidal phase when ADCP data were recorded. The March-April climatologically hydrographic fields were used to set up the initial conditions of temperature and salinity (Figure 3). The climatological dataset was downloaded from the Dartmouth College Ocean Modeling Group website (http://www-nml.dartmouth.edu, accessed on 24 November 2023) [28,29]. The data were processed by investigators at the Bedford Institute of Oceanography. The wind stress was calculated using the meteorological data recorded on the NOAA environmental buoy #44011 east of GB (Figure 1 and Figure 4). Surface heat flux was calculated using an empirical formula based on the air-sea temperature difference recorded on buoy #44011. Relative humidity and cloud coverage were estimated using the data from the NOAA Climate Diagnostic Center (http://www.cdc.noaa.gov, accessed on 24 November 2023) [30,31]. The amplitudes and phases of the observed and simulated tidal currents for each tidal constituent were calculated using a MATLAB R2014 harmonic analysis program developed by [32].

To investigate processes in the coastal ocean using the adjoint method, the model was first run in the forward mode only for three months (from 19 March to 16 June 1999), then rerun with data assimilation. This strategy helped us first determine the level of accuracy for a 2-D model in predicting tidal and subtidal currents on GB and then evaluate if adding the adjoint data assimilation method could improve the model accuracy. It also helped identify some key issues in conducting a data assimilation experiment in a coastal ocean environment with strong tidal currents.

3. Forward Simulation Results for Tidal and Subtidal Currents

3.1. Tidal Simulation

The time series of model-simulated currents at each grid point were harmonically reconstructed into the ten tidal constituents: M₂, S₂, N₂, K₁, P₁, O₁, L₂, K₂, MU₂, and MSF. The same method was used for the observed currents recorded at the seven ADCP moorings. The model data comparisons were made for each tidal constituent and the time series of total tidal currents (the sum of ten tidal constituents). For example, station 567 was at the 62 m isobath. At this site, the vertical profile of the simulated M₂ tidal current was in good agreement with the observation (Figure 5), except for an underestimation in the minor axis in a ~25 m water column above the bottom and a slight overestimation in the major axis near a depth of ~35 m. Similar results were found in S₂, N₂, K₁, P₁, O₁, L₂, K₂, MU₂, and MSF, but their amplitudes were one order of magnitude smaller. The vertically averaged absolute values of the differences between model-simulated and observed currents were 0.68 and 1.05 cm/s in major and minor axes, 1.0° and 1.17° in orientation and phase for M₂. Statistics for the other nine tidal constituents are given in

Table 1. Difference between model-predicted and observed tidal current ellipse parameters.

	Depth(m)	$\Delta {\mathit{U}}_{\mathit{m}\mathit{a}\mathit{j}}$ (cm/s)	$\Delta {\mathit{U}}_{\mathbf{min}}$ (cm/s)	$\Delta \mathit{O}\mathit{r}\mathit{i}\mathit{e}\mathit{n}$ (Deg G)	$\Delta \mathit{P}\mathit{h}\mathit{a}\mathit{s}\mathit{e}$ (Deg G)
M2	10	0.5 ± 1.2	−1.0 ± 1.2	0 ± 3	2 ± 3
	20	1.2 ± 0.9	−0.1 ± 1.1	−1 ± 3	2 ± 3
	30	1.6 ± 1.0	0.0 ± 1.2	−1 ± 3	1 ± 2
	40	0.5 ± 0.8	−1.6 ± 0.9	−2 ± 2	1 ± 2
	50	0.2 ± 0.7	−1.9 ± 0.8	−1 ± 2	1 ± 2
	60	0.1 ± 0.8	−1.7 ± 0.8	−1 ± 2	0 ± 2
S2	10	0.9 ± 1.2	0.4 ± 1.2	−1 ± 26	8 ± 24
	20	0.2 ± 1.0	−0.1 ± 1.3	4 ± 23	0 ± 22
	30	0.2 ± 0.9	0.1 ± 1.1	6 ± 21	−2 ± 19
	40	0.3 ± 0.7	0.0 ± 0.9	9 ± 18	−2 ± 16
	50	0.4 ± 0.7	0.0 ± 0.8	10 ± 16	−5 ± 16
	60	0.1 ± 0.8	−0.2 ± 0.8	8 ± 18	−4 ± 19
N2	10	0.4 ± 1.2	−0.4 ± 1.3	−9 ± 13	6 ± 14
	20	0.0 ± 1.0	−0.8 ± 1.3	0 ± 14	−2 ± 12
	30	0.3 ± 0.9	−0.4 ± 1.1	2 ± 12	−2 ± 11
	40	0.6 ± 0.7	0.1 ± 0.9	0 ± 11	0 ± 10
	50	0.6 ± 0.7	0.0 ± 0.9	−2 ± 10	−1 ± 10
	60	0.2 ± 0.7	−0.4 ± 0.8	−1 ± 11	−3 ± 12
K1	10	0.8 ± 1.1	0.7 ± 0.9	15 ± 22	2 ± 26
	20	1.2 ± 0.9	0.8 ± 0.6	17 ± 20	2 ± 23
	30	0.8 ± 0.8	0.4 ± 0.5	18 ± 20	4 ± 24
	40	0.7 ± 0.6	0.6 ± 0.5	15 ± 14	8 ± 17
	50	0.9 ± 0.8	0.7 ± 0.5	13 ± 16	9 ± 17
	60	0.7 ± 0.6	0.5 ± 0.5	10 ± 17	13 ± 17
P1	10	0.3 ± 1.0	0.2 ± 0.9	15 ± 54	2 ± 59
	20	0.4 ± 0.7	0.2 ± 0.6	17 ± 58	2 ± 58
	30	0.3 ± 0.8	0.1 ± 0.5	18 ± 57	4 ± 69
	40	0.2 ± 0.6	0.2 ± 0.5	15 ± 40	7 ± 46
	50	0.3 ± 0.4	0.2 ± 0.4	13 ± 52	8 ± 52
	60	0.2 ± 0.5	0.2 ± 0.4	10 ± 55	13 ± 58
O1	10	0.1 ± 1.3	0.7 ± 1.1	0 ± 24	−14 ± 27
	20	0.6 ± 0.8	0.3 ± 0.7	−5 ± 27	−8 ± 33
	30	0.3 ± 0.8	0.2 ± 0.7	11 ± 21	−20 ± 25
	40	0.5 ± 0.7	0.4 ± 0.5	11 ± 34	−19 ± 38
	50	0.6 ± 0.6	0.4 ± 0.6	−13 ± 44	4 ± 43
	60	0.5 ± 0.5	0.3 ± 0.6	−20 ± 50	9 ± 49
L2	10	1.0 ± 1.1	1.3 ± 1.1	−27 ± 30	16 ± 31
	20	0.3 ± 1.0	0.7 ± 1.2	−14 ± 27	18 ± 27
	30	0.9 ± 0.8	0.9 ± 0.9	1 ± 30	6 ± 30
	40	0.8 ± 0.7	0.6 ± 1.0	1 ± 23	8 ± 23
	50	0.7 ± 0.8	0.4 ± 0.8	−2 ± 26	18 ± 24
	60	1.0 ± 0.6	0.7 ± 0.8	−2 ± 34	23 ± 33
K2	10	0.2 ± 1.4	0.1 ± 1.4	−1 ± 95	8 ± 80
	20	0.0 ± 1.2	0.0 ± 1.4	4 ± 87	0 ± 68
	30	0.1 ± 1.1	0.0 ± 1.1	6 ± 91	−2 ± 82
	40	0.1 ± 0.9	0.0 ± 1.1	9 ± 72	−2 ± 63
	50	0.1 ± 0.8	0.0 ± 0.9	10 ± 72	−5 ± 65
	60	0.0 ± 0.8	−0.1 ± 0.8	8 ± 92	−3 ± 80
MU2	10	0.3 ± 1.1	0.9 ± 1.3	−31 ± 45	3 ± 45
	20	0.5 ± 0.9	0.4 ± 1.0	−28 ± 74	12 ± 68
	30	0.3 ± 0.9	0.1 ± 0.8	−17 ± 93	−8 ± 89
	40	0.2 ± 0.6	0.0 ± 0.8	11 ± 95	−43 ± 82
	50	0.1 ± 0.7	−0.1 ± 0.7	−9 ± 87	−32 ± 83
	60	0.1 ± 0.7	−0.1 ± 0.7	−17 ± 94	−35 ± 88
MSF	10	0.9 ± 2.5	0.1 ± 2.5	−58 ± 105	−94 ± 122
	20	0.1 ± 2.2	0.0 ± 2.2	−30 ± 79	−76 ± 89
	30	0.7 ± 2.3	−0.3 ± 1.5	−21 ± 41	−78 ± 176
	40	0.9 ± 2.3	−0.1 ± 0.9	−7 ± 31	−57 ± 199
	50	0.7 ± 1.9	0.1 ± 0.5	13 ± 21	−43 ± 137
	60	0.3 ± 1.4	0.1 ± 0.6	17 ± 26	−41 ± 113
Overal deviation		0.6	0.6	14	24

Note: $\Delta U_{maj}$, $\Delta U_{\min }$, $\Delta Orien$, and $\Delta Phase$ are differences between model-predicted and observed major axis, minor axis, orientation, and phase of the tidal current ellipse at site 567, respectively.

. These differences were within the estimated uncertainty of the ADCP measurements at this site. However, the overall simulation errors in the along- and cross-bank tidal velocity components over the seven ADCP sites were 3.81 cm/s and 2.33 cm/s, respectively. These errors are about 1.0–1.5 cm/s over the estimated overall uncertainty of the ADCP current measurements.

3.2. Simulation of Subtidal Currents

A priori requirement to ensure success in data assimilation is that the model must capture the right physics. This means that the forward model (without data assimilation) should be robust enough to reproduce the significant tidal and subtidal features observed in the natural ocean. For this reason, in addition to tidal forcing, the wind stress and heat flux were also considered in numerical experiments over 19 March–16 June 1999. The simulated total (tidal plus subtidal) and subtidal flows were compared with observed currents at measurement depths of seven moorings.

The forward model reasonably simulated the total velocity measured at the ADCP sites. An example is given in Figure 6, in which the model-predicted and observed along- and cross-bank total velocities approximately agreed with each other on both semidiurnal (short-term) and weekly-to-monthly (long-term) time scales. However, the model overestimated the flow magnitude, especially in the upper 40 m, where wind forcing was more critical. The simulated and observed velocity difference was much larger than the tidal-induced errors, implying that such a difference was partially due to the inaccurate wind stress and heat flux forcing conditions. It was supported by a comparison between simulated and observed subtidal currents shown in Figure 7. The model significantly overestimated the low-frequency variation of the along-bank velocity, especially during the strong southward wind period in April. In addition, the model also underestimated the wind-induced, low-frequency variation of the cross-bank velocity near the bottom. The observed velocity showed the time-coherent variation near the bottom and surface, but the model did not capture these features.

Despite these deficiencies, the forward model did capture the primary pattern of the low-frequency currents in response to the varying wind stress at ADCP measurement sites. It suggested that the model was robust enough to resolve the lowest order of the correct physics (tidal plus surface forcing) controlling current variability on the southern flank of GB during the simulation period. Meanwhile, the 2-D model did not resolve the basin-scale pressure field set up by the larger-scale wind stress [34]. Although errors due to the basin-scale pressure gradient forcing and the wind stress spatial variability could cause the subtidal model-data differences, they were secondary to surface and boundary forcing conditions on the cross-bank moored ADCP transect.

4. Assimilation Results for Tidal and Subtidal Currents

Turning on adjoint assimilation significantly improved the simulation accuracy of the total and subtidal currents at the ADCP measurement sites. It was evident by comparing Figure 8 with Figure 6. This comparison showed that the assimilated along-bank velocity quickly converged to the observed velocity after 12 iterations within a 6 h integration window, even though no significant improvement was found for the cross-bank velocity. Without data assimilation, the overall model errors in total along- and cross-bank velocities were 7.56 and 9.91 cm/s at the ADCP measurement points. After the adjoint assimilation was applied, the overall uncertainties were reduced to 7.10 and 6.63 cm/s (

Table 2. The standard deviation of the difference between simulated/assimilated and observed along-bank and cross-bank currents.

Cases	Simulation		Assimilation
	$\Delta \mathit{u}$ (cm/s)	$\Delta \mathit{v}$ (cm/s)	$\Delta \mathit{u}$ (cm/s)	$\Delta \mathit{v}$ (cm/s)
Total current error	7.56	9.91	7.10	6.63
Subtidal current error	3.75	7.58	2.77	4.16

Note: $\Delta u$ and $\Delta v$  are the standard deviation of the along- and cross-bank current differences, respectively.

). A significant improvement was also evident for the subtidal current, shown in Figure 9. In the case without assimilations, overall simulation errors in along- and cross-bank subtidal velocities were 3.75 and 7.58 cm/s, respectively, at the ADCP measurement points. These errors were reduced to 2.77 and 4.16 cm/s in the assimilation case. The overall improvement at these comparison sites was 26 and 45% in the cross-and along-bank directions.

The assimilated subtidal velocity exhibited the same spatial distribution patterns as observations at seven ADCP sites. In Figure 10, for example, the data showed a subsurface along-bank velocity maximum at a depth of 20–30 m on 11 June 1999. This feature was reproduced in the model solution with assimilation. A relatively strong vertical shear of the cross-bank subtidal velocity was found in the upper 20 m on the same day, which also appeared in the vertical distribution of the assimilated subtidal flow.

5. Discussion

Our simple model experiments show that the adjoint assimilation method can minimize model data differences in a least-square sense for current assimilation. However, two critical issues arise regarding this method: (1) computational efficiency and (2) convergence. We had difficulty using the $\mathsf{\epsilon }$-criterion to control the number of iterations required for a long-term (three-month) assimilation since the time required to reach a specific criterion was too long. Therefore, based on the time required to run an ocean forecast model, the adjoint data assimilation approach used in these experiments seems inefficient and practical enough to be used for the forecast application in the coastal ocean without a high-performance MPI cluster computer.

Convergence is sensitive to the assimilation window and iterations. Unlike the open ocean, many coastal regions are characterized by dominant semidiurnal tidal motion. Previous current measurements on GB revealed that the independent time scale of the tidal flow is about one hour [35]. This means that a semidiurnal tidal cycle only contains about 12 independent data records. Using a short time integration window in an adjoint model could reach a faster convergence rate toward the observation, but it is accompanied by a sacrifice of statistical confidence level. Determinations of a proper integration window with a rapid convergence rate and sufficient statistical confidence are the critical issues that should be considered when the adjoint data assimilation method is applied to the tidal-dominant coastal oceans.

Our numerical experiments showed that the error convergence rate in the assimilated velocity was sensitive to the forward/backward integration windows. When the window length was close to the M₂ tidal period, the adjoint model failed to reach a convergence solution. When a 12-h window was chosen, the standard deviation error between assimilated and observed along-bank velocities quickly decreased in the first three iterations, but gradually increased after the 6th iteration (Figure 11). As the window length was reduced to 6 h (about half a M₂ tidal period), the convergence was significantly improved. The standard deviation error exhibited a rapid decrease rate in the first ten iterations and gradually converged to the observation as iterations increased (Figure 11). The deviation error was reduced to 2.0 cm/s after 40 iterations.

In our experiments, no attempts have been made to assimilate the temperature and salinity because sufficient hydrographic data were unavailable across the entire study area during the ADCP measurement period. It is still being determined whether the convergence rate can change as more control parameters are added to this M₂ tidal-dominant coastal region.

6. Summary

The ECOM-si adjoint model was developed and applied to assimilate the March–June 1999 flow field across the southern flank of GB. The current measurement data in the assimilation were from seven ADCP tripod sites centered at the 60 m isobath. Numerical experiments were carried out to understand critical issues affecting the assimilation efficiency and accuracy. Without assimilation, the model was robust enough to reproduce the temporospatial variability of tidal and subtidal flows observed at seven ADCP measurement sites. However, the tidal currents were underestimated near the bottom and overestimated near the surface, with an overall error of 2–4 cm/s. Adding wind stress and surface heat flux in the model captured the variability of the low-frequency flow, but the subtidal flow was significantly overpredicted during relatively strong wind periods. The maximum simulation error was up to ~20 cm/s. The adjoint assimilation significantly improved tidal and subtidal current predictions, especially during the storm periods in April and May. The assimilation accounted for a 45% improvement for the along-bank current and 26% for the cross-bank current.

The convergence rate of the assimilated currents toward observed currents depends on the assimilation window and iteration number. When the window length was close to the M₂ tidal period, the assimilated flow could not converge to the observed flow as the iteration number increased. Shortening the window length provided a faster convergence rate, but the statistical confidence level was reduced. As a practical compromise, a 6-+h window was recommended for the semidiurnal tidal dominant region. In our experiments, this window length led to a rapid convergence rate in the first ten iterations. After 40 iterations, the model reached a numerical solution that satisfied a given $\mathsf{\epsilon }$-criterion.

Our assimilation experiments focused only on using observed currents. Since no hydrographic measurements were made near the ADCP array during the study period, no efforts were made to assimilate temperature and salinity data.