Benchmarking a new open-source 3D circulation model (ELCIRC)

1791

Benchmarking a new open-source 3D circulation model (ELCIRC) Yinglong Zhang and Antdnio M. Baptista a* aCenter for Coastal and Land-Margin Research, Oregon Health & Science University, 20000 NW Walker Road, Beaverton, OR 97006 Released recently as an open source code, ELCIRC (Eulerian-Lagrangian Circulation) solves the primitive shallow-water Navier-Stokes equations with turbulence closure submodels. Numerically it uses a semi-implicit finite-difference/volume method on unstructured horizontal grids and structured grids in the unstretched vertical direction. An Eulerian-Lagrangian method (ELM) is used to treat the advection, and wetting and drying is a natural part of the algorithm. The model has low-order accuracy, but is very flexible, computationally efficient, and robust. Overall, it has shown excellent ability to address complex river-to-ocean systems, and is currently used as the computational engine for our observation and forecasting system for the Columbia River estuary and plume. As a part of the development of ELCIRC, we carefully assessed its performance against a wide set of controlled benchmark problems: wave propagation on a slope, geostrophic flow in a straight channel, and adjustment under gravity, etc. In this paper we report on these benchmark studies, which provide very useful insights on the capabilities and limitations of the model. Using carefully defined error metrics, convergence studies are carried out. Compared against well-established higher-order models (in particular, ADCIRC and 1ROMS), ELCIRC has the capability of compensating for its low-order accuracy through inexpensive high resolution. Benchmarks will be made available electronically before the publication of the conference proceedings. 1. I N T R O D U C T I O N

Recently the authors proposed ELCIRC (Eulerian-Lagrangian Circulation) [13] as a lower-order alternative to other benchmark ocean circulation models like POM [1], QUODDY [7], ROMS [4], etc. The major differences include" the use of hybrid triangular/quadrangular elements in the horizontal grid, unstretched z-coordinates in the vertical, Eulerian-Lagrangian treatment of the nonlinear advection term, and no mode splitting between external surface gravity and internal modes with the semi-implicit method. The result is an efficient and accurate, robust, and mass conservative model that is well suited for cross-shelf scale applications where processes of contrasting temporal and spatial scales are important. It has been successfully applied to the complex *The research was funded by NOAA (NA17FE1486, NA17FE1026, NA87FE0405), the U.S. Fish and Wildlife Service (133101J104), the Office of Naval Research (N00014-00-1-0301, N00014-99-1-0051), and the National Science Foundation (ACI-0121475).

1792 baroclinic circulation in Columbia River estuary and plume, and, via web dissemination (http'//www.ccalmr.ogi.edu/CORIE/modeling/elcirc/) and beta users, to a number of other field scale problems. While we have carried out extensive set of benchmark tests with ELCIRC before [13], the complexity of the model warrants more thorough exams, especially with regard to the following issues: the choice of triangular versus quadrangular elements, orthogonal versus non-orthogonal grids, optimal time steps in baroclinic applications, etc. The current paper is aimed at addressing some of these important issues. We briefly review essential aspects of the ELCIRC formulation in Section 2, and proceed in Section 3 to reporting the results of 3 targeted benchmark tests. The advantages/disadvantages of the formulation are summarized in the last section. 2. E L C I R C F O R M U L A T I O N 2.1. P h y s i c a l formulation The model solves the conventional set of Reynolds averaged N avier-Stokes equations with hydrostatic and Boussinesq approximations, together with the transport and conservation of salt and heat. In a Cartesian coordinate system (x, y, z) with the x-axis pointing to the east, y-axis to the north, and z-axis upward, the system of equations can be written as: (;Qw

V.u+

Oz = 0 ,

Orlo__{+ v . / HR+~ JHR-h

D UtD -- - f k

(1)

- o,

(2)

x u - V g ( rl - a ~ ) + -~o - -p-o

DS O(KhOS ) Dt = Oz ~ ,

DT O(OT) Dt=O

V p d z + -~z K , ~ -~z

,

(3)

(4)

Q (5)

where h(x, #) is the bathymetry, HR is the z-coordinate of the Mean Sea level (MSL), r/ is the free surface elevation, u is the horizontal velocity, w is the vertical velocity, f is the Coriolis constant, ~ is the earth tidal potential as defined in Reid [8], c~ is the effective earth elasticity, p is the water density (with a reference value of p0) which is related to the salinity S and temperature T via the equation of state of ISE80 standard, Pa is the atmospheric pressure, K,~ and Kh~ are the vertical eddy viscosity and diffusivity respectively, (~ is the rate of absorption of solar radiation, and Cp is the specific heat of water. We neglect the horizontal diffusion in the momentum and transport equations (3)-(5), primarily because the numerical diffusion may exceed its physical counterpart The equation system is complete with appropriate initial and boundary conditions, as well as turbulence closure schemes that determine the eddy viscosity and diffusivity. The vertical boundary conditions for the momentum and transport equations require the balance between internal and external turbulent fluxes, and horizontal boundary conditions

1793 can be chosen from an assortment of essential, natural and radiation conditions [13]. A variety of turbulence closure schemes are available in ELCIRC, ranging from simplest zero-equation to two-and-a-half-equation GLS model recently proposed by Umlauf and Burchard [10]. 2.2. N u m e r i c a l m e t h o d

Only an overview of the method is given here. A combination of horizontal triangular/quadr~ngular elements and vertical layers renders the domain into a series of prisms, the basic discretization unit used in ELCIRC. Note that the prisms may be cut off near the bottom and surface to account for unevenness there, resulting in a staircase representation. A semi-implicit finite volume method is applied to the integrated continuity equation (2) (and as a result, volume is conserved within each element and globally as well), and a semi-implicit finite difference is applied to the momentum equation (3) along each side of an element, with the elevation gradient being treated semi-implicitly and the total derivative approximated by an Eulerian-Lagrangian method (ELM) via backtracking along particle paths (or characteristic lines). Due to the semi-implicit ELM, the only remaining CFL restriction on the time step comes from the baroclinic term, and is very mild, resulting in greater efficiency. Substituting the inverted momentum equation into the integrated continuity equation leads to a single equation for the elevations, whose matrix is positive definite, sparse, and symmetric, and thus amenable to efficient solvers like Jacobian Conjugate Gradient method. After the elevations are found, the horizontal velocity is then solved from the momentum equation, and the vertical velocity from the continuity equation (1) using a finite volume method to conserve volume locally. A fully implicit finite difference method is used to solve the transport equations (4,5) with the total derivative again approximated by the ELM. A split of elements at the foot of a characteristic line before interpolation helps reduce numerical diffusion [13]. A similar technique is applied to the GLS closure equations. The advantage of the methodology adopted in ELCIRC includes the volume conservation, automatic treatment of wetting and drying, and ability to circumvent the CFL restriction commonly found in other benchmark models. The disadvantage lies in its lower-order nature and numerical diffusion associated with the ELM. The latter can be partially mitigated with grid refinement; the use of unstructured grids and large time steps bestows upon it a competitive edge over other models even under grid refinement. We shall demonstrate this in the next section with numerical examples. Due to the finite-difference attribute of the method, strictly speaking orthogonal grids are required, and circumcenters should be used inside each element for evaluating the elevation gradient. Second-order accuracy is achieved when uniform (triangular or quadrangular) grids are used, while the method is only first-order accurate for non-uniform orthogonal grids [2]. However, generating orthogonal grids is by no means a trivial endeavor, and so far there is no robust technique for accomplishing this, even for simple domains; see the discussions in Frey and George [3] for example. Furthermore, a more serious problem associated with non-uniform orthogonal grids is that the distance between two circumcenters of neighboring elements may be small, as in the case of two nearly right triangles. This leads to nearly singular matrix, and unpredictable results. In such a case, the two triangles should be combined into one quadrangle. We shall demonstrate in the

1794 following section that non-uniform orthogonal triangular grids in generally do not lead to better accuracy as compared to non-orthogonal grids using centroids as element centers, but rectangular elements (being orthogonal) are highly desirable in any application. Our final recommendation in terms of choosing a grid in ELCIRC is to use a combination of rectangular and general triangular elements, and to use centroids as anchor points for evaluating the elevation gradient. The ELM nature of the model bestows upon it a peculiar quality: the accuracy does not necessarily improve with the use of smaller time steps. It has been shown that large time steps are optimal for ELM alone, as this reduces numerical diffusion [13]. Therefore, for advection-dominated phenomena (such as the transport of salt and heat), the best numerical accuracy is achieved at a moderate time step due to the opposing effects of the time step on the advection and semi-implicit schemes. In other words, there exists an optimal time step for baroclinic problems. This will be shown in Section 3.3. 3. T E S T

RESULTS

3.1. W a v e r u n u p o n a l i n e a r s l o p e

While this problem has been well studied [13,11], here we are mainly concerned with the choice of orthogonal versus non-orthogonal triangular grids. The domain is a quarter annulus with inner and outer radii of 60.96 km and 152.4km, respectively. An M2 tide of amplitude 0.3048m is imposed at the outer boundary, while all other boundaries are perfectly reflective. The bottom is a frictionless linear slope that varies from 10.02m to 25.05m (from inner to outer boundary), and the problem is radially symmetric. The non-orthogonal grid we used divides the radial direction into 6 equal segments and the azimuthal direction into 8 segments (Fig. la). With the aid from a grid optimizer AGrid [12], an orthogonal grid is created based on the 6 x 8 grid (Fig. lb), with maximum internal angle being about 85 ~. The two grids have roughly same resolution. It is noteworthy that generating an orthogonal grid for this simple domain is by no means a trivial task; manual intervention was needed after AGrid was applied, which was often laborious. The 5-day runs were allowed to spin up for 2 days before a quasi-steady state is established, which is then compared to the linearized solution of Lynch and Gray [6]. For illustration purpose, we examine the errors in the predicted elevations at a point on the inner boundary (Fig. 2). Centroids are used as element centers for non-orthogonal grid, while circumcenters are used as element centers for the orthogonal grid. It is clear that the solution from the orthogonal grid is worse than that of the non-orthogonal grid. Refining the orthogonal grid by splitting each element into 4 generates another orthogonal grid of similar quality (Fig. lc), with slightly improved accuracy as compared to the non-orthogonal grid solution (Fig. 2). Therefore non-uniform orthogonal grids, with circumcenters being used as element centers, do not offer any advantage over general non-orthogonal grids. 3.2.

Geostrophic

flow in a c h a n n e l

In the absence of friction, a rapidly rotating fluid is geostrophic, with the Coriolis acceleration balanced by the pressure gradient, and isobars coinciding with streamlines. In the special case of a uniform geostrophic flow in a straight channel of constant depth, a lateral surface slope is developed perpendicular to the mean flow direction. The numerical

1795

,\\x \\\\

'\

(~)

(b)

,

\///

(c) Figure 1. Grids for a quarter annulus domain. (a) 6 x 8 non-orthogonal grid; (b) orthogonal grid; (c) orthogonal grid generated by splitting each element of the grid in (b) into 4.

1796

0,15 .....................

, ............

1

..........

'

....................

0.10

v

0,05

/~,~ ',

t m ',

"'"~

'/.'
'

j

.

2

9~

0.00

;:'

,

~",,/~~<~:"

:k

',

:/m ',

" 'G ',

:i.

',

/"2

',., ./

....:::

', /..?.y"

,

',

1)

m -0,05

10 L

|

-0.15

m.

Non-orthogonal grid (96 elements)

.... O~hogoooJr ...........

5,0

(90 ~i~.t.~) i 4,0

.j

I ...................

9 5,0

Time, (days)

Figure 2. Elevation errors at a point on the inner boundary of the quarter annulus domain.

model is set up with a channel 5kin long, 1kin wide and 50m deep. Uniform inflow is imposed at the left end with a total discharge of 49000m a/s, and elevation there is allowed to adjust by itself with the aid of a radiation boundary condition. The flow starts from rest and ramps up during the first day, shortly after which a steady state is established. The flow should be uniform everywhere with u = 0.98m/s; any departure from this value is taken as a numerical error. We consider two grids here. The first grid consists of uniform quads with Am = Ay = 100rn, and the second is an orthogonal triangular grid adapted from the uniform grid (Fig. 3), with the circumcenters being used as element centers. The two grids have roughly the same resolution. The best effort was made to make the orthogonal grid as uniform as possible, but generating a uniform orthogonal triangular grid even for this simple domain is extremely difficult, if at all possible. Table 1 shows the mean and standard deviation of the computed velocity for the two grids, and it is obvious that the uniform (orthogonal) quads did a better job. This is hardly surprising as second-order accuracy is achieved with a uniform orthogonal grid. Uniform quads are indeed highly desirable in general and provide an efficient way of covering a large open area with minimum number of elements. They are also more natural candidates in representing channels with steep side banks, as long as their sides are aligned along the channel boundary, because ELCIRC assumes constant depth for each side.

3.3. Baroclinic adjustment under gravity Ceostrophic adjustment under gravity of two fluids of different densities has been studied by a number of models [5]. As an effort to foster inter-model comparison, here we took as an example the problem from the ROMS web site [9]. The domain is a closed box of

1797

Figure 3. Orthogonal grid for a straight channel.

Table 1 Comparison of velocity errors for quadrangular and triangular grids for a geostrophic channel flow exact soln. Quads Orthogonal triangles Average velocity (m/s) 0.981 0.9803 0.9796 Standard deviation (m/s) 0 0.0016 0.0022

1798 Table 2 Specification of numerical parameters for adjustment under gravity Duration A x - Ay Az Tt=o St=o K,~ 12 hours 0 5km; 1kin 0 5m 4~ 3 12511 - tanh( x-x~ 0 9

"

"

1 0 0 0

Khv

f

Cd

0

0

0

Table 3 Errors in calculated internal wave speed with At - 5rnin Relative error (%) A z - A y - 0 . 5 k m Az-Ay-lkm Quad 15.10 20.00 Triangles 15.22 21.69

64krn x 20krn x 20rn, with two fluids separated by a vertical wall and each occupying half of the domain initially. After the wall is removed, two density fronts travel in opposite directions at the internal wave speed. The numerical parameters used in this problem are summarized in Table 2. Different combinations of grid sizes and time steps were attempted; the only results shown here relate to differences in the calculated internal wave speed. The internal wave speed is taken as the mean velocity at which the two fronts (or more precisely, the intersections between the two fronts and the surface/bottom respectively) travel during the 12 hour run. Both quads and triangles (generated from splitting each quad) are used, but the former yields better results (Table 3), again confirming the superior performance of the quads. We also tried an orthogonal grid, but the resulting salinity field exhibited substantial ~noise" as compared to the sharp fronts obtained by the previous two grids. Note that the problem is 2D in nature. The noise makes it difficult to evaluate the internal wave speed. It was pointed out earlier that there exists an optimal time step for baroclinic problems, due to the conflicting requirements for time step in ELM and semi-implicit schemes. This is verified in Fig. 4 for the quadrangular grid runs. The optimal step for this grid is about 6 rain. Note the exceptionally large time steps used in this and previous problems, where the Courant number is on the order of one or larger. The robust stability of ELCIRC is crucial in many field scale applications, where the grid size often must be kept small in certain regions. 4. C O N C L U S I O N S

The finite difference nature of ELCIRC requires orthogonal grids, which are often difficult to generate especially with triangles. We showed in this paper that non-uniform orthogonal triangular grids with circumcenters being used as anchoring points for evaluation of elevation and density gradients lead to inferior results as compared with general non-orthogonal grids using centroids as anchoring points. On the other hand, rectangular elements lead to high accuracy and therefore should be used as much as possible. Hybrid grids with rectangles and triangles are most suitable for most applications. The choice of optimal time step seems to be controlled by CFL or similar criterion associated with the baroclinic term, and is related to the maximum internal wave speed of

1799

\\ v

218

\

\

"\

\\

2:> Y, 517

.\

\ ',,,\.

16

//Z

%,,,

ka

\'x,. 15

f

2

3

,

t

4

,

I ........

,

5

1

6

.......

I....

,

.....

7

Time step (minutes)

Figure 4. Numerical errors as a function of time steps used for adjustment under gravity.

the system. Unfortunately this is often difficult to evaluate and a trial-and-error approach may be easier to follow. REFERENCES 1. A.F. Blumberg and G. L. Mellor, A description of a three-dimensional coastal ocean circulation model. Three-Dimensional Coastal Ocean Models. N. Heaps. Washington, D.C., AGU. 4:1-16 (1987). 2. V. Casulli and P. Zanolli, A Three-Dimensional Semi-Implicit Algorithm for Environmental Flows on Unstructured Grids. Institute for Computational Fluid Dynamics Conference on Numerical Methods for Fluid Dynamics VI, 1998. 3. P.J. Frey and P.-L. George, Mesh generation, Hermes Science 2000. 4. D.B. Haidvogel, H. G. Arango, K. Hedstrom, A. Beckmann, P. Malanotte-Rizzoli and A. F. Shchepetkin, Model evaluation experiments in the North Atlantic Basin: simulations in nonlinear terrain-following coordinates. Dyn. Atmos. Oceans 32 (2000) 239-281. 5. D.R. Lynch and A.M. Davies (eds.), Quantitative skill assessment for coastal ocean models, American Geophysical Union, 1995. 6. D.R. Lynch and W. G. Gray, Analytic solutions for computer flow model testing. ASCE Journal of the Hydraulics Division 104(HY10) (1978) 1409-1428. 7. D.R. Lynch, J. T. Ip, C. E. Naimie and F. E. Werner, Comprehensive coastal circulation model with application to the Gulf of Maine. Cont. Shelf Res. 16(7) (1996) 875-906.

1800 8. R.O. Reid, Waterlevel changes. Handbook of Coastal and Ocean Engineering. J. Herbich. Houston, Texas, Gulf Publishing, 1990. 9. ROMS web site: h t t p : / / m a r i n e . r u t g e r s . e d u / p o / t e s t s / g r a v / i n d e x . h t m l 10. L. Umlauf and H. Burchard, A generic length-scale equation for geophysical turbulence models. J. Mar. Res. 6(12) (2003) 235-265. 11. J.J. Westerink, R.A. Luettich, J.K. Wu and R.L. Kolar, The influence of normal flow boundary conditions on spurious modes in finite element solutions to the shallow water equations. Int. J. Num. Meth. Fluids 18 (1994), 1021-1060. 12. Y.-L. Zhang and A.M. Baptista, An efficient adaptive editor for unstructured grid. 7th International Conference on Numerical Grid Generation in Computational Field Simulations, Whistler, B.C., Canada, 395-404, 2000. 13. Y.-L. Zhang and A.M. Baptista, A cross-scale model for 3D baroclinic circulation in estuary-plume-shelf systems: I. Formulation and skill assessment. Cont. Shelf Res. (2004; accepted).