Picard iterations for a finite element shallow water equation model

Ocean Modelling 10 (2005) 316–341 www.elsevier.com/locate/ocemod

Picard iterations for a finite element shallow water equation model Julia C. Muccino *, Hao Luo Department of Civil and Environmental Engineering, Arizona State University, ECG 252, P.O. Box 875306, Tempe, AZ 85287-5306, United States Received 22 July 2004; received in revised form 16 September 2004; accepted 20 September 2004 Available online 2 November 2004

Abstract ADCIRC, a finite element circulation model for shelves, coasts and estuaries, will be used for variational data assimilation. The nonlinear Euler–Lagrange (EL) problem will be solved using the iterated indirect representer algorithm. This algorithm makes such large, nonlinear but functionally smooth optimization problems feasible by iterating on linear approximations of the nonlinear problem (Picard iterations) and by making preconditioned searches in the ‘‘data subspace’’ at each iterate. Before solving the nonlinear EL using such Picard iterations, it essential that the iteration scheme be carefully examined within the framework of the nonassimilative or forward problem. The purpose of this paper is (1) to detail a Picard iteration procedure for ADCIRC, including the problematic bottom friction term; (2) to examine the ability of the iteration scheme to recover the nonlinear forward solution from deficient background fields; and (3) to present a study of different interpolation methods for reducing the memory/disk requirements of the iteration scheme. The iteration scheme is shown to be quite robust in its ability to recover the nonlinear solution from a variety of deficient background fields. A new cubic Hermitian interpolation method is shown to be a more effective alternative to standard linear interpolation for reducing memory/disk requirements, especially for high frequency overtides.
2004 Elsevier Ltd. All rights reserved.

*

Corresponding author. Tel.: +1 480 965 0598. E-mail addresses: [email protected] (J.C. Muccino), [email protected] (H. Luo).

1463-5003/$ - see front matter
2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ocemod.2004.09.006

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

317

1. Introduction Data assimilation serves several purposes. The ostensible purpose is preparing analyses of ocean circulation (e.g., gridded fields of elevation and velocity) based on observations and model dynamics. These analyses aid the scientific study of the dynamics of the ocean, and may also be used to initialize the operational ocean forecasts which are increasingly needed for marine resource management, for search and rescue, for weather and climate prediction and for national defense. Other assimilation products include: analyses of errors in prior information such as internal parameterizations (e.g., Bennett et al., 1998); covariances of errors in the analyses of circulation and forcing fields (e.g., Bennett and Thorburn, 1992); quantitative assessments of the efficiency of the observing system that collected the data (e.g., McIntosh, 1987); and hypothesis tests (e.g., Muccino et al., 2004). All of these products can be generated by fixed-interval smoothing of nonlinear dynamics and nonlinear observing processes using the iterated indirect representer algorithm. The dynamics may be weak constraints; as a result, the number of computational degrees of freedom may be very large (e.g., Bennett et al., 1996, 1997; Xu et al., submitted). The representer algorithm makes such large, nonlinear but functionally smooth optimization problems feasible by iterating on linear approximations of the nonlinear problem using linear Picard iterations (e.g., Bennett and Thorburn, 1992; Chua and Bennett, 2001; Bennett, 2002) and by making preconditioned searches in the ‘‘data subspace’’ at each iterate. (e.g., Egbert et al., 1994; Chua and Bennett, 2001; Bennett, 2002). This algorithm is quite intricate and complex, deterring many scientists and engineers from implementing it. However, much of the algorithm is model independent and a new modular implementation requires only that the user supply a linear iteration scheme along with the forward and adjoint codes associated with a particular ocean model. This Inverse Ocean Modeling system (IOM) uses functional programming to custom generate the remainder of the iterated indirect representer algorithm for each ocean model, including a preconditioned conjugate gradient solver, efficient convolutional algorithms, diagnostic tools, hypothesis tests and measurement operators. 1 ADCIRC (Luettich et al., 1991; Westerink et al., 1993; Kolar et al., 1994; Luettich and Westerink, 2000) is one of several models that will interface with the IOM. A finite element circulation model for shelves, coasts and estuaries, ADCIRC is widely used by the university research community, the US Geological Survey and the Army Corps of Engineers as well as for matters of national security by the Naval Research Laboratory and the US Naval Oceanographic Office. Typical applications include modeling tidally and wind driven circulation in coastal waters, forecasting hurricane storm surge and flooding, dredging feasibility and material disposal studies and larval transport studies (Luettich and Westerink, 2000). The purpose of this paper is to detail a linear Picard iteration scheme for ADCIRC, including the problematic bottom friction term. An extensive investigation of the ability of the iteration

1

The measurement operators are generated by the IOM, which composes the specific instrument kernel (supplied by the observer) with a routine that interpolates from model coordinates to observational coordinates (supplied by the modeler).

318

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

scheme to recover the nonlinear solution from a variety of deficient background fields is also presented. Finally, owing to the fact that this and other iteration schemes are very expensive in terms of both memory and disk requirements, two approaches to reducing this requirement are evaluated, including a new Hermitian cubic interpolation method. 2. ADCIRC: the governing equations The governing equations of ADCIRC (in two-dimensional Cartesian coordinates) are the ‘‘primitive’’ shallow water equations; specifically these are the continuity equation og þ r
ðVH Þ ¼ 0 ot and the momentum equation (in conservative form) Lðg; VÞ ¼

ð1Þ


 oðH VÞ ps ss sb þ gðg  anÞ  Er2 ðH VÞ  þ ¼ 0 M ðg; VÞ ¼ þ r
ðH VVÞ þ H f  V þ H r ot q0 q0 q0 C

ð2Þ where g(x, y, t) is the free surface displacement h(x, y), V(x, y, t) is the depth-averaged velocity vector with scalar Cartesian components u(x, y, t) and v(x, y, t), H(x, y, t) = h + g is the total depth, ^ is the Coriolis parameter, X is the angular velocity of the earth, / is the latitude, f ¼ 2X sin /k g is the acceleration of gravity, ps(x, y, t) is the atmospheric pressure at the free surface, q0 is the reference density of water, a is the effective Earth elasticity factor, n(x, y, t) is the Newtonian equilibrium tidal potential, E is the eddy coefficient for horizontal diffusion, ss(x, y, t) is the wind stress vector applied at the water surface sb(x, y, t) and is the bottom stress vector. Appropriate boundary conditions (either no-flow, specified elevation, or some combination) and initial conditions must also be supplied. Solutions of the shallow water equations, (1) and (2), using finite element methods with linear triangular elements are characterized by artificial, short wavelength spatial oscillations (e.g., Kinnmark, 1986). These oscillations are suppressed in a physically correct way by implementation of the second order generalized wave continuity equation formulation (Lynch and Gray, 1979; Kinnmark, 1986) oL  r
MC þ s0 L ¼ 0 ð3Þ ot where s0 is a numerical parameter reflecting the relative balance in the generalized wave equation between the primitive continuity equation and the wave continuity equation. ADCIRC solves (3) with the advective terms in nonconservative form W ðg; VÞ ¼

W ðg; VÞ ¼

o2 g og þ s0  VH
rs0 þ r 2 ot ot 
  og ps og ss sb þ gðg  anÞ  Er  þ þ s0 VH
V  H V
rV  H f  V  H r ot ot q0 q0 q0

¼0

ð4Þ

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

319

It has been shown that local mass conservation is improved by consistent treatment of the advective terms in the wave and momentum equations. Thus, in conjunction with (4), ADCIRC solves the nonconservative momentum equation Mðg; VÞ ¼ ðMC  VLÞ=H


 oV ps E ss sb þ V
rV þ f  V þ r þ ¼0 ¼ þ gðg  anÞ  r2 V  ot H q0 q0 H q0 H

ð5Þ

3. An iteration scheme for ADCIRC The iterated indirect representer algorithm assumes nonlinear dynamics and nonlinear observing systems and thus solves a nonlinear Euler–Lagrange (EL) problem (e.g., Bennett and Thorburn, 1992; Bennett et al., 1996; Ngodock et al., 2000). It uses Picard iterations to reduce the nonlinear EL problem to a sequence of linear EL problems, each of which is solved indirectly by a data space search. One of many approaches to generating an appropriate iteration scheme is Tangent Linearization (TL) (e.g., Le Dimet and Talagrand, 1986; Le Dimet et al., 1997; Ngodock et al., 2000; Uboldi and Kamachi, 2000); the nonlinear terms are expanded in Taylor series about a background field, and all second order and higher terms are neglected. For example, consider the local advection term in the generalized wave continuity equation (4)  n og ogn1 ogn1 oðgn  gn1 Þ 2 þ Oðgn  gn1 Þ þ ðVn  Vn1 Þ þ Vn1 ¼ Vn1 V ot ot ot ot ogn1 ogn ogn1 þ Vn1  Vn1 ð6Þ ot ot ot where the superscript n represents the current iteration and the superscript n  1 represents a background field. All nonlinear terms in the ADCIRC equation (3) and are linearized using the TL scheme, with the notable exception of the bottom friction terms in both equations; these terms require special consideration. Consider first the bottom friction term in the generalized wave continuity equation using a quadratic parameterization  n sb ¼ cf An Vn ð7Þ q0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 where An ¼ ðun Þ þ ðvn Þ and cf is a dimensionless friction coefficient. The TL of (7) is  n sb cf   cf An1 Vn þ n1 un1 ðun  un1 Þ þ vn1 ðvn  vn1 Þ Vn1 ð8Þ q0 A  Vn

This iteration scheme is problematic because, as An1 approaches zero, the second term on the RHS of (8) becomes undefined. Of course, other iteration schemes are possible. For example, one might choose to simply neglect the problematic second term on the RHS of (8); unfortunately the sequence of solutions thus obtained does not converge (i.e., the difference between successive

320

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

iterations increases with time until eventually computational overflow errors are encountered). As an alternative approach, the original expression for the bottom friction is slightly reformulated  n sb ¼ cf Bn Vn ð9Þ q0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 where Bn ¼ ðun Þ þ ðvn Þ þ b2 and b is a small positive constant. Note that although Bn will be nonzero for a resting fluid, the entire bottom friction term will be zero in that case, regardless of the value of Bn, since Vn = 0; thus the reformulated friction term does not induce motion in a resting fluid. The absolute value of the error introduced by the modification in (9) is limited by approximately E ¼ cf ðBn  An Þ j Vjn < cf 0:5b2

ð10Þ

Throughout this paper, b = 0.0001 m/s; for typical values of cf, which range from 0.001 to 0.01, E is substantially less than roundoff error. The TL of the reformulated bottom friction term (9) is  n sb cf   cf Bn1 Vn þ n1 un1 ðun  un1 Þ þ vn1 ðvn  vn1 Þ Vn1 ð11Þ q0 B where now the denominator of the second term on the RHS has a minimum value of b (rather than zero) and is therefore always defined. Similarly, the quadratic formulation of the bottom friction term in the nonconservative momentum equation is reformulated similarly, and the TL of that term is used in the iteration scheme. The fully linearized equations are provided in the Appendix A. These linear equations are formally consistent with the nonlinear equations (4) and (5), with sb/q0 defined as in (9), as n ! 1. Thus, if the solution of the nonlinear equations is unique and if the sequence of linearized solutions converges, the converged solution must be the solution of the nonlinear equations. Oliger and Sundstrom (1978) show that if the nonlinear shallow water equations (like (3) and (5) modulo rotation, bottom friction and diffusion terms) have at least one smooth solution for a given geometry and forcing, then it can have only one such solution (see Bennett, 2002). It is elementary to extend their proof to include rotation and bottom friction, although, to our knowledge, inclusion of conventional diffusion has yet to be accomplished. However, none of the experiments presented here include diffusion. 2 Thus, these solutions to the nonlinear shallow water equations are unique and, if the sequence of linearized solutions converges, it must converge to the solution of the original nonlinear equations and not to some other field. In some cases, particularly if the background field is sufficiently poor, the linearized solutions may not converge. In our experience with linearized ADCIRC, lack of convergence is manifest by differences between subsequent linearized solutions that increase with time until eventually computational overflow errors are encountered. Other researchers (e.g., Bennett and Thorburn, 1992;

2

In practice, diffusion is often not needed in ADCIRC simulations, e.g., in the western North Atlantic and Gulf of Mexico (Kolar et al., 1994) and the Neuse River estuary (Luettich et al., 2002). However, sometimes modest diffusion is necessary for stability, e.g., in the Beaufort Inlet (Hench and Luettich, 2003) and Cape Lookout Shole (McNinch and Luettich, 2000).

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

321

Bennett et al., 1993) have encountered the situation where the linearized solutions do not converge, but rather cycle around apparent limit points every three or four iterations. We have not observed such cycling in our linearized ADCIRC solutions.

4. Dataless iterations The iteration scheme described here is based on the TL operator, which is schematically defined as follows. Let D½a ¼ F

ð12Þ

where D is a nonlinear operator that operates on a and F is the forcing. Also D½a þ Da ¼ F þ Df

ð13Þ

where Df is a perturbation to the forcing. The nonlinear perturbation operator would be: D0nl ½Da; a ¼ D½a þ Da  D½a

ð14Þ

where D0nl operates on Da and a is now a parameter. A linear perturbation operator, D0l , would be D0nl modulo its nonlinear terms; indeed this operator is identical to the TL operator just described. There are any number of approaches to testing the iteration scheme which is not the same as testing the boundedness of the TL operator (D0l ). In the latter case, one might let: ðD0nl Þ1 ½Df ; a ¼ Danl

and ðD0l Þ1 ½Df ; a ¼ Dal

ð15Þ

where these inverse operators operate on Df, with a as a parameter, yielding Danl and Dal, respectively. If the nonlinear perturbation solution (Danl) and the linear perturbation solution (Dal) are of the same order, then the two operators have roughly the same sensitivity. This general numerical analysis gives insight into the stability and accuracy of the TL operator. Indeed, although we do use the TL operator here, any number of other linear operators could be used to accomplish the linear Picard iterations. Our objective here is to test the specific Picard iteration scheme, including the linear operator and its nonhomogeneous terms (e.g., the last term on the second line of (6)). An effective way to accomplish this is the so-called ‘‘dataless iterations.’’ It is important to note that the iteration scheme described in Section 3 has not been developed to solve the nonlinear governing equations; their solution can be found more efficiently by solving them directly (e.g., ADCIRC). Rather, the iteration scheme has been developed to solve the nonlinear EL equations, of which the forward model is a member. Moreover, generation of the sequences of EL equations without data (i.e., the forward model) has been found to be an invaluable preparatory test that gives considerable insight into the convergence behavior of the algorithm in a data assimilation context (Bennett and Thorburn, 1992; Ngodock et al., 2000; Muccino and Bennett, 2002), although such iterations are not a comprehensive mathematical test of the stability and accuracy of the linearized operator. In fact, dataless iterations are actually a test of the ability of the iteration scheme to apply incremental perturbations at each iterate such that, ultimately, the nonlinear solution is recovered. Convergence of dataless iterations is considered to be a practical prerequisite for convergence of the iteration scheme within a data assimilation framework.

322

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

Dataless iterations proceed with the following steps: 1. First, a background field (n = 0) is specified for all space and time to begin the iterations. The background field may be, for example, a resting field (i.e., all dependent variables are zero), the solution to a simplified model, or a field inferred from climatological data; 2. Then, using the linearized equations, a solution is generated over all space and time for n = 1, noting the pointwise difference between this solution and the background field. A convergence criterion is applied, such as j g1 ðx; y; tÞ  g0 ðx; y; tÞ j< eg

ð16aÞ

j u1 ðx; y; tÞ  u0 ðx; y; tÞ j< eu

ð16bÞ

j v1 ðx; y; tÞ  v0 ðx; y; tÞ j< ev

ð16cÞ

for all (x, y, t) where eg,eu and ev are specified tolerances; if the convergence criterion is satisfied, the solution is deemed converged. If the pointwise differences increase with time and computational overflow errors are encountered, the iterations have failed to converge. In either case, the iterations stop. If the convergence criterion is not satisfied and no overflow errors are encountered, n is increased to 2 and iterations continue with Step 3; 3. This step is much like Step 2, except that the solution from the previous iteration is used as the background field. The nth linearized solution is computed and the convergence criterion is applied j gn ðx; y; tÞ  gn1 ðx; y; tÞ j< eg

ð17aÞ

j un ðx; y; tÞ  un1 ðx; y; tÞ j< eu

ð17bÞ

j vn ðx; y; tÞ  vn1 ðx; y; tÞ j< ev

ð17cÞ

in the same way; again, if the convergence criterion is satisfied or if overflow errors are encountered, the iterations stop. If the convergence criterion is not satisfied and no overflow errors are encountered, n is increased by one and iterations continue by repeating Step 3. Of course, in these dataless iterations, one could replace (17) with a convergence criteria that compares the nth iteration with the known nonlinear solution. However, within the more general data assimilation framework, the optimal solution to which the Picard iterations converge is not known. In preparation for the latter, we use the criterion in (17) for these dataless iterations. Examples of a sequence of solutions that does converge, as well as an example of a sequence of solutions that does not converge are shown in Section 5. Representative results from two highly nonlinear domains are presented in the following two sections.

5. Bight of

ABACO

(bottom friction terms dominate)

The nonlinear bottom friction term is0 the most problematic term in the iteration scheme. The physics of the Bight of Abaco (2630 N, 77W) is dominated by this term owing to the very

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

323

Fig. 1. Bathymetry (m) of the Bight of Abaco. The depth at Node 1 is 3.0 m and the depth at Node 2 is 6.75 m. Term balances at these nodes are shown in Fig. 2.

shallow bathymetric depth (h is everywhere less than 10 m, as shown in Fig. 1), and so it is used as a relatively stringent test of the iteration scheme. The domain is bounded by land along the southern, eastern and northern boundaries by the islands of Great Abaco, Little Abaco and Grand Bahama. The flow is driven by O1, K1, N2, M2, and S2 tidal constituents along the western boundary with amplitudes (phases) of 0.075 m (194.8), 0.095 m (2.06.3), 0.100 m (340.0), 0.395 m (0) and 0.060 m (43.0), respectively. The simulation is spun up from rest. For five simulation days, the simulation is ‘‘ramped’’; i.e., the forcing is multiplied by a hyperbolic tangent function. Then the simulation is run using the full forcing for another 5 days, for a total simulation length of 10 days. The time step is 219 s. The finite element grid (not shown) has 926 nodes and 1696 elements with a relatively constant node spacing of approximately 2.7 km. There is no horizontal eddy viscosity, the nonlinear bottom friction coefficient, cf, is 0.009, and the Coriolis parameter, f, is ^ 5.9 · 105 s1 k. Time series of term balances in the nonlinear solution at two locations within the domain are presented in Fig. 2. At Node 1 (see Fig. 1 for node location) where the water depth is very shallow (3.0 m), the surface gradient term and the bottom friction term roughly balance one another in both the wave equation and the momentum equation, with other terms playing relatively minor roles. Even in deeper regions, as represented by Node 2 (depth of 6.75), the bottom friction is still a substantial contributor to the term balance with a magnitude of approximately 40% or more of the leading surface gradient and time derivative terms. The iteration scheme requires a background (n  1) field for each iteration. For the first iteration, when n = 1, a background field must be specified by the modeler; for all subsequent iterations, the background field is the solution from the previous iteration. Hereafter, ‘‘background field’’ refers specifically to the background field that must be specified for the first iteration. The background field used here is the solution to the linear generalized wave equation

324

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

Scaled term magnitude for wave equation

Node 1 (depth = 3.0 m)

Node 2 (depth = 6.75 m)

1

1

1

0.8

0.8

0.6

0.6

0.4

0.4

1

0.2

0.2

0

0

0

0.2

0.2

0.4

0.4

0.6

0.6

0.8

0.8

-1 1

1 7.5

0

8

8.5

9

9.5

10

Scaled term magnitude for x-momentum equation

1

1 7.5

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

0

0.2

0.2

0.4

0.4

0.6

0.6

8.5

9

9.5

10

-1 1

1 7.5

9.5

10

8

8.5

9

9.5

10

8.5

9

9

9.5

10

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

0.2

0.2

0.4

0.4

0.6

0.6

0

0

0.8

1 7.5

9

0.8

8

1

-1

8.5

0

0.8

-1 1

8

1

0.8

1 7.5

Scaled term magnitude for y-momentum equation

-1 1

0.8

8

8

8.5

9

9 Time (days)

9.5

10

10

Surface gradient term

-1

1 7.5

8

8

10

Time (days) Finite amplitude term

Bottom friction term

τ 0 terms

Time derivative term

Coriolis term

Fig. 2. Scaled term balances for the Bight of Abaco at a shallow Node 1 (left) and deeper Node 2 (right). See Fig. 1 for node locations. Only the relative magnitude of terms within a panel are significant. Note the time derivative in the wave equation is o2g/ot2 and the time derivative in the momentum equations is oV/ot. For clarity, within each figure, terms with magnitudes less than 10% of the largest term in that figure are not shown.

o2 g og þ s0  hV
rs0  r 2 ot  ot
  ps og ss sb þ gðg  anÞ  Er  þ þ s0 Vh ¼ 0
hf  V  hr ot q0 q0 q0

ð18Þ

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

325

and the linear momentum equation
 oV p E ss sb þ f  V þ r s þ gðg  anÞ  r2 V  þ ¼0 ot h q0 q0 h q0 h

ð19Þ

and where the bottom friction terms are approximated using the linear relationships sb sb ¼ cfl hVn ¼ cfl Vn q0 q0 h

ð20Þ

where cfl is a linear bottom friction coefficient. Note, only the background field is computed using these linear equations. The Picard iterations use the linearized equations given in Appendix A. The magnitude of these linear bottom friction terms should be consistent (at least within an order of magnitude) with the magnitude of the analogous quadratic bottom friction terms. Thus, in general, cfl 5 cf, but rather Bn ð21Þ h Here, based on the velocity and depth scales of 0.5 m/s and 10 m, respectively, cfl = 0.05cf. The convergence of the elevation solution is swift, as shown in Fig. 3; this convergence behavior is representative of the convergence of elevation and velocity solutions throughout the domain. As mentioned earlier, if the background field is sufficiently poor, the iteration scheme is unable to recover the nonlinear solution; the sequence of linear solutions does not converge. One example cfl ¼

0.025

0.025

0.02

0.15 0.015

0.01

0.005

0.1 0.1

0

Elevation difference (m)

0

0.005

0.01

0.05

0.015

0.02

0.025 7.5

00

0.025 7.5

0.05

0.1 -0.1 00

8

8.5

9

9.5

10

10

0.002 2

x 10

3

1.5

1

2

3

4

5

6

5 Time (days)

7

8

9

10

10

0.001 1

0.5

0 0

0.5

-0.001 1

Difference between Iterations 0 & 1 Difference between Iterations 1 & 2 Difference between Iterations 2 & 3

9.4

9.5

9.5

9.6

9.7

9.8

9.9

10

10

Difference between Iterations 3 & 4 Difference between Iterations 4 & 5 Difference between Iterations 5 & 6

Fig. 3. Convergence of iteration scheme at Node 2 using a linear background field. The differences between consecutive elevation solutions are plotted. Clockwise: (a) all six iterations over the entire time domain; (b) an enlarged of the region indicated in panel (a); (c) an enlarged view of the region indicated in panel (b). For clarity, the difference between iteration 0 and 1 is not shown in panels (b) and (c) and the difference between iterations 1 and 2 is not shown in panel (c). The difference between Iterations 4 and 5 and 5 and 6 are not distinguishable from zero at these scales.

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341 Absolute value of elevation difference (m)

326

10 10 120

120

10 10 100

100

10 8010

80

10 6010

60

10 4010

40

10 2010

20

10 010

0

00

1

2

3

4 6 55 Time (days)

Difference between Iterations 0 & 1

7

8

9

10 10

Difference between Iterations 1 & 2

Fig. 4. Lack of convergence of iteration scheme at Node 2 using a resting background field. The absolute values of the differences between consecutive elevation solutions are plotted using a log scale in the vertical. The difference between Iterations 1 and 2 increases with time after approximately 1.3 days, leading eventually to computational overflow errors which force the iterations to cease.

of a sufficiently poor background field for this problem is a ‘‘resting’’ solution (elevation and velocity equal to zero for all space and time). As shown in Fig. 4, the difference between consecutive iterations increases with time; the linear solutions do not converge to any solution. Typically, iteration schemes require more iterations for longer time domains. However, this particular scheme seems to be quite robust in terms of simulation length. In fact, for the Bight of Abaco, at most one additional iteration is required when the simulation length is increased from 10 days to 300 days, regardless of convergence criterion used (see Table 1). [Throughout, the same numerical value is used for each of eg,eu and ev in (17) and (18); e.g., e = 1 · 103 implies that eg = 1 · 103 m and eu = ev = 1 · 103m/s.] The results in Table 1 also indicate that the convergence rate increases as the number of iterations increases. For example, for a 10 day simulation, the e = 1 · 102 criterion is satisfied in five iterations. But, one additional iteration satisfies both the e = 1 · 103 and 1 · 104 criteria, and a final iteration satisfies the, e = 1 · 105, e = 1 · 106, Table 1 Bight of Abaco: number of iterations required for convergence for varying simulation lengths and a range of convergence criteria Convergence criterion 2

1 · 10 1 · 103 1 · 104 1 · 105 1 · 106 1 · 107

Simulation length (days) 10

50

100

150

200

250

300

5 6 6 7 7 7

5 6 6 7 7 7

5 6 6 7 7 7

6 6 7 7 7 8

6 6 7 7 7 8

6 6 7 7 7 8

6 6 7 7 7 8

Cell blocks are color coded; the warmer the color, the more iterations required. The same color code is used in reporting all results to ease comparison of results between experiments.

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

327

and e = 1 · 107 criteria. Throughout the experiments presented here, if convergence is reached, the rate of convergence typically increases with the number of iterations. Consider now the impact of the linear bottom friction parameter on the iteration scheme. Recall, the magnitude of the linear bottom friction term should be the same order of magnitude as the nonlinear bottom friction term; thus the linear bottom friction parameter is not necessarily equal to the nonlinear bottom friction parameter, but rather the two bottom friction parameters are related by a scale factor a: cfl ¼ acf

ð22Þ

n

where a = B /h. The number of iterations required for convergence using a range of values of the friction scale factor a and a range of convergence criteria is shown in Table 2. As mentioned earlier, a reasonable value for a in this case, based on depth and velocity scales of the Bight of Abaco, is about 0.05; Table 2 shows that this choice results in the minimum number of iterations regardless of the convergence criterion. The iteration scheme fails to converge for a 6 0.01 of a P 2.0. In order to gain some insight into this behavior, time series of the background elevation fields (at Node 2) computed with a range of a values is plotted in Fig. 5; the nonlinear field and the resting background field are also included for comparison purposes. As expected, the amplitude of the background amplitude decreases as a increases, and in fact, the background elevation is closest to the nonlinear solution when a = 0.05 and a = 0.1; these yield background elevation amplitudes that are, respectively, roughly increased by 30% and roughly decreased by 30% when compared to the nonlinear elevation amplitudes. The smallest value of a that leads to convergence is 0.02, which yields background elevation amplitudes that are increased by roughly 100% when compared to the nonlinear solution; when a is further reduced to 0.01, the elevation amplitudes increase by roughly 200% when compared to the nonlinear solution and the iteration algorithm fails to converge. Similarly, the largest value a of for which convergence is reached is 1.0; the background elevation amplitude is reduced by roughly 80%; when a is increased again to 2.0, the background elevation amplitude is reduced by 90% and the iteration scheme fails to converge. Thus, from these experiments on the Bight of Abaco, it seems that the convergence occurs when the amplitude of the background field is between roughly 80% and +100% of the amplitude of the nonlinear field. Table 2 Bight of Abaco: the number of iterations required for convergence using a range of values of the friction scale factor a and a range of convergence criteria Convergence criterion, e

a 0.01

0.02

0.05

0.1

0.2

0.5

1.0

2.0

1 · 101 1 · 102 1 · 103 1 · 104 1 · 105 1 · 106 1 · 107

NC NC NC NC NC NC NC

6 7 7 8 8 9 9

4 5 6 7 7 7 7

4 5 6 7 7 7 8

5 7 7 8 8 9 9

7 8 9 9 10 10 10

8 9 10 11 12 12 12

NC NC NC NC NC NC NC

‘‘NC’’ indicates that convergence is not reached (See Fig. 4 for an example of the behaviour of the sequence of solutions when convergence is not reached).

328

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341 1

0.8 b) 0.8

0.6

1 0.8 0.8

0.4

a)

0.4

0.2

0.0

Elevation amplitude (m)

0

0.6

0.2

0.4 0.4

-0.4

0.4

0.6

0.2

-0.8 1 0.8

0.00

-0.4 0.4

0.6

1

2

3

4 6 55 Time (days)

7

8

9

1.2

1.4

1.6

1.8

2

2.2

2.4

2.6

2.8

Solution to nonlinear eqns Solution to linear eqns, α = Solution to linear eqns, α = Solution to linear eqns, α = Solution to linear eqns, α = Solution to linear eqns, α = Solution to linear eqns, α = Resting solution

0.2

0.8 -0.8 00

1

10 10

3

3

0.01 0.02 0.05 0.1 1.0 2.0

Fig. 5. Time series at Node 2 of background field computed using linear equations (18) and (19) with a range of linear friction scaling coefficients: (a) the entire time series, while (b) an enlarged view of the time series from day one through day three. The nonlinear field and the resting background field are also included for comparison purposes. The amplitude of the background field decreases as a increases.

Consider now the case when the background field is determined by solving the linear equations (18) and (19) with the amplitude of the forcing differing from the forcing amplitudes used in the dataless iterations (but the phases remain the same as those in the dataless iterations and a = 0.05). The results, shown in Table 3, indicate that the scheme is quite tolerant of small amplitude differences. For example, at most one additional iteration is required for amplitudes increased by as much as 60% or reduced by as much as 60%, regardless of the convergence criteria. However, for amplitudes that are increased by 140% or more or reduced by 100% (i.e., resting), no converged solution is reached. These results are consistent with the findings associated with the bottom friction experiments just presented. Consider now the case when the background field is determined by solving the linear equations with phases that differ from the phases used in the dataless iterations (but the amplitudes remain

Table 3 Bight of Abaco: number of iterations required for convergence for background fields with amplitude difference from dataless iterations and a range of convergence criteria Convergence criterion, e 0

1 · 10 1 · 101 1 · 102 1 · 103 1 · 104 1 · 105 1 · 106 1 · 107

Amplitude difference 100

80

60

40

20

0

20

40

60

80

100

120

140

NC NC NC NC NC NC NC NC

2 6 7 8 8 9 9 9

2 4 6 7 7 7 8 8

2 4 5 6 6 7 7 7

2 4 5 6 7 7 7 7

2 4 5 6 7 7 7 7

2 4 5 6 7 7 7 8

2 5 6 6 7 7 8 8

2 5 6 7 7 8 8 8

3 5 6 7 7 8 8 8

3 5 6 7 8 8 8 8

4 6 7 7 8 8 8 9

NC NC NC NC NC NC NC NC

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

329

Table 4 Bight of Abaco: number of iterations required for convergence for background fields with phase difference from dataless iterations and a range of convergence criteria Convergence criterion, e 1

1 · 10 1 · 102 1 · 103 1 · 104 1 · 105 1 · 106 1 · 107

Phase difference (degrees) 90

40

20

0

20

40

90

180

7 8 9 9 10 10 10

5 7 8 8 9 9 9

5 6 7 8 8 9 9

4 5 6 7 7 7 7

4 6 6 7 7 8 8

5 7 8 8 9 9 10

7 8 9 10 10 11 11

8 9 10 11 11 11 11

the same as those in the dataless iterations). The number of iterations as a function of both phase difference and convergence criteria is shown in Table 4. Although a phase difference increases the number of iterations required for convergence for a given convergence criterion, convergence is eventually reached regardless of how large the phase difference is. Next, consider the case when the background field is computed by solving the linear equations with both phases and amplitudes that differ from the phases and amplitudes used in the dataless iterations. The results are shown in Table 5. As expected, the number of iterations required increases with differences in both amplitude and phase. When both amplitude and phase differences are both large, the iteration scheme does not converge. In fact when the background field is exactly out of phase (phase difference of 180), the iteration scheme, in most cases, cannot recover the nonlinear solution if the amplitude difference is large as well. The exception is the one case where the amplitude difference is 80%; the fact that this case converges (albeit in 11 iterations) demonstrates the somewhat unpredictable behavior of the iteration scheme. These tests were repeated with an additional 25 knot (12.8 m/s) wind forcing from the east. The same forcing (all tidal constituents and wind) is included both in the background field and the Picard iterations. The results (not shown) are qualitatively similar to those just described, although typically one additional iteration is required. These tests were repeated on the Bight of Abaco with a background field that is forced only by the M2 tide while the Picard iterations are forced by all tidal constituents. The results (not shown) are similar to those just presented for background field forced with all of the constituents used in Table 5 Bight of Abaco: number of iterations required for convergence for background fields with phase and amplitude difference from dataless iterations, with a convergence criterion of e = 1 · 103 Phase difference (degrees) 0 20 40 90 180

Amplitude difference (%) 0

20

40

60

80

100

120

6 6 8 9 10

6 7 8 9 NC

6 7 8 9 NC

7 7 8 10 NC

7 7 8 10 11

7 8 8 10 NC

7 8 9 NC NC

330

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

the dataless iterations, although one additional iteration is sometimes required to reach the same convergence criterion. However, when both phase and amplitude differences are large, the required iterations are somewhat unpredictable. For example, when the background field is forced only by the M2 tide, and the phase difference is 180, the iteration scheme converges (in 11 iterations) for 100% and 120% amplitude differences, whereas these same experiments failed to converge when the background field was forced by all constituents (Table 5).

6. Idealized inlet (advective terms dominate) When advective terms dominate the flow physics, ADCIRC typically requires a very small time step owing to the explicit treatment of that term and the resulting numerical stability limitations. Here, the robustness of the iteration scheme is examined using an Idealized Inlet domain with very large advective terms. The geometry of the Idealized Inlet is shown in Fig. 6a. All boundaries are land (no-flow) boundaries except for the southern boundary, which is an open ocean boundary forced by an M2 tide with an amplitude of 0.15 m. The bathymetry is relatively shallow: the depth within the

Fig. 6. The Idealized Inlet: (a) bathymetry in meters and (b) finite element grid with nodal spacing ranging from 0.1 km to 1 km. The inset shown in (a) will be used in Fig. 10. In (b), Node 1 is located in the middle of the constriction, and Node 2 is located just inside and to the right of the inlet opening. The depth at both nodes is 5.0 m. Term balances at these nodes are shown in Fig. 7.

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

331

inlet is 5 m, and the depth increases linearly from 5 m at the mouth of the inlet to 14 m at the open boundary. The finite element grid (Fig. 6b) has 11,768 nodes and 23,104 elements, with the distance between nodes ranging from 0.1 km to 1.0 km. The simulation is spun up from rest. For two simulation days, the simulation is ‘‘ramped’’; i.e., the forcing is multiplied by a hyperbolic tangent function. Then the simulation is run using the full forcing for another 1.6225 days, for a total simulation length of 3.6225 days. The time step is 10.916 s. There is no horizontal eddy

Scaled term magnitude for wave equation

Node 1 (depth = 5.0 m)

Node 2 (depth = 5.0 m)

1

1

1

0.8

0.8

0.6

0.6

0.4

0.4

1

0.2

0.2

0

0

0

0

0.2

0.2

0.4

0.4

0.6

0.6

0.8

0.8

-1 1 1

2.6

2.7

2.8

2.9

3

3.1

3.2

3.3

3.4

3.5

3.6

Scaled term magnitude for x-momentum equation

1

1

2.6

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

0

0.2

0.4

0.4

0.6

0.6

0.8

0.8

-1 1

-1 1

2.7

2.8

2.9

3

3.1

3.2

3.3

3.4

3.5

3.6

1

1

2.6

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

3.1

3.2

3.3

3.4

3.5

3.6

2.7

2.8

2.9

3

3.1

3.2

3.3

3.4

3.5

3.6

2.8

2.9

3

3.1

3.2

3.3

3.4

3.5

0

0.2

0.2

0.4

0.4

0.6

0.6

0.8

0.8

2.6

3

0

0

1

2.9

1

0.8

-1

2.8

0

0.2

2.6

2.7

1

0.8

1

Scaled term magnitude for y-momentum equation

-1 1

2.7

2.8

8

2.9

3

3.1

3.2

9 Time (days)

3.3

3.4

3.5

3.6

10

-1 1

2.6

2.7

8

9

3.6

10

Time (days)

Surface gradient term

Finite amplitude term

Bottom friction term

Advection term

Time derivative term

Coriolis term

Fig. 7. Scaled term balances for the Idealized Inlet at a shallow Nodes 1 and 2. See Fig. 6 for node locations. Only the relative magnitude of terms within a figure are significant. Note the time derivative in the wave equation is o2g/ot2 and the time derivative in the momentum equation is oV/ot. For clarity, within each figure, terms with magnitudes less than 10% of the largest term in that figure are not shown.

332

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

viscosity, the nonlinear bottom friction coefficient, cf, is 0.0025, and the Coriolis parameter, f, is ^ The background field is the solution to the linear governing equations (18) and 8.59 · 105 s1 k. (19), with cfl = 0.05cf. Time series of term balances in the nonlinear solution at Nodes 1 and 2, shown in Fig. 7, demonstrate that advection is the dominant term in the wave equation at both nodes, and is predominantly balanced by the surface gradient term at Node 1 (at the center of the inlet) and by the surface gradient and bottom friction at Node 2 (just inside and to the right of the inlet). Advection is also significant in the momentum equations at Node 2, with nonlinear bottom friction also making an important contribution to the term balance. The iteration scheme converges to the nonlinear solution with e = 1 · 104 in six iterations. The impact of the background field on convergence is studied further in Section 9.

7. Implications of ‘‘dataless iterations’’ These results are encouraging, particularly given that these are especially stringent tests. Experience shows that iteration schemes often converge more swiftly within the context of the data assimilation algorithm (that is, with data) than in dataless iterations because the iterations are weakly constrained to come close to the data and thus the data can guide the iterations (e.g., Ngodock et al., 2000). One might well expect the performance of a particular iteration scheme to improve once it is introduced into the assimilation algorithm with data. Given these encouraging results, good convergence behavior is anticipated for inversions with data. 8. Memory/disk requirements The iteration scheme described thus far requires that the entire solution field at every time step associated with the n  1 iteration be either saved in memory or written to disk so that it is available when computing the nth iteration. For large space–time domains, this requirement becomes expensive. In particular, for highly nonlinear flows, the time step constraint required for stability may result in a prohibitively large number of time steps. In the Idealized Inlet, for example, more than 4000 time steps are required per M2 period when the largest allowable time step is used. As an alternative to saving the entire solution associated with the n  1 iteration, one might choose to save only every Nth time step (where N is a positive integer). Then, the n  1 solution at intermediate time steps can be interpolated from the saved values when required by the Picard iterations. Two potential techniques for interpolating the intermediate values are illustrated in Fig. 8. The first is straightforward piecewise linear interpolation (LI, in blue) (e.g., Weaver et al., 2003). While piecewise linear interpolation leads to a continuous approximation to the n  1 solution, the derivatives are not continuous. A second approach involves Hermitian cubic interpolation (HCI, in red) and yields a continuous approximation with continuous derivatives f ðxÞ ¼

2 X j¼1

fj u0j ðxÞ þ

df j u ðxÞ dx j 1j

ð23Þ

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

333

f ( x)

1

0

-1 0 0

0.1

0.2

0.2

0.3

0.4

0.4 x

0.5

0.6

0.6

0.7

0.8

0.8

Fig. 8. Illustration of memory/space-saving algorithms. The black line is a sine wave, f = sin(2p(x + 0.15)), the blue line is the linear interpolation of the sine wave with interpolation points (blue circles) separated by Dx = 0.2, and the red line is the Hermitian cubic interpolation of the sine wave with interpolation points (red dots) separated by 2Dx = 0.4. Both interpolation schemes have the same memory/disk requirements (For interpretation of colour the reader is referred to the web version of this article).

where subscript 1 indicates the beginning of the interpolation interval, subscript 2 indicates the end of the interpolation interval and the Hermitian cubic bases are (e.g., Celia and Gray, 1992) u01 ðxÞ ¼ u02 ðxÞ ¼

ðx2  xÞ2 ðx2  x1 Þ3 ðx  x1 Þ2 ðx2  x1 Þ3

u11 ðxÞ ¼ 

½3ðx2  x1 Þ  2ðx2  xÞ ½2ðx2  xÞ þ ðx2  x1 Þ

ðx2  xÞ2 ðx  x1 Þ ðx2  x1 Þ2

;

u12 ðxÞ ¼

ð23aÞ ðx  x1 Þ2 ðx2  xÞ ðx2  x1 Þ2

ð23bÞ

The derivatives df/dx are approximated using a second ordered centered difference of function values at the n  1 iteration. Note that HCI requires that both the function and its derivative be saved periodically, while LI requires only the function values. Thus, for the same interval N, HCI requires twice the storage than does LI. In all comparisons made here, the HCI interpolating interval is double the LI interpolating interval so that the storage requirements are the same. The errors introduced by LI and HCI are examined here. [A third approach to reducing memory/disk requirements, called ‘‘checkpointing’’ has been widely used within the variational assimilation framework. In general terms, the time domain is divided into a number of windows, such that the solution within one window fits in memory. First, the entire solution is determined by time stepping, and the model state is saved to disk at time steps coinciding with the beginning of each window, but not at other time steps. Then, the solution can be recomputed within each window, as needed. Within a tangent linearization framework, checkpointing obviates the need to save the entire n  1 solution but requires that it be computed twice: first to periodically save the model state, and then again to recompute intermediate values as required to compute the nth iteration, thus doubling the CPU cost. The benefit is that no error is introduced. See Marotzke

334

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

et al. (1999); Griewank and Walther (2000); Ngodock (2004) and Charpentier (2001) for examples of checkpointing within the variational assimilation framework. The Idealized Inlet is used to illustrate the errors associated with the LI and HCI methods. Maximum elevation error is defined as gki  gki j; Eg;max  max j ~ ~ gki

i ¼ 1; NN ;

k ¼ 1; NK

ð24Þ gki

where is the elevation field computed with LI or HCI, is the elevation field computed by saving every time step (that is, N = 1), NN is the number of nodes in the finite element grid and Nk is the number of time steps. Maximum errors for u and v are defined analogously. The maximum error for each of g,u and v for both LI and HCI is plotted vs. N in Fig. 9. The HCI method has significantly smaller maximum errors, particularly for large N. It is also interesting to point out that, for all values of N, LI requires seven iterations, but HCI requires only six. Thus, even though the Hermitian cubic method requires nominally more computation than the linear method, the additional expense is easily offset because it requires fewer iterations. Furthermore, within each iteration, there is less reading and writing to disk, and the testing of the convergence criteria is twice as fast because there are half as many saved function values to examine. The result is that, in our tests, HCI actually required less overall computation time than did LI.

Fig. 9. Maximum error for g, u and m vs. N for LI (blue, using bottom axis) and HCI (red, using top axis) interpolation methods. Eg,max in solid lines (m), Eu,max in dashdot lines (m/s) and Em,max in dotted lines (m/s). Memory/CPU requirements for LI and HCI are the same along vertical lines. The maximun errors are relatively close for small N, but for large N, the HCI method has significantly smaller maximum errors (For interpretation of colour the reader is referred to the web version of this article).

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

335

Table 6 Tidal constituents included in harmonic analysis, and number of points per period for LI and HCI methods Constituent

Steady

M2

M4

M6

M8

M10

Period (h) Interpolation intervals (LI, N = 100) Interpolation intervals (HCI, N = 200)

–

12.42 41.0 20.5

6.21 20.5 10.3

4.14 13.7 6.8

3.11 10.3 5.1

2.48 8.2 4.1

Next, the last day of the Idealized Inlet solution was harmonically analyzed (e.g., Emory and Thomson, 1998; Godin, 1972) for the M2 tide and its overtides. The constituents included in the harmonic analysis are listed in Table 6, along with the periods and the number of interpolation intervals per period for both LI and HCI. Results (not shown) indicate that both LI and HCI are quite accurate for the constituents with longer periods, although HCI always outperforms LI by about an order of magnitude for both elevation (amplitude and phase) and velocity (major-semi axis, major-semi axis phase, eccentricity and major-semi axis direction). However, for overtides with relatively short periods (e.g., M10), and therefore a relative small number of interpolation intervals per period (roughly eight for LI and 4 for HCI), the errors in LI become

Fig. 10. Errors resulting for M10 overtide using LI and HCI schemes for the region around the opening to the inlet (inset region indicated in Fig. 6). Top two panels are relative elevation amplitude error (%) and the bottom two panels are elevation phase error. The left panes errors resulting from the LI method (N = 100) and the right two panels are errors resulting from the HCI method (N = 200). The LI method introduces errors in relative elevation amplitude ranging from 2.1% to 2.6% (a), while the HCI method introduces errors ranging from 0.020% to 0.051% (b). The LI method introduces errors in elevation phase ranging from 2 to 0.42 (c), while the HCI method introduces errors ranging from 0.024 to .014 (d).

336

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

large enough that they could be significant, especially in the region near the opening to the inlet, as shown in Fig. 10. The errors in the HCI solution in this region are about two orders of magnitude less than those in the LI solution. Errors in the velocity components (not shown), even for overtides with short periods, are small enough to be insignificant even for LI (e.g. O(101) percent for relative major-semi axis, O(101) degrees for major-semi axis phase, O(102) for eccentricity and O(101) degrees for major-semi axis direction). Although HCI requires nominally more computation for the interpolation LI, our experience shows that this additional work is typically offset by fewer iterations required to reach a particular convergence criterion. Given the improvement in elevation amplitude and phase results for short period constituents, we recommend the HCI method when short period constituents are important components of the solution. We expect that these results should hold for a wide variety of ocean models other than SWE models.

9. Tests of background field for idealized inlet using HCI method In Section 5, the ability of the iteration scheme to recover the nonlinear solution from various background states was examined within the context of the Bight of Abaco. Owing to memory/ CPU limitations, these tests could not be repeated for the Idealized Inlet without some kind of interpolation method, and based upon the results just presented, we use HCI with N = 200. Given that the maximum elevation and velocity errors increase with increasing N (see Fig. 9), the minimum value of convergence criterion that we implement is 1 · 105, as compared to the Bight of Abaco where the minimum value of e was 1 · 107. First, reconsider the impact of the linear bottom friction parameter on the iteration scheme. Recall that cfl = acf, where a = Bn/h. The number of iterations required for various values of a and a range of convergence criteria is presented in Table 7. Compared to the analogous results for the Bight of Abaco (Table 4), the number of iterations required is roughly the same when the estimate of a is relatively good, but when the estimate of a is poor, the Idealized Inlet simulations require substantially fewer iterations. More striking is the fact that convergence is reached for the Idealized Inlet simulation regardless of the values of a. The more robust behavior of the Table 7 Idealized Inlet: the number of iterations required for convergence using a range of values of the friction scale factor a and a range of convergence criteria Convergence criterion, e 1

1 · 10 1 · 102 1 · 103 1 · 104 1 · 105

a 0.01

0.02

0.05

0.1

0.2

0.5

1.0

2.0

4(NC) 5(NC) 6(NC) 6(NC) 6(NC)

4(2) 5(2) 5(2) 6(2) 6(2)

4(0) 5(0) 5(1) 6(1) 6(1)

4(0) 5(0) 5(1) 6(1) 6(1)

4(1) 5(2) 6(1) 7(1) 7(1)

5(2) 6(2) 6(3) 7(2) 7(3)

5(3) 6(3) 6(4) 7(4) 7(5)

5(NC) 6(NC) 6(NC) 7(NC) 7(NC)

Number in parenthesis is the difference between the number of iterations required for the Bight of Abaco (see Table 4) and the number of iterations required for the Idealized Inlet. ‘‘(NC)’’ indicates the analogous Bight of Abaco simulation did not converge.

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

337

iteration scheme for the Idealized Inlet in the face of a poor estimate of the bottom friction most likely owes to the fact that the magnitude of the bottom friction term is not as large in the Idealized Inlet as in the Bight of Abaco, particularly in the wave equation (compare the top two panels in Figs. 2 and 7). Consider now the case when the background field is computed using the linear governing equations (18) and (19), but with tidal forcing that differs from that used in the dataless iterations. While the Bight of Abaco was forced with five tidal constituents, the Idealized Inlet is forced only with the M2 tide, with overtides generated by the nonlinearities. When there is amplitude difference in the background field (as compared to the amplitude used in the dataless iterations), the number of iterations required for the Idealized Inlet simulations are typically the same or one fewer than the number required for the Bight of Abaco, as shown in Table 8. However, the Bight of Abaco simulations converged with amplitude differences as large as 120%, whereas the Idealized Inlet simulations failed to converge for amplitude differences of 80% and larger. On the other hand, the Idealized Inlet iterations converge when the amplitude is decreased by even 100%; the solution in this case is the resting solution. When there is a phase lag or lead in the background field, the iteration scheme again converges, regardless of the size of the phase difference, and convergence is reached with one or two less iterations than required for the same phase difference and convergence criterion (results not shown). Finally, the background field is computed by solving the linear equations with both phases and amplitudes that differ from the phases and amplitudes used in the dataless iterations, and the results are presented in Table 9. These simulations also fail to converge for amplitude differences of Table 8 Idealized Inlet: number of iterations required for convergence for background fields with amplitude difference from dataless iterations Convergence criterion 1 · 101 1 · 102 1 · 103 1 · 104 1 · 105

Amplitude difference (%) 100

80

60

40

20

0

20

40

60

80

5(NC) 6(NC) 7(NC) 7(NC) 7(NC)

4(2) 6(1) 6(2) 7(1) 7(2)

4(0) 5(1) 6(1) 6(1) 6(1)

4(0) 5(0) 5(1) 6(0) 6(1)

3(1) 4(1) 5(1) 6(1) 6(1)

4(0) 5(0) 5(1) 6(1) 6(1)

4(0) 5(0) 6(0) 6(1) 6(1)

4(1) 5(1) 6(0) 6(1) 7(0)

4(1) 5(1) 6(1) 6(1) 7(1)

NC NC NC NC NC

Table 9 Idealized Inlet: number of iterations required for convergence for background fields with phase and amplitude difference from dataless iterations, with a convergence criterion of e = 1 · 103 Phase difference (degrees) 0 20 40 90 180

Amplitude difference (%) 0

20

40

60

80

5(1) 6(0) 7(1) 7(2) 9(1)

6(0) 6(1) 7(1) 8(1) 9(NC)

6(0) 6(1) 7(1) NC 10(NC)

6(1) 6(1) 9(+1) NC NC

NC NC NC NC NC

338

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

80% or more. Experiments (not shown) indicate that the lack of convergence here does not owe to the large interpolation interval. Otherwise, these results are quite similar to those obtained for the Bight of Abaco, with more iterations required as both the phase and amplitude difference increase. 10. Conclusions We have found that, in a variety of compelling situations, the linear Picard iteration scheme based upon the tangent linear operator converges to the solution of nonlinear governing equations. We are therefore encouraged to pursue 4D variational assimilation using this iteration scheme. The tangent linearization fails for the bottom friction term owing to the magnitude of the velocity appearing in the denominator of one of the linearized bottom friction terms. At a stagnation point, this term becomes undefined. However, if the denominator of that problematic term is given a small, positive minimum value, the term remains finite, even at stagnation points. This small modification causes no distinguishable change in the nonlinear solution. When nonlinear terms dominate the dynamics, experimental results show that a good background field in the first iteration is prerequisite to convergence. The solution to the linear governing equations is successfully used here as a background field in the first iteration, as well as several variations of that field. In the Bight of Abaco, where bottom friction is one of the leading terms in the governing equations, it is found that a reasonable value must be used for the linear bottom friction parameter in order for the sequence of linear equations to converge to the nonlinear solution. Otherwise, the sequence fails to converge in the sense that the difference between successive iterations increases dramatically with time, eventually leading to computational overflow errors. Experiments also demonstrated that the iteration scheme is typically able to recover from a background field that is forced differently than the nonlinear equations, and therefore the dataless iterations, and is particularly adept when there is a phase lag or lead between the background field and the nonlinear equations, although the bigger the phase lag or lead, the more iterations are required. Differences in the amplitude of the forcing were somewhat more problematic, as the iteration failed to converge when the forcing used to generate the background field was either too large, or too small, with the exact limits being somewhat domain dependent. It is recognized that the memory/disk requirements for this iteration scheme (indeed, any iteration scheme) may be prohibitively large, and two modifications which save the solution periodically are investigated. When the intermediate solution values are interpolated from the saved values as they are needed, the errors introduced by the approximation are relatively small. Two interpolation methods are examined here: the well-known linear interpolation method and a new Hermitian cubic method. Overall, the latter had smaller errors, typically by one or two orders of magnitude. These errors were particularly significant for short period overtides, such as M10. Although in these experiments the errors were not particularly large, they could be significant when highly accurate solutions are desired. The nominal additional computation required for the higher degree interpolation is compensated in these experiments by faster convergence and therefore fewer iterations. Regardless of the interpolation method used, one must be cautious to ensure that an appropriate interval between saved model states is used so that higher frequency dynamics are resolved adequately. The appropriate interval will ultimately be a compromise

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

339

between the importance of higher frequency dynamics and the availability of computational resources.

Acknowledgments The authors thank Richard A. Luettich, Jr. and Andrew F. Bennett for their valuable contributions to this paper. This material is based upon work supported by the National Science Foundation under Grant no. OCE-0121315 and by the Office of Naval Research under Grant no. N000140410439.

Appendix A. Tangent linear equations The linearized wave equation is o2 gn ogn þ s  ½ðVn  Vn1 ÞH n1 þ Vn1 H n 
rs0 0 ot2 ot 
 n1 n n n1 og n1 og þV þ r
ðV  V Þ ot ot  n1 n  H ðV  Vn1 Þ
rVn1 þ ðH n  H n1 ÞVn1
rVn1 þ H n1 Vn1
rVn

  n1 ps n n1 n1 n n n1  H f  ðV  V Þ þ H f  V þ gðg  anÞ  H r q0  og ss n1 n1 n1 n H rðgg Þ þ H rðgg Þ  Er þ  cf Bn1 Vn ot q0 
  n cf n1 n1 n1 n n1 n n1 n1 n n1 n1  n1 u ðu  u Þ þ m ðm  m Þ V þ s0 ðV  V ÞH þ V H ¼0 B

vW L ðg;VÞ ¼

and the linearized nonconservative momentum equation is


 oVn  n ps n1 n1 n1 n n n Mðg; VÞ ¼ þ ðV  V Þ
rV þ V
V þ f  V þ r þ gðg  anÞ ot q0 " # E E n1 n1 n n1 n1 n n n1 n1 ðH  H ÞD½V H   n1 D½ðV  V ÞH þV H  þ H ðH n1 Þ2 " #
ss 2 Hn cf Bn1 ðH n  H n1 Þ n1  þ  V þ q0 H n1 ðH n1 Þ2 H n1 H n1  1 n1 n1 n n n1 n1 n n1 ðA:1Þ þ n1 2 ½ðu  u Þu þ ðm  m ÞV V þ V ðB Þ

340

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

References Bennett, A.F., Thorburn, M.A., 1992. The generalized inverse of a nonlinear quasi-geostrophic ocean circulation model. Journal of Physical Oceanography 22, 213–230. Bennett, A.F., Leslie, L.M., Hagelberg, C.R., Powers, P.E., 1993. Tropical cyclone prediction using a barotropic model initialized by a generalized inverse method. Monthly Weather Review 121, 1714–1729. Bennett, A.F., Chua, B.S., Leslie, L.M., 1996. Generalized inversion of a global numerical weather prediction model. Meteorology and Atmospheric Physics 60, 165–178. Bennett, A.F., Chua, B.S., Leslie, L.M., 1997. Generalized inversion of a global numerical weather prediction model, II: analysis and implementation. Meteorology and Atmospheric Physics 62, 129–140. Bennett, A.F., Chua, B.S., Harrison, D.E., McPhaden, M.J., 1998. Generalized inversion of Tropical AtmosphereOcean data and a coupled model of the Tropical Pacific. Journal of Climate 11, 1768–1792. Bennett, A.F., 2002. Inverse Modeling of the Ocean and Atmosphere. Cambridge University Press, United Kingdom, 234 pp. Celia, M.A., Gray, W.G., 1992. Numerical Methods for Partial Differential Equations. Prentice Hall, NJ. Charpentier, 2001. Checkpointing schemes for adjoint codes: application to the meteorological model MESO-NH. SIAM Journal on Scientific Computing 22, 2135–2151. Chua, B.S., Bennett, A.F., 2001. An inverse ocean modeling system. Ocean Modelling 3, 137–165. Egbert, G.D., Bennett, A.F., Foreman, M.G.G., 1994. TOPEX/POSEIDON tides estimated using a global inverse method. Journal of Geophysical Research 106, 22,475–22,502. Emory, W.J., Thomson, R.E., 1998. Data Analysis Methods in Physical Oceanography. Pergamon, Oxford. Godin, G., 1972. The Analysis of Tides. University of Toronto Press, Toronto. Griewank, A., Walther, A., 2000. Algorithm 799: an implementation of checkpointing for the reverse or adjoint mode of computational differentiation. ACM Transactions on Mathematical Software 26, 19–45. Hench, J.L., Luettich, R.A., 2003. Transient tidal circulation and momentum balances at a shallow inlet. Journal of Physical Oceanography 33, 913–932. Kinnmark, I., 1986. The Shallow Water Wave Equations: Formulation, Analysis and ApplicationLecture Notes in Engineering, vol. 15. Springer-Verlag, Berlin, 187 p. Kolar, R.L., Gray, W.G., Westerink, J.J., Luettich Jr., R.A., 1994. Shallow water modeling in spherical coordinates: equation formulation, numerical implementation and application. Journal of Hydraulic Research 32, 3–24. Le Dimet, F., Talagrand, O., 1986. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects. Tellus Series A 38, 97–110. Le Dimet, F.X., Ngodock, H.E., Luong, B., Vernon, J., 1997. Sensitivity analysis in variational data assimilation. Journal of the Meteorological Society of Japan 75, 245–255. Luettich, R.A., Westerink, J.J., Scheffner, N.W., 1991. ADCIRC: an advanced three-dimensional circulation model for shelves, coasts and estuaries, Report 1: Theory and Methodology of ADIRC-2DDI and ADCIRC-3DL. Dredging Research Program, Technical Report DRP-92-6, Department of the Army, Washington, DC, 150p. Luettich, R.A., Jr., Westerink, J.J., 2000. ADCIRC: A (parallel) advanced circulation model for oceanic, coastal and estuarine waters, Consolidated and updated theory and users manual. Available from: . Luettich Jr., R.A., Carr, S.D., Reynolds-Fleming, J.V., Fulcher, C.W., McNinch, J.E., 2002. Semi-diurnal seiching in a shallow, micro-tidal lagoonal estuary. Continental Shelf Research 22, 1669–1681. Lynch, D.R., Gray, W.G., 1979. A wave equation model for finite element tidal computations. Computers and Fluids 7, 207–228. Marotzke, J., Giering, R., Zhang, K.Q., Stammer, D., Hill, C., Lee, T., 1999. Construction of the adjoint MIT ocean general circulation model and application to Atlantic heat transport sensitivity. Journal of Geophysical Research 104 (C12), 29,529–29,547. McIntosh, P.C., 1987. Systematic design of observational arrays. Journal of Physical Oceanography 17, 885– 902. McNinch, J.E., Luettich, R.A., 2000. Physical processes around a capsulate foreland: implications to the evolution and long-term maintenance of cape-associated shoal. Continental Shelf Research 20, 2367–2389.

J.C. Muccino, H. Luo / Ocean Modelling 10 (2005) 316–341

341

Muccino, J.C., Hubele, N.F., Bennett, A.F., 2004. Significance testing for variational assimilation. Quarterly Journal of the Royal Meteorological Society 130, 1815–1838. Muccino, J.C., Bennett, A.F., 2002. Generalized inversion of the Korteweg-de Vries equation. Dynamics of Atmospheres and Oceans 35, 227–263. Ngodock, H.E., Chua, B.S., Bennett, A.F., 2000. Generalized inversion of a reduced-gravity primitive equations ocean model and tropical atmosphere ocean data. Monthly Weather Review 128, 1757–1777. Ngodock, H.E., 2004. Efficient implementation of covariance multiplication for data assimilation with the representer method. Ocean Modelling, in press. Oliger, J., Sundstrom, A., 1978. Theoretical and practical aspects of some initial boundary value problems in fluid dynamics. SIAM Journal on Applied Mathematics 35, 419–446. Uboldi, F., Kamachi, M., 2000. Time–space weak-constraint data assimilation for nonlinear models. Tellus 52A, 412– 421. Weaver, A.T., Vialard, J., Anderson, D.L.T., 2003. Three- and four-dimensional variational assimilation with a general circulation model of the tropical Pacific Ocean. Part1: Formulation, internal diagnostics and consistency checks. Monthly Weather Review 131, 1360–1378. Westerink, J.J., Luettich Jr., R.A., Scheffner, N.W., 1993. ADCIRC: an advanced three-dimensional circulation model for shelves, coasts and estuaries. Report 3: Development of a Tidal Constituent Data Base for the Western North Atlantic and Gulf of Mexico, DRP Technical Report DRP-92-6, Department of the Army, US Army Corps of Engineers, Waterways Experiment Station, Vicksburg, MS, 154p. Xu, L., Rosmond, T., Daley, R., 2004. Development of NAVDAS-AR: formulation and initial tests of the linear problem. Tellus.