elcom_methods
43 Pages

elcom_methods

Course Number: SITE 2006, Fall 2009

College/University: University of Texas

Word Count: 12859

Rating:

Document Preview

Numerical Techniques in CWR-ELCOM (code release v.1) Ben R. Hodges Centre for Water Research The University of Western Australia Nedlands, Western Australia, AUSTRALIA 6907 CWR manuscript WP 1422 BH March 21, 2000 Numerical Techniques in CWR-ELCOM, March 2000 i copyright 2000 The Centre for Water Research, University of Western Australia Numerical Techniques in CWR-ELCOM, March 2000 ii Acknowledgments...

Unformatted Document Excerpt
Coursehero >> Texas >> University of Texas >> SITE 2006

Course Hero has millions of student submitted documents similar to the one
below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.

Course Hero has millions of student submitted documents similar to the one below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.

Techniques Numerical in CWR-ELCOM (code release v.1) Ben R. Hodges Centre for Water Research The University of Western Australia Nedlands, Western Australia, AUSTRALIA 6907 CWR manuscript WP 1422 BH March 21, 2000 Numerical Techniques in CWR-ELCOM, March 2000 i copyright 2000 The Centre for Water Research, University of Western Australia Numerical Techniques in CWR-ELCOM, March 2000 ii Acknowledgments This research was conducted with the support of the Centre for Environmental Fluid Dynamics (CEFD) at the University of Western Australia. Additional funding and scientific cooperation have been provided by the Western Australia Estuarine Research Foundation, Western Australia Waters and Rivers Commission, Kinneret Limnological Laboratories and the Sydney Catchment Authority. Abstract A description of a 3D estuary and lake computer model (ELCOM) is provided. The model solves the unsteady Reynolds-averaged Navier-Stokes (RANS) equations using a semi-implicit method with quadratic Euler-Lagrange discretization of momentum advection and a conservative ULTIMATE QUICKEST approach for scalar transport. A one-dimensional mixed-layer model is extended to 3D for turbulence closure of vertical Reynolds stress terms. Contents 1 Introduction 1 2 Governing Equations 2.1 2.2 2.3 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Hydrodynamic equations and models Surface thermodynamics . . . . . . . . . . . . . . . . . . . . . . . . 3 3 3 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Computational constraints 3.1 3.2 3.3 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Flooding and drying of grid cells . . . . . . . . . . . . . . . . . . . . . . . . . . . Time step limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 10 11 11 4 Numerical method 4.1 4.2 4.3 4.4 4.5 4.6 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Semi-implicit formulation for momentum . . . . . . . . . . . . . . . . . . . . . . Horizontal diffusion discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . Baroclinic discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Free surface discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Advective discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6.1 4.6.2 4.6.3 Hybrid momentum advection scheme . . . . . . . . . . . . . . . . . . . . . Linear semi-Lagrangian methods . . . . . . . . . . . . . . . . . . . . . . . Quadratic Euler-Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . ii 14 14 14 16 16 17 17 17 18 19 Numerical Techniques in CWR-ELCOM, March 2000 iii 4.7 Scalar transport 4.7.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 23 24 24 25 28 28 29 30 Scalar horizontal diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Mixing and vertical diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8.1 4.8.2 4.8.3 4.8.4 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mixing energy budgets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Discrete mixing operator . . . . . . . . . . . . . . . . . . . . . . . . . . . Advantages and limitations of the mixing model . . . . . . . . . . . . . . 4.9 Wind momentum model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Sidewall and bottom boundary conditions . . . . . . . . . . . . . . . . . . . . . . 5 Future work 32 Chapter 1 Introduction While the terms "simulation" and "model" are sometimes used interchangeably in the literature, we would like to distinguish between the two concepts applied to numerical solution of temporal evolution equations. The former is taken to be the solution to a set of equations based on first principles, while the latter is the solution of a set of equations that combines first principles with empirical observations and/or scaling laws. While the distinction may seem only of pedagogical importance, it does play a role in how a numerical scheme is developed, tested and applied. The Centre for Water Research Estuary and Lake Computer Model (ELCOM) described in this report is designed as a model of physical processes rather than a simulation of fluid flow. Simulation methods might be considered to encompass Direct Navier-Stokes (DNS) and Large-Eddy Simulations (LES) techniques and a subset of the Reynolds-Averaged Navier-Stokes (RANS) methods presented in the literature. A simulation method will solve the evolution equations to the designed degree of accuracy given a sufficiently fine grid. The accepted test of "sufficient" resolution requires a grid resolution study with at least one numerical simulation performed with twice the resolution as the base grid scale (i.e. the grid used for analysis of simulation results). This typically requires an increase in computational time and memory of at least one order of magnitude. If the simulation results change dramatically with the grid resolution, then the base grid is considered to be under-resolved. Since maximum practical resolution is a function of computational power, the accepted procedure in numerical simulation is to reduce the Reynolds and/or Froude numbers of the flow under consideration until a wellresolved flow can be simulated with the available computer. In effect, simulation methods require re-scaling the problem to fit the computer. The fundamental tenet of flow simulation is that a well-resolved solution to the evolution equations is invariant to further refinement of the grid or change of numerical method. In contrast, three-dimensional (3D) models of geophysical flows with topographic effects are often not suitable for rescaling without distorting the physics that are under investigation. If the flows cannot be rescaled, the attainable resolution is strictly a function of the available computational power. It is rare to see a numerical model of a geophysical flow over complex topography with a grid resolution study most models are run with all available computational power at the maximum possible grid resolution. Grid refinement studies are simply not 1 Numerical Techniques in CWR-ELCOM, March 2000 2 practical for most investigators. One might suspect that many of these models are, strictly speaking, "under-resolved" in the solution of the evolution equations. It might be somewhat less derogatory to label this practice "coarse-grid" modeling. Solutions obtained with 3D coarsegrid models may be highly dependent on the grid resolution and the numerical method. The purist might argue that such solutions are simply nonsense and should be ignored. However, carefully conducted coarse-grid 3D models are a necessary bridge between highly-resolved threedimensional simulation of idealized flows (presently obtainable through DNS, LES and RANS for limited Reynolds and Froude numbers), and 1D models of lake dynamics and 2D models of estuaries. Indeed, it is our contention that 3D modeling over complex topography will remain the province of coarse-grid numerical models for the forseeable future. The challenge of coarse-grid modeling is to develop numerical methods that are sufficiently accurate and robust for attainable grid resolutions while providing reasonable representations of the physics. That better results might be attained with finer grids or a higher-order numerical method cannot be allowed to deter the use of methods that are sufficiently accurate methods for the purposes at hand. It must be held in mind that a geophysical numerical model is often forced with field data of an uncertain quality. In particular, available wind data for lakes and estuaries is often limited to a single value of wind speed and direction that may be located some distance from the water. Thus, an extremely accurate solution of the flow physics on a fine grid may not provide a more accurate model of the physics of a system given the uncertainty in the forcing data. Our objective is to produce a numerical model that can reproduce the first-order 3D physical response of a lake to environmental forcing on a coarse grid with low CPU time requirements so that the method can be used to develop greater understanding of the variables which drive mixing dynamics and scalar transport in geophysical systems. For coarse-grid modeling, the evolution equations derived from first principles (as used in simulations) must be carefully reconsidered in conjunction with the behaviour of the numerical methods applied. The a priori decision to use coarse-grid methods for a model of a geophysical flow might seem to imply that using numerical schemes higher than 1st order in accuracy would be, in the words of J. Koseff, "measuring with a micrometer after cutting with an axe". However, for coarse-grid modeling, the character of the truncation error of the numerical methods is may be more important than the magnitude of the error or the formal order of the scheme. As a result, suitable methods may be required to be higher order than one might first expect. In ELCOM, the advancement in time for free-surface evolution is a first-order backwards-Euler method with second-order discretization used for most spatial derivatives. Scalar transport, however, uses a conservative third-order scheme to minimize unphysical oscillations, numerical diffusion, and non-monotonicity in the scalar fields that are a feature of lower-order schemes. Chapter 2 Governing Equations 2.1 Introduction The governing equations and fundamental models used for three-dimensional transport and surface thermodynamics in ELCOM are summarized in figures 2.1 through 2.4. The transport equations are the unsteady Reynolds-averaged Navier-Stokes (RANS) and scalar transport equations using the Boussinesq approximation and neglecting the non-hydrostatic pressure terms. The free surface evolution is governed by and evolution equation developed by a vertical integration of the continuity equation applied to the Reynolds-averaged kinematic boundary condition. Surface thermodynamics are governed by bulk transfer models. 2.2 Hydrodynamic equations and models The unsteady RANS equations are developed by filtering the unsteady Navier-Stokes equations over a time period that is long relative to sub-grid scale processes, but small relative to the unsteady grid-scale processes that are of interest. In an unsteady RANS numerical method the time scale of the averaging is the time step used in advancement of the evolution equations. Thus, the maximum time step for a given grid resolution is fundamentally limited by the grid scale physics, regardless of the numerical method. This is an important point that is lost when "unconditionally stable" numerical methods are discussed: the stability of the method at a large time step is irrelevant if the size of the time step implies a temporal averaging scale that encompasses grid-scale motions. In most numerical approaches, the sub-time scale nonlinear terms (ui uj )are represented by an eddy viscosity (jk ) defined in a tensor form as: jk ui ui jk - ui uj xk xk (2.19) where is the molecular viscosity, uk are the sub-time-scale velocity fluctuations, and the overbar indicates the Reynolds-averaging filter. The off-diagonal eddy-viscosity terms are generally set 3 Numerical Techniques in CWR-ELCOM, March 2000 4 Hydrodynamic equations Transport of momentum: U t + Uj + Continuity Uj = 0 xj Boundary conditions on momentum free surface: (2.2) U 1 = -g + xj x 0 x U U 1 + 2 x1 x1 x2 x2 dz z (2.1) 3 U x3 - f U + x3 U = 0 x3 Ui = 0 (2.3) bottom and sides: (2.4) Transport of scalars: C + (C Uj ) = t xj x1 Boundary conditions on scalars C = 0 xj Free-surface evolution = - t x (u ) = C10 Momentum input by wind 2 1 C x1 + x2 2 C x2 = x3 3 C x3 +S (2.5) (2.6) U dz b (2.7) Free-surface wind shear 1 (air) (W W ) 2 W (water) (2.8) (u ) dU = dt h 2 (2.9) Figure 2.1: Numerical Techniques in CWR-ELCOM, March 2000 5 Nomeclature in hydrodynamic equations i, j, k, m , U 0 g f h ij three component space (e.g. i = 1, 2, 3) horizontal two component space (e.g. = 1, 2) Reynolds-averaged velocity Reynolds-averaged free surface height reference density density anomaly (i.e. the difference between the in situ density and the reference density) gravitational constant Coriolis constant height of wind-mixed layer Kronecker delta (0 for i = j, 1 for i = j) two-component permutation tensor molecular viscosity eddy-viscosity tensor scalar source term turbulent Schmidt number (or Prandtl number for temperature) scalar concentration x3 value at bottom of domain vector wind speed in direction bulk wind stress coefficient for wind values at 10 meters wind shear velocity in direction Figure 2.2: km Q Sc c b W C10 (u ) Numerical Techniques in CWR-ELCOM, March 2000 6 Surface thermodynamics equations Temperature/internal energy relation T = Longwave radiation emitted: QR(emitted) = - absorbed: QR(absorbed) = Evaporative heat loss QW = L Sensible heat flux QH = CH cp(air) (air) T(air 10) - T(water) Short wave radiation penetrating free surface (d = 0): 2 QS (0) = Q(surf ace) 1 - 0.65C(cloud) (1 - Rsw ) (air) 2 1 + 0.17C(cloud) (water) Q (water) V cp (2.10) 273.2 + T(water) 4 (2.11) 273.2 + T(air 2) 4 (1 - Rlw ) (2.12) 0.622 CW u(wind) (air) e(sat) (Rh - 1) P(atm) (2.13) (2.14) (2.15) distribution with depth ( - z): QS (d) = QS (0) exp {-e d} Emissivity of air (air) (2.16) = CE 273.2 + T(air 2) 2 (2.17) Saturation vapor pressure log10 e(sat) = 0.7859 + 0.03477T -3 1 + 0.00412T (2.18) Figure 2.3: Numerical Techniques in CWR-ELCOM, March 2000 7 Nomenclature for surface thermodynamics equations T QR QW QH QS Q(surf ace) temperature longwave radiant heat flux evaporative heat flux sensible heat flux shortwave radiant heat flux shortwave radiation reaching water surface emissivity d e(sat) (air) C(cloud) L Rlw Rsw Rh CH CW u(wind) cp depth (i.e. - z) saturated vapor pressure air density fractional cloud cover latent heat of evaporation Stefan-Boltzmann constant reflectivity of surface to incoming longwave radiation reflectivity of surface to incoming shortwave radiation relative humidity bulk transfer coefficient for sensible heat bulk transfer coefficient for evaporative heat wind speed specific heat at constant pressure Figure 2.4: Numerical Techniques in CWR-ELCOM, March 2000 8 to zero and the molecular viscosity is ignored so that we obtain: 1 ui = ui u1 x1 ; 2 ui = ui u2 x2 ; 3 ui = ui u3 x3 (2.20) In ELCOM, eddy-viscosity is used to represent the horizontal turbulence closure. In the vertical direction, ELCOM can apply either a vertical eddy viscosity or a mixed-layer model discussed in 4.8. The free-surface evolution is governed by vertical integration of the continuity equation for incompressible flow from the bottom (b) of the water column to the free surface () applied to the kinematic boundary condition (e.g. Kowalik and Murty, 1993). To derive the Reynoldsaveraged free-surface evolution, equation (2.7) the Reynolds-averaging filter must be applied to the kinematic boundary condition, resulting in = U3 - U - u h t x x (2.21) Note that the sub-time scale correlation term (u h ) has not been discussed in the literature as the Reynolds-averaging filter has not been rigorously applied to boundary conditions; the equivalent term in large-eddy simulation is discussed in Hodges and Street (1999). The modeling of this term remains an unexplored area of free-surface flows; thus we are forced to assume that it is either small or only of local importance. In the former case it is readily neglected. In the latter case, the present use of implicit discretization of the free-surface solution is numerically diffusive at high barotropic CFL conditions, so local barotropic effects are not well-resolved. Thus, neglecting the subtime-scale nonlinear effects is justified by our focus on basin-scale baroclinic motions. However, we can speculate that the nonlinear correlation u h, might be both large and non-local where nonlinear sub-basin-scale surface waves enter shallow or restricted waters. The neglect of this term is predicated on the assumption that sub-time scale nonlinear effects do not significantly effect the local mean free-surface height. In general, this is a good assumption for the resolution typically achieved in simulations of lakes and estuaries, and will be used herein. The free surface boundary condition is an approximate form of the dynamic boundary condition that neglects: dynamic pressure, local horizontal variation of the wind, updrafts, downdrafts and surface tension. This is appropriate for lake simulations using the hydrostatic approximation. 2.3 Surface thermodynamics Heat exchange through the water's surface in ELCOM is governed by standard bulk transfer models found in the literature (e.g. Amorocho and DeVries, 1980; Imberger and Patterson, 1981; Jacquet, 1983) and provided in Figures 2.3 and 2.4. The energy transfer across the free surface is separated into non-penetrative components of longwave radiation, sensible heat transfer and evaporative heat loss, complemented by penetrative shortwave radiation. Non-penetrative effects are introduced as sources of temperature in the surface-mixed layer, while penetrative effects are introduced as source terms in one or more grid layers based on and exponential decay and an extinction coefficient (Beer's law). Numerical Techniques in CWR-ELCOM, March 2000 9 The simplest approach to introducing non-penetrative effects is to assume that all the surface heat transfer is absorbed in the uppermost grid cell in the domain. However, this approach results in the effective model of the surface heat transfer changing as a function of the grid resolution at the surface. In ELCOM, an approach that is grid-independent for vertical grid resolution of one meter or less is employed. The surface heat transfer is modeled as occurring over the first one meter below the free surace with an exponential decay such that 98% of the surface heat transfer is absorbed in this region. Naturally, if the vertical grid resolution near the free surface is coarser than 1 m, then the entire surface heat transfer occurs into the upper grid cell. In shallow or very clear waters, shortwave radiation may penetrate to the bottom of the lake or estuary. There is an open question as to how best to model the radiant energy that reaches the bottom. A comprehensive heat budget model would require: (1) absorption and reflection of the short wave radiation by the sediments, (2) long wave radiation emission from the sediments, and (3) conduction and convection models at the bottom boundary. Data to develop and validate a comprehensive model is presently not available. As a simpler approach, we will consider that any short wave radiation that reaches the bottom boundary is treated by a model similar to that used for the surface heat transfer near the free surface. We assume the shortwave radiation is converted into longwave radiation and/or sensible heat transfer that is absorbed in the first meter above the bottom using exponential decay. Further details on the ELCOM thermodynamic modeling can be found in a CWR technical report (Hodges, 2000). Chapter 3 Computational constraints 3.1 Introduction The ELCOM numerical method takes its basic structure from the TRIM scheme of Casulli and Cheng (1992) with adaptations for accuracy, scalar conservation, numerical diffusion and implementation of a mixed-layer model. Other adaptations of TRIM can be found in Casulli and Cattani (1994); Casulli (1997); Gross et al. (1998, 1999). Hereinafter these are collectively referred to as the TRIM method, with distinctions made between adaptations only where relevant. The methods in ELCOM are extended beyond TRIM by including: (1) a hybrid advection scheme for momentum; (2) an energy-based mixing model for vertical diffusion; and (3) conservative advection of scalars using an third-order explicit method, in contrast to the implicit scheme in the adaptation of TRIM by Gross et al. (1998). The solution grid uses rectangular Cartesian cells with fixed x and y (horizontal) grid spacing while the vertical z spacing may vary as a function of z but is horizontally uniform. The uniform horizontal spacing and rectangular cells enables application of a simple, efficient finite-difference/finite-volume scheme on a staggered grid. Vertically varying layer thickness allows grid layers to be concentrated where the greatest resolution is required. The free surface height in each column of grid cells is computed from the free surface evolution equation and is allowed to move up and down through the grid layers as required. The grid stencil is the Arakawa C-grid: velocities are defined on cell faces with the free-surface height and scalar concentrations on cell centers. The free-surface height in each column of grid cells moves vertically through grid layers as required by the free-surface evolution equation. The numerical scheme for computing the temporal evolution of velocity is a semi-implicit solution of the hydrostatic Navier-Stokes equations using a hybrid discretization of advection terms. The hybrid scheme locally alters the numerical discretization method used for the advective terms as a function of the local CFL condition. Semi-implict solution of the hydrostatic Navier-Stokes equations involves discretizing the free surface evolution in an implicit manner while retaining explicit source terms for advection, horizontal diffusion, and baroclinic terms (e.g. Casulli and Cheng, 1992; Kowalik and Murty, 1993). 10 Numerical Techniques in CWR-ELCOM, March 2000 11 3.2 Flooding and drying of grid cells In solving the Navier-Stokes and transport equations with a free surface on a fixed grid, the discretization of cells that may be "wet" at one time step and "dry" at another presents problems for semi-implicit methods (i.e. any method that solves the free surface evolution in an implicit formulation with a mix of explicit and implicit discretizations for advective, diffusive and baroclinic terms). With fully-explicit methods, all terms are discretized in a consistent manner for the time `n' solution space. Similarly, a consistent fully-implicit method can be written (although its solution is more complicated) using all terms discretized in the time `n + 1' solution space. However, for a semi-implicit technique on a fixed-grid with a moving free-surface, there is a fundamental inconsistency between the time `n' and time `n + 1' solution spaces; i.e. they may not be exactly equal, so some free-surface cells may have only time `n' values while other cells may have only time `n + 1' values. The distribution of these cells is unknown a priori. It follows that the entire solution space of the governing equations for a free-surface flow cannot be formally discretized in a consistent manner using a semi-implicit technique and a fixed grid. The inconsistency can be eliminated by ensuring there is a sufficiently thick upper grid cell that encapsulates all possible free surface motion. Unfortunately, many lakes undergo changes in surface level that make it impractical to specify limits on the surface position. The TRIM method side-steps this conundrum by discretizing the momentum equations against the background of the time `n' free surface. Using this approach, numerical solution of the semi-implicit equations is not attempted in a cell that exists only in time `n + 1' space. Thus, newly wet cells must be initialized with velocity and scalar data from surrounding cells that were previously wet. Newly dry cells must be cleared of data since the solution in these cells (while calculated) is meaningless. The approach in TRIM has been adopted in the present work as it provides a simple framework around which to build an efficient semi-implicit method. We find it convenient to define our solution space for the advance from time `n' to time `n + 1' as the intersection of the set of (i, j, k) cells at time `n' with the set of (i, j, k) cells at time `n + 1'. Hereafter this will be referred to as the `n ' set of cells. 3.3 Time step limitations Both TRIM and ELCOM are unconditionally stable for purely barotropic flows; that is, they will produce stable numerical results for any size of time step. However, for stratified flows, both methods use explicit discretization of the baroclinic terms in Equation 2.1 leading to a time step constraint based on the internal wave Courant-Friedrichs-Lewy condition (CFLb ) such that CFLb < 2 is required, or 1 t < 2 (3.1) (g D) 2 x The left-hand side is defined as the baroclinic CFL number (CFLb ), where g is the reduced gravity due to stratification (g-1 ), the effective depth is D, and g D is an approximation 0 of the wave speed of an internal wave. This baroclinic stability condition is generally the most restrictive condition in a density-stratified flow. The importance of the internal wave speed can Numerical Techniques in CWR-ELCOM, March 2000 12 be seen in a simple scaling analysis: Typical lake stratifications provide a modified gravity of g O 10-2 ms-2 while lake depths (D) are O (10) m to O 102 m. Internal waves propagate 1/2 at C (g D) , giving wave speeds of O (1) ms-1 . Maximum horizontal water velocities in a lake are typically O 10-1 ms-1 , and the desired horizontal grid size is O 102 m. The maximum allowable time step for a limiting CFL condition is: t < CFL x U (3.2) If U is taken as the horizontal water velocity, the maximum allowable time step for a limiting CFL of one is O 103 s. If U is taken as the internal wave speed, the maximum allowable time step for a limiting CFL of one is O 102 s. Thus, in the horizontal direction, the internal wave speed rather than the flow velocity controls the maximum allowable time step and high CFL capability in the horizontal direction is generally unnecessary. However, in the vertical direction, there is a definite advantage to a numerical scheme that is stable for CFL > 1. Practical vertical grid resolutions are typically O 10-1 m to O (10) m depending on the lake morphology and available computational power. Internal wave motions can produce vertical velocities of O 10-2 ms-1 , so fine grid resolutions with a CFL limit of one can result in an unacceptable time step limitations of O (10) s. Scalar transport in ELCOM uses an explicit approach has an advective Courant-FriedrichsLewy condition (CFLa ) such that u tx-1 < 1 is required. This condition does not effect the time step of the momentum solution (t), but instead is used to compute sub-time steps (t) for the scalar transport solution. The implicit scheme of Gross et al. (1998) does not have CFL limitation for vertical velocity and is therefore preferable where fine grid resolution is used in the vertical direction and high vertical CFL values are expected. At coarse vertical grid resolutions (low vertical CFL values) the present explicit approach should prove more computationally efficient. A final stability constraint for semi-implicit schemes with explicit horizontal diffusion (e.g. TRIM and ELCOM) is the viscous stability condition (derived for homogenous flows in Casulli and Cattani, 1994); x2 y 2 t (3.3) 2 (x2 + y 2 ) This is typically at least an order of magnitude less restrictive than the baroclinic stability condition. A constraint that becomes important when using large time steps in a geophysical model is whether the velocity field can be considered Lipschitz at the grid and time scales applied: i.e. is the field sufficiently smooth for a numerical approximation (Iserles, 1996). In Smolarkiewicz and Pudykiewicz (1992) the numerical Lipschitz constant B is defined with the condition that it must be less than unity: v t < 1 (3.4) B = x This was demonstrated to be a necessary condition in a multi-time-level semi-Lagrangian method to prevent trajectories from intersecting as they are traced back in time and space. While the present Euler-Lagrange method tracks trajectories back only in space, the above Lipschitz Numerical Techniques in CWR-ELCOM, March 2000 13 condition heuristically applies as a fundamental statement of the necessary behavior of the velocity field to allow reasonable approximation by a numerical model. As demonstrated in Smolarkiewicz and Pudykiewicz (1992), a numerical method may remain stable at a high B, but the results will not be accurate. This places a fundamental limit on the allowable time step that may be reasonably used in a model. The importance of this point is that the maximum time step that can be used may be a function of the physics, rather than the stability of the numerical method. Chapter 4 Numerical method 4.1 Introduction The governing equations are discretized on a Cartesian solution grid in a staggered formulation where single velocity components are defined on each face and scalars are defined at the cell centers. In discrete equations, the cell faces are represented by subscripts such as i + 1/2, while the centers are represented with integer (i, j, k) values. The notation of Casulli and Cheng is used in the following description. For the discrete form of equations, we will use subscripts to represent the position in discrete (i, j, k) space. Let Un+1 represent the water column vector of i,j time `n + 1' velocity values at position (i, j), that exist in the time n solution space for all k that satisfy bi,j kmax m=1 n zi,j,m i,j (4.1) where bi,j is the height of the bottom of the domain at point (i, j), i,j is the height of the free surface, and kmax is the maximum number of grid cells in the vertical direction. Similar definitions are applied for other vector quantities in this report. 4.2 Semi-implicit formulation for momentum ELCOM is designed to operate with various solution options in the momentum discretization. The fundamental semi-implicit evolution of the velocity field can be discretized on the n solution space in a manner similar to the TRIM approach as: Un+1,j = A-11 ,j Gn 1 ,j - g i+ i+ 1 i+ 2 2 2 t n+1 n+1 1 i+1,j - i,j x t n+1 n+1 1 i,j+1 - i,j y n n + (1 - 1 ) i+1,j - i,j n n + (1 - 1 ) i,j+1 - i,j (4.2) (4.3) n+1 Vi,j+ 1 = A-1 1 Gn 1 - g i,j+ i,j+ 2 2 2 14 Numerical Techniques in CWR-ELCOM, March 2000 15 where G is an explicit source term vector, 1 is the "implicitness" of the free surface discretization1 . The default semi-implicit scheme in ELCOM is a backwards-Euler discretization (i.e. = 1) of the free surface evolution that is formally 1st order accurate in time. It has been demonstrated (Casulli and Cattani, 1994) that the backwards Euler method for solution of the hydrostatic momentum equations can be extended to the general two-level scheme of Equations (4.2) and (4.3) that is formally second-order accurate (when = 0.5). However, in coarse grid modeling, this increase in numerical accuracy does not always result in an increase in model skill. In general, many lake an estuary simulations are conducted where the barotropic mode will be solved with CFL conditions that may range from 5 to 10 or greater. Under these conditions, semi-implicit discretizations may be stable, but the "fidelity" of the representation of the flow physics is a function of the type of flow under consideration The character of the truncation error is critical to understanding the performance of the method. With the 1st order method, the lead error term is 2nd order and results in damping of waves on the free surface. With the 2nd order method, the leading error term is dispersive and results in numerical waves on the free surface which propagate accross the domain; typically causing a linear barotropic wave to evolve into a steep-fronted bore, causing high velocities in localized areas as the surface wave is affected by topography. Thus, the first order method results in good representation of the shape of the free surface and the local barotropic velocities, but shows excessive damping of the inertial response of the free surface. In contrast, the second order method maintains the energy in the surface wave form with minimal numerical dissipation, but has a poor representation of the wave form. In a hydrostatic solution, the dispersive waves cause spurious local forcing of throughout the water column and are detrimental to the skill of the solution. In contrast, the excessive damping of the surface wave causes a decrease in the large scale motions associated with the barotropic response when the wind relaxes, i.e. the "ringing" of the barotropic mode is damped. In general, a strongly forced system is better modeled with the backwards Euler scheme as the wave energy beyond two or three periods is often irrelevant to the first-order physics. Using an any two-level implicit discretization (e.g. Casulli and Cheng, 1992), or any explicit discretization technique, the A matrix can be represented as: bn + k cn 0 0 0 an-1 bn-1 cn-1 0 0 bn-2 cn-2 0 0 an-2 A = . . . . . . . . . . . . . . . b2 c2 0 0 a2 0 0 0 a1 b1 + 1 The terms in A are set by boundary conditions (see 4.10), while the a, b, and c terms are: bk = -ak + zk - ck ak = -2 3 t z k+ 1 2 (4.4) (4.5) 1 The generalized implicitness option for 0.5 1.0 is coded in ELCOM version 1, but has not been fully 1 tested. Numerical Techniques in CWR-ELCOM, March 2000 16 ck = -2 3 t z (4.6) k- 1 2 The 2 coefficient is determined by the choice of numerical discretization techniques. For 2 = 1 the vertical viscous term is discretized using a backwards Euler technique. For 2 = 0, the discretization is explicit and A is zero for all terms off the main diagonal. ELCOM selects 2 = 1 when the user chooses to model vertical turbulent diffusion with an eddy-viscosity. For the mixed-layer model, ELCOM sets 2 = 0. The source terms G in Equations (4.2) and (4.3) can be represented as: ~ Gn 1 ,j = L Ui+ 1 ,j i+ 2 2 ~ - t Bn 1 ,j + Dx U i+ 2 i+ 1 ,j 2 ~ + Dy U ~ + Dy V i+ 1 ,j 2 ~ - f Vi+ 1 ,j 2 ~ + f Ui,j+ 1 2 (4.7) (4.8) ~ Gn 1 = L Ui,j+ 1 i,j+ 2 2 ~ - t Bn 1 + Dx V i,j+ 2 i,j+ 1 2 i,j+ 1 2 The L ( ) operator represents advective discretization (see 4.6) , B() represents baroclinic discretization (4.4), D() represents horizontal turbulent diffusion discretization (4.3). ELCOM allows vertical diffusion to be computed either using the eddy viscosity terms or a vertical mixing ~ model. When the eddy-viscosity form is used, U = U n , while the mixing model is represented as an operator such that: n ~ Ui,j,k = M Ui,j,k (4.9) The mixing operator M () is described in detail in 4.8. 4.3 Horizontal diffusion discretization While the original TRIM approach (Casulli and Cheng, 1992) applied the discretization of the horizontal diffusive terms at the pathline origin, the additional complexity was not found to provide any significant advantage in the accuracy of the present method. Therefore, the horizontal diffusion terms (Dx , Dy ) from Equations 4.7 and 4.8 are discretized using a secondorder stencil such that: Dx n i,j,k = n - 2n + n i,j,k i-1,j,k x2 i+1,j,k (4.10) 4.4 Baroclinic discretization The baroclinic term B in the x direction is discretized as: Bn 1 ,j,k i+ 2 g = 0 x F F n i+1,j,m - m=k m=k i,j,m (4.11) where k = F is the cell containing the free surface. Similar expressions for diffusion and baroclinic terms are obtained for the y direction. The elimination of the vertical diffusion terms in Equations 2.1 and 2.5 allows the present scheme to dispense with the tridiagonal matrix inversion for each horizontal velocity component and transported scalar required for each (i, j) water column in the TRIM scheme. Numerical Techniques in CWR-ELCOM, March 2000 17 4.5 Free surface discretization The free surface evolution is governed by the solution of Equation 2.7, which can be discretized as: n+1 n i,j = i,j - 1 t t T T x (Zn ) Un+1 - y (Zn ) Vn+1 x y t t T T x (Zn ) Un - y (Zn ) Vn - (1 - 1 ) x y (4.12) where Z is a vector of the vertical grid spacing, and the operators x and y indicate discrete differences, e.g. x () i+ 1 - i- 1 (4.13) 2 2 Substitution of equations (4.2) and (4.3) into equation (4.12) provides a pentadiagonal system of equations for the time `n + 1' free surface height that is readily solved using a conjugate gradient method. The approach adopted for conjugate gradient solution in the present work is similar to TRIM, as detailed in Casulli and Cheng (1992). 4.6 Advective discretization ELCOM has several different advective discretization approaches available: quadratic or linear Euler-Lagrange, upwind, QUICKEST, centered finite differences and a hybrid scheme2 . 4.6.1 Hybrid momentum advection scheme One of the difficulties of numerical modelling at geophysical scales is the wide range of flow conditions that may be present. In particular, internal waves may intermittently produce strong vertical motions over relatively small regions. Where resolution of internal waves is deemed important, the choice of numerical method is driven by the need for accuracy and stability in a small portion of the overall flow field. While many explicit spatial discretization methods are stable for CFL < O (1), their accuracy in 3D computations for CFL > O (0.5) is generally poor when the flow direction is not aligned with the grid. The drawback of the quadratic method is that it is computationally expensive. However, for low CFL regions (CFL < 0.1), the solution of a quadratic semi-Lagrangian method is dominated by terms in the same seven-point upwind stencil that is obtained by using quadratic upwind discretization. This similarity can be exploited to reduce the computational requirement in regions of low CFL without significantly reducing the accuracy of the overall solution method. The concept of applying different schemes in different regions can be generalized into the concept of a "hybrid" numerical method. A general hybrid method is one where different 2 Not all schemes have been thoroughly tested in ELCOM version 1. It is recommended that users apply the quadratic Euler-Lagrange or the hybrid scheme. Numerical Techniques in CWR-ELCOM, March 2000 18 solution schemes are applied in different flow regions based upon a criteria measured in the flow field. For our purposes, the criteria is the CFL, and we apply one discretization technique for low CFL, and a different technique for high CFL. The present hybrid method has been tested to two levels, using quadratic semi-Lagrangian discretization for regions where (0 < CF L < 2). The upper limit cutoff is the maximum CFL that can be computed using a quadratic SL method without requiring the additional computational expense of adjusting the stencil region. As a practical matter, the accuracy of computations for CF L > 2 is questionable, so the upper limit is not an unreasonable requirement. For regions with CF L > 2, the default ELCOM method applies linear semi-Lagrangian discretization (i.e. the approach used in TRIM) to minimize the computational effort of repositioning the stencil. ELCOM has the capability3 to use a third discretization method (e.g. QUICK) for very low CF L regions (i.e. CF L < 0.1). 4.6.2 Linear semi-Lagrangian methods For a stratified flow, semi-Lagrangian methods are desirable in that they are both accurate and stable in the regime where 0.1 < CFL < 1. As a further advantage, they retain their stability in higher CFL regimes, although their accuracy suffers at high CFL values unless the flow streamlines are well resolved by the grid. However, in stratified flows, the ability to provide stable solutions for CFL > 1 in the horizontal direction is of secondary importance in the selection of a numerical method. The maximum time step is generally limited by either by (1) the baroclinic wave propagation speed in the region of strongest stratification, or (2) the maximum acceptable numerical diffusion in scalar transport. The former requires wave CFL < 1 for stability, while the latter limitation is a matter of grid resolution and the length of time over which model results are required. The semi-Lagrangian form of advection is obtained by finding the approximate point in continuous space (the "Lagrange" point) which would be advected to a discrete point (i, j, k) by the velocity field (U, V, W ) over the time step t. That is, the particle position (i, j, k) is numerically marched back along the streamlines represented by the velocity field U, V, W . The U, V, W field may be obtained from single or multiple time levels, depending upon the accuracy and computational complexity desired. A detailed review of semi-Lagrangian techniques is found in Staniforth and C^t (1991). In a linear, single-time level semi-Lagrangian method, the oe Lagrange point is found using x = -U t y = -V t z = -W t (4.14) (4.15) (4.16) The value of a variable at the Lagrange point is found using a trilinear upwind interpolation: L = x y z i-1,j-1,k-1 + x y z 1- x x y z i,j-1,k-1 y z 3 Three-level hybridization is not fully tested in ELCOM version 1. Numerical Techniques in CWR-ELCOM, March 2000 19 + + + x x x x z x y y z i-1,j,k-1 + 1- i-1,j-1,k y z x y z y x y z z 1- 1- i-1,j,k + 1 - 1- y z x y z z x x y y z 1- i,j,k-1 + 1 - 1- 1- 1- x y z x y z 1- i,j,k (4.17) The above constitutes an eight-point stencil for a semi-Lagrangian methods with linear interpolation. As the flow field approaches a uniform 1D flow, the linear semi-Lagrangian technique reduces to a linear upwind method for CF L 1. < Indeed, it is perhaps easier to think of the semi-Lagrangian method as a 3D linear upwind stencil that is successively repositioned for high CFL conditions. It follows that linear semi-Lagrange method exhibits the unacceptable levels of numerical diffusion that are characteristic of linear upwind methods. This can be ameliorated by using a quadratic method for interpolation of values at the Lagrange point (see 4.6.3) . 4.6.3 Quadratic Euler-Lagrange The quadratic semi-Lagrange method extends the 8-point upwind stencil with trilinear interpolation (used in TRIM) to a 27-point upwind stencil using quadratic Lagrange polynomial interpolation. In the notation of Casulli and Cheng (1992), the operator interpolates the velocity field (on the x face) to the position i + 1 - a, j - b, k - d , where a, b and d are real 2 numbers that represent displacements from the i + 1 , j, k position. This is an estimate of the 2 origination of the pathline of a particle moving through the time n velocity field that ends at the position i + 1 , j, k after time t. A similar notation is used for the y faces at j + 1 . It 2 2 is convenient to denote the pathline origination point by a superscript (p) so that the advective operators L() in Equations (4.2) and (4.3): ~ L Ui+ 1 ,j 2 ~ ~ = Ui+ 1 -a, j-b, k-d = U (p) 2 (4.18) The process of computing the pathline origin in 2D for bilinear interpolation is discussed in Casulli (1990) and is illustrated in Figure 4.1 for 2D quadratic Lagrange interpolation. For Lagrange interpolation on non-uniform grids it is convenient to consider the physical space coordinates (x, y, z) that correspond to the (i, j, k) computational indices. Three-dimensional quadratic Lagrange interpolation to find the value of U (p) at any point x(p) , y (p) , z (p) is computed in three sweeps illustrated in Figure 4.2, requiring 9 vertical interpolations of the form ^ Ui+, j+ = L0 j+ Ui+, j+, k + L1 j+ Ui+, j+, k1 + L2 j+ Ui+, j+, k2 (4.19) i+, i+, i+, where the (i, j) subscripts denote the position of the vertical line to be interpolated, while and are {0, 1, 2}, with the sign determined by the upwind direction of the stencil (as is the sign for the k increment in the U position subscripts). The Lagrange polynomial coefficients, Lm for each (i, j) line are computed from the standard Lagrange coefficient formula (e.g. Al-Khafaji and Tooley, 1986) 2 z (p) - zk Lm = : m = 0, 1, 2 (4.20) zkm - zk =0 =m Numerical Techniques in CWR-ELCOM, March 2000 20 B A -VA dt C -UA dt D b y x a Figure 4.1: Two-dimensional illustration of Euler-Lagrange streamline computation using quadratic Lagrange interpolation. Velocity vector A, with components UA and VA , is used to track particle path from position (i, j) back to the base of vector B using the momentum sub-time step dt. Velocity vector B is computed from the nine grid nodes upwind of the velocity vector at position (i, j). Vector B is used to track the particle path back to the base of velocity vector C, which is again interpolated from the surrounding nine grid nodes. This is repeated n times where n dt = t. If a vector base position is not contained within the nine upwind grid nodes, the upwind stencil must be repositioned. In the present code, linear interpolation is used for the rare instances when repositioning is necessary. The final vector is the result of the Euler-Lagrange operator in Equations 4.7 and 4.8, i.e. L (Ui,j ) = Ui-a, j-b . The number of sub-time steps (n) may be set arbitrarily, with high values providing greater accuracy and higher computational cost. Note that n = 1 everywhere corresponds to quadratic upwind discretization 1. As a general rule, the minimum n is set as and has poor accuracy characteristics unless the CFLa a local function of the grid and flow field such that U dt x-1 < 1 Numerical Techniques in CWR-ELCOM, March 2000 21 (i-2, j, k) (i,j,k) (i-2, j-2, k) d z y x b (i-2, j-2, k-2) a (i, j-2, k-2) interpolation points in k (9 lines) interpolation points in j (3 lines) interpolation points in i (1 line) pathline origin {x(p), y(p), z(p)} (i, j, k-2) Figure 4.2: Figure 4.2. Schematic of three-dimensional quadratic Lagrange interpolation with successive interpolations in the k, j, then i directions. For clarity, this illustration shows interpolation for a uniform grid, however the method can be applied to non-uniform grids without any further adaptation. where z p is the vertical coordinate of the interpolated point and the sign is determined to obtain an upwind stencil. The vertical interpolations are followed by 3 horizontal interpolations in the y direction of the form ^ ^ ^ Ui+ = L0 j Ui+, j + L1 j1 Ui+, j1 + L2 j2 Ui+, j2 i+, i+, i+, Finally, a single interpolation in the x direction is applied as U (p) = L1 Ui + L2 Ui1 + L3 Ui2 (4.22) (4.21) The Lagrange coefficients in Equations 4.21 and 4.22 are computed using Equation 4.20 with y or x substituted for z as appropriate. The quadratic stencil used for the Euler-Lagrange interpolation is advantageous as it reduces artificial damping of internal waves that occurs with an 8-node linear stencil; thus improving the ability of the method to resolve the free motions of a stratified basin. While this improvement is necessary for modeling stratified lakes, it is probably of lesser importance in estuarine modeling where forced motions dominate the flow physics. The extra computational requirements of the quadratic stencil are somewhat ameliorated by the ability to compute source terms for flows in the range 1 < CFL < 2 without repositioning the stencil. Numerical Techniques in CWR-ELCOM, March 2000 22 4.7 Scalar transport Scalar transport is (arguably) the most critical piece of the hydrodynamic numerical algorithm for strongly stratified flows. If the scalar transport is not sufficiently accurate, we cannot obtain correct evolution of the density field and internal wave motions. In ELCOM, a conservative third-order scalar transport method is applied. Use of conservative methods has been found critical to stratified lake and estuary models as the distribution of density feeds back into the momentum equation through the baroclinic term. Non-conservation results in loss of momentum in baroclinic forcing such that internal waves are dissipated too rapidly and strong gradients that drive underflows may simply cease to exist. Using non-conservative methods, loss of mass and deterioration of sharp gradients at salinity and temperature fronts degrades the skill of the model. A three-stage numerical algorithm for transport of a scalar concentration C can be defined as: ~ C C C n+1 = M (C n ) + S (c) ~ ~ C Uj = C - t xj C = C + x x (4.23) (4.24) + O (t) 2 (4.25) As with the momentum mixing and source terms, the mixing operator in Equation 4.23 represents vertical mixing by the Reynolds stress term in Equation 2.5 and is discussed in detail in a subsequent section on the vertical mixing model. The S (c) in Equation 4.23 represents scalar sources (e.g. heat transfer across the free surface into the wind-mixed layer). Equation 4.24 is advection of the scalar field by the resolved flow field, and Equation 4.25 is horizontal diffusion by turbulent motions. For clarity in the above and the following, advection, i.e. Equation 4.24, is defined over time step t. However, when MAX (CFLa ) > 1, the sub-time step t is used, where m t = t and Equation 4.24 is iterated m times. In the following derivations, substitution of t for t and n + m t for n + 1 makes the equations reflect the sub-time step iteration process. In differential notation, a conservative form of the transport equation is: t C d + Ax C U dAx + Ay C V dAy + Az C W dAz = S (c) (4.26) where is a control volume, Ax , Ay , Az are surface areas of the control volume faces. For clarity in the discrete form, it will be useful to omit the (i, j, k) subscript notation for any centered variable so that Ci,j,k is written simply as C, and Ci,j+ 1 ,k is written as Cj+ 1 . The discrete 2 2 advected scalar concentration field (C ) written as a flux-conservative advection over the n set of cells is: n t ~ z C = C - (4.27) x (Q) + y (Q) + z (Q) z n+1 x y z n+1 cell faces, defined on the n cell set for the i + 1/2 face as: Qi+ 1 = 2 ~ y z n+1 U n+1 C i+ 1 2 where operators of the form x are defined in Equation 4.13. Q is the scalar flux through the (4.28) Numerical Techniques in CWR-ELCOM, March 2000 23 Similar definitions apply on the j +1/2 and k+1/2 faces. There is no flux through the uppermost face of cell containing the free surface, so Qi,j,F +1/2 = 0, where k = F is the cell containing the n+1 n+1 n free surface. It follows that zi,j,F - zi,j,F = tWi,j,F +1/2 and Ci,j,F + 1 = Ci,j,F . For any 2 cell (i, j, F ) in the n set: Ci,j,F n zi,j,F n+1 zi,j,F = Ci,j,F - t n+1 1 Ci,j,F + 1 n+1 W 2 zi,j,F i,j,F + 2 t ~ z W n+1 C z n+1 (4.29) Thus, for all cells in the n set (including free surface cells): t ~ ~ x U n+1 C C = C - x - t ~ y V n+1 C y - (4.30) Since the scalar concentrations are updated at the cell centers, it is necessary to define a method of interpolation for cell face values such as x (U C). The ULTIMATE flux-limiting filter applied with third-order QUICKEST interpolation (Leonard, 1991) performs particularly well in maintaining monotonic scalar fields while limiting numerical diffusion. Implementation in 2D and demonstration of its effectiveness estuarine flows has been documented by Lin and Falconer (1997). The conservative ULTIMATE QUICKEST method is limited to CFLa < 1 in all coordinate directions. The present semi-implicit scheme remains stable at higher CFL (providing CFLb does not exceed unity), so the ULTIMATE QUICKEST algorithm must be computed successively over sub-time steps such that the maximum CFLa is less than one in each sub-time step. In practice, coarse resolution models of strongly stratified basins will have the defining time step limitation based on the baroclinic mode in the momentum solution and CFLa > 1 may never occur. The recent evaluation of transport schemes by Gross et al. (1999) did not directly examine the ULTIMATE limiter applied to a QUICKEST approach; however, their results did show that effective results for a 2D model of South San Francisco Bay could be obtained with either the QUICKEST approach or Roe's superbee limiter applied to a Lax-Wendroff method. As Lax-Wendroff is a second-order discretization, it is likely that the ULTIMATE limiter applied to third-order QUICKEST would have been at least as effective as the schemes they tested. The present use of an explicit discretization approach constrasts with the vertical-implicit approach taken in the recent adaptation of TRIM by Gross et al. (1998). The advantage of the latter is that it is stable for scalar transport solutions at CFL > 1, whereas the present method requires sub-time step iterations of the scalar transport routines when high CFL's are encountered. Thus the approach of Gross et al. (1998) is preferable for a fine vertical grid resolution while the present method should prove more computationally efficient at coarse grid resolution. 4.7.1 Scalar horizontal diffusion The horizontal diffusion terms in Equation 4.25 are discretized to obtain the time `n + 1' scalar field over the solution space `n ': n+1 ~ ~ Ci,j,k = Ci,j,k + Dx Ci,j,k ~ + Dy Ci,j,k (4.31) where Dx and Dy are finite-difference operators for the second derivative. As with the velocity solution, any time n + 1 locations that are not in the time n solution space are updated using Numerical Techniques in CWR-ELCOM, March 2000 24 the concentration of neighbor cells. 4.8 4.8.1 Mixing and vertical diffusion Introduction Parameterization of the cumulative effect of small-scale mixing events as a function of largescale (i.e. resolvable) processes remains a challenge for modeling of stratified flows. In most multi-dimensional numerical models, mixing across stable density gradients is parameterized by an eddy diffusivity term in the vertical transport equations. The mixing by sub-grid scale and sub-time scale fluctuations (i.e. turbulence) is characterized by the resolved scale strain rates using relations such as: c (4.32) x1 Turbulence closure modeling using this approach is the process of specifying the j diffusivities, and is one of the main sources of error in numerical modelling of stratified flows. However, the -cu = 1 eddy diffusivity approach is not the only possible approach to modelling the effect of mixing events. In ELCOM, one option for vertical diffusion we take a different approach that is derived from the mixing energy budgets used in 1D lake modelling. The methodology for 1D modeling is presented in Imberger and Patterson (1990), Spigel et al. (1986) and Imberger and Patterson (1981). Mixing in a stratified fluid can be characterized as two types of events: (1) convective mixing that decreases the potential energy of the fluid and releases turbulent kinetic energy (TKE) through mixing of unstable density gradients, and (2) stable mixing that uses turbulent kinetic energy to mix stable density gradients, thereby increasing the potential energy of the fluid. Mixing is computed by comparing the available mixing energy (EA ) from convective overturns, shear production and wind stirring to the potential energy increase (ER ) required to mix a lower layer up into a well-mixed region. Where EA > ER , complete mixing between the lower layer and the upper well-mixed region occurs and EA is decremented by ER . If mixing occurs, the sweep downward through the water column continues with the shear between the mixed region and the subsequent lower layer being added to EA prior to comparison with the new ER . When the downward sweep reaches the end of a well-mixed region, the remaining EA is reduced by dissipation (E ) and is stored in the lowest layer of the well-mixed region for later transport as scalar; thus mixing energy can accumulate with time. When the bottom of the wind-mixed layer is found, i.e. EA (k) < ER (k - 1), the downward sweep through the water column continues with EA computed only from the local velocity shear combined with the transported mixing energy remaining from previous time steps. The mixing algorithm proceeds in the following sequence: 1. 2. mix scalars and momentum in regions with unstable density gradients compute the TKE released by unstable mixing Numerical Techniques in CWR-ELCOM, March 2000 25 3. 4. compute additional TKE available for mixing due to wind stirring mix scalars and momentum from the free surface to depth of mixed layer (i.e. continue mixing als long as Ea > Er ) 5. compute additional TKE available due to shear production for any stable density gradient mix scalars and momentum where (Ea > Er ) decrease available TKE through dissipation (E ) advect TKE using transport equation 6. 7. 8. In the discrete implementation, vertical mixing is modelled as a process by which cells (i, j) on the discrete layer k with a volume-averaged concentration Ci,j,k are either mixed with the volumes on the k - 1 layer below or remain unmixed. This approach to vertical mixing substitutes episodic mixing "events" for the differential equation used in eddy-diffusivity models to represent vertical gradient of c u3 . Where no mixing occurs, our model neglects molecular diffusion in the vertical direction. This can be readily added either in an implicit or explicit manner, but will generally have no effect on the resulting computations. Numerical diffusion in the solution of the advective equation and horizontal diffusion across angled isotherms will dominate the vertical molecular diffusion. In discretizations of momentum and transport equations, vertical turbulent diffusion terms in Equations 2.1 and 2.5 are replaced by a mixing operator M ( ) in Equations 4.9 and 4.23, using an approach derived from the mixing energy budgets developed for 1D lake modeling (Imberger and Patterson, 1981; Spigel et al., 1986; Imberger and Patterson, 1990). In ELCOM, the mixing process is modeled on a layer-by-layer basis through each (i, j) water column as describe below. Further detail on the mixing model can be found in Hodges et al. (2000). 4.8.2 Mixing energy budgets To determine the mixing energy budget, the vertical mixing model requires computation of the energies available for mixing (EA ), required for mixing (ER ), and dissipated (E ). Derivation of the mixing energies is based upon the 3D turbulent kinetic energy transport equation: E E + Uj t xj = -Rij Ui - xj xj Euj + uj p o - - g u3 0 (4.33) + 2 2 Rij 2E + xj xj xi xj where E ui ui /2 is the turbulent kinetic energy, Rij ui uj is the Reynolds stress tensor and is the dissipation, defined as: 2 ui uj + (4.34) 2 xj xi Numerical Techniques in CWR-ELCOM, March 2000 26 If we neglect the viscous transport terms and horizontal gradients (appropriate simplifications for modeling at coarse grid resolutions), we obtain: E Ui E + Uj = -Ri3 - t xj x3 x3 Eu3 u3 p + 2 0 - - g u3 0 (4.35) which adds advective transport terms to the presentation of 1D TKE in Spigel et al. (1986). A well-mixed region of a water column (e.g. the surface wind-mixed layer, the benthic boundary layer, or a shear layer undergoing mixing) may defined in a general continuous form as a z b for some x, y, or in discrete form as ka k kb for some i, j. In a 3D system with interior shear mixing and benthic boundary mixing, there may be multiple well-mixed regions in a water column and the distribution of the regions may vary horizontally. Thus ka = ka (i, j, k) and kb = kb (i, j, k) and all k cells in a single well-mixed region of an (i, j) water column have identical values of ka and kb . For unmixed cells we have ka = kb = k. The vertical gradients in Equation 4.35 are negligible in the interior of a well-mixed region, so the right-hand side can be written as the sum of the changes through dissipation and buoyant production within the region and the production through shear, turbulent pressure work and turbulent transport at the boundaries of the mixing region: E E + Uj = - Qb + Qa - t xj wp 0 b a dz - g 0 b w dz a (4.36) Following Spigel et al. (1986), we neglect the turbulent "leakage" at the lower boundary Ew + 0 a (4.37) Then Qa and Qb are the fluxes of turbulent kinetic energy through the lower and upper boundaries of a well-mixed region: Qb Qa = uwU + vwV + Ew + wp 0 (4.38) b = {uwU + vwV }a (4.39) where U and V are the velocity differences in the vertical direction between two unmixed regions for the x and y components of velocity respectively. Thus, at the lower boundary of a mixing region (ka ), U is defined as the difference between the velocity on the ka and ka - 1 layer: U{k} = U{ka } - U{ka -1} A similar expression defines V in the y direction. The last term in Equation 4.36 is a buoyancy term that can act either as a source or a sink of TKE. Where there is an unstable density gradient, this term is a source of TKE produced in convective overturns. Where there is a stable density gradient, this is term represents a sink of TKE due to the energy required to mix heavier water upwards (i.e. increasing potential energy). This term is modeled as the buoyancy scale w based on the change in density from time t (i.e. ) to time t + t (i.e. ) due to mixing: ~ (w ) 3 (4.40) = g b a w dz = - g t ~ b a z dz - a+b 2 b dz a (4.41) Numerical Techniques in CWR-ELCOM, March 2000 27 where the density after complete mixing is given as: = ~ 1 b-a b dz a (4.42) It will be convenient to define wu as the unstable buoyancy term (TKE source due to unstable convection) such that wu = w wu = 0 for w < 0 for w 0 In a similar fashion, ws as the stable buoyancy term (TKE sink due to mixing) ws = w ws = 0 for w > 0 for w 0 This approach is taken as wu is always created by the presence of unstable density gradients, while ws is only a sink of TKE if there is sufficient mixing energy available. Equation 4.36 may be reformulated in terms of energy available for mixing, energy required for mixing, and dissipation: D Dt The sources of TKE are EA = - t uwU + vwV + Ew + wp 0 3 + {uwU + vwV }a - wu b E dz = a ER E EA - - t t t (4.43) (4.44) b The production due to wind stirring at the free surface (-Qb ) is parameterized as: 3 CN u3 = - 2 uwU + vwV + Ew + wp 0 (4.45) b where CN is an empirically determined coefficient and u is the wind shear velocity scale. The shear production between two unmixed layers ka and ka - 1 is parameterized with the coefficient CS zka -1 U 2 + V 2 = uw U + vw V 2 t so the energy available for mixing is obtained from: EA = 3 CN 3 CS u t + 2 2 3 U 2 + V 2 zka -1 - wu t CS : (4.46) (4.47) The energy required for mixing from layer k into a layer k - 1 with stable stratification is: ER = (ws ) t Finally, the dissipation of TKE is E = t b 3 (4.48) dz a (4.49) This is modeled with the coefficient CE and the available mixing energy (Spigel et al., 1986) such that CE 3 E 2 t E = (4.50) 2 A Numerical Techniques in CWR-ELCOM, March 2000 28 Using the above representations of EA , ER and E , the right-hand-side of the mixing energy evolution Equation 4.43, is computed. After sweeping vertically through the domain to determine the mixing, the advection of available mixing energy is computed using scalar transport Equation 4.26. 4.8.3 Discrete mixing operator In summary, the mixing algorithm can be thought of as a process for determining the lower (ka ) and the upper (kb ) grid cells that bound each well-mixed region in a water column. The mixing operator M ( ) modifies both the velocity and scalar fields and can be written as operating on some field as M n i,j,k = 1 n+1 i,j,k i,j,k kb (i,j,k) ( z)i,j,m m=ka (i,j,k) n : k k k (4.51) where k is the layer containing the free surface, k is the layer containing the lake bottom, and i,j,k is the total thickness of the well-mixed region that includes cell (i, j, k), computed from kb i,j,k = m=ka n zi,j,m : k k k (4.52) The new density, n+1 , is determined by a discrete form of Equation 4.42: n+1 = i,j,k 1 i,j,k kb n n zr i,j,r r=ka : k k k (4.53) Note that if ka = kb then i,j,k = zi,j,k and n+1 = n so that M () = n . The wind mixed layer is a special case of a well-mixed region whose depth (h) can be defined as: kb (i,j,k ) hi,j = m=ka (i,j,k ) zi,j,m (4.54) The mixing model requires three coefficients. These are set to empirical values available in the literature that capture the efficiencies of the above processes. An extensive discussion is found in Spigel et al. (1986), where we have taken the values: Cn = 1.33, Ce = 1.15 and Cs = 0.2. 4.8.4 Advantages and limitations of the mixing model The mixing model outlined herein is a first attempt at a 3D mixing model for lakes and therefore has some limitations in its derivation and implementation which will be subject to later improvement. The primary advantage of the mixing model is its ability to capture the correct depth of the wind-mixed layer at coarse vertical grid resolutions, and thereby obtain the first-order dynamic forcing of the thermocline4 . A secondary advantage of the present mixing model is that 4 It is not uncommon for a lake model to be limited to five or ten grid cells in the wind-mixed layer and three or four cells in the metalimnion, despite our a priori knowledge that strong gradients and curvature in profiles Numerical Techniques in CWR-ELCOM, March 2000 29 it eliminates the need for the solution of the vertical diffusion equation, thus forgoing multiple tridiagonal inversions in each water column for momentum and scalar transport required by implicit vertical diffusion methods (e.g. Casulli and Cheng, 1992). Algorithm limitations of the present mixing model include the use of a simple downward sweep through the water column to determine mixing layers. This leads to a bias in the direction of mixing that may not be appropriate for shear layers or the benthic boundary layer. In particular, this artificially limits the benthic boundary layer to a thickness of two grid cells; a problem that can be addressed by inverting the mixing algorithm for an upward sweep from the benthic boundary. Additionally the algorithm does not presently include the possibility for partial mixing that may occur as a result of billows or slow entrainment rates. Limitations of the model derivation include the neglect of the entrainment time and a characterization of shear that is fundamentally grid-dependent. The former problem leads to the assumption that the wind is capable of mixing the surface layer to the mixed-layer depth in a single time step without explicit reference to the mixed-layer depth at the previous time step; thus, the mixing dynamics are not invariant to the size of the model time step. These drawbacks do not appear to significantly affect the ability of the method to capture the basin-scale internal waves (Hodges et al., 2000). The downward bias of the mixing algorithm in shear regions and the griddependence of the shear term can be considered small errors compared to the cumulative effect of numerical diffusion that artificially thickens the metalimnion. Furthermore, the code appears to excessively damp Kelvin wave propagation around the edges of a basin so the additional frictional effects of an upward-sweeping mixing model in the benthic boundary would merely exacerbate an existing problem. The lack of temporal dynamics in the entrainment means that the code effectively sees a deeper wind-mixed layer than might be otherwise predicted during the onset of the wind. However, given the uncertainties in the spatial variation of the wind field and the numerical errors associated with coarse-grid solutions, it is difficult to argue that this is of first-order importance in modeling summer stratifications. It appears that numerical diffusion and damping are critical problems to address in this model, after which further refinements to the mixing model can be considered. 4.9 Wind momentum model The momentum input of the wind is typically modeled (e.g. Casulli and Cheng, 1992) using a stress boundary condition at the free surface: u x = u2 z= (4.55) where is an eddy viscosity and u2 is the wind stress. In ELCOM, this form of wind stress can be applied with eddy-viscosity formulation for the vertical diffusion by adding tu2 to the source of density and velocity will occur in these regions. Given the comparative scales of the physics and the grid, it is questionable as to whether one can adequately discretize partial differential equations for turbulent diffusion and transport. Thus, we should not be surprised if k - or eddy-viscosity closure schemes that prove competent in idealized test cases (with high resolution) perform poorly when used in real systems where the system size requires coarse grid resolution. Numerical Techniques in CWR-ELCOM, March 2000 30 term Gi, j, M , where k = M is the grid cell containing the free surface. As demonstrated by Glorioso and Davies (1995), the formulation chosen for computing eddy viscosity has a dramatic influence on the development of pressure-driven upwind flows in a wind-forced homogenous system. Using an eddy-viscosity/diffusivity approach in a stratified system, the depth, downwind velocity and velocity shear in the wind-mixed layer are functions of the values used for eddy viscosity and diffusivity above the thermocline. The resulting prediction of mixed-layer depth using a coarse vertical grid resolution is generally unsatisfactory. Even at fine resolutions, the capability of k - and algebraic stress models may be suspect based on the 1D results of Martin (1985). The vertical mixing model introduced in Hodges et al. (2000) eliminates the vertical diffusion term from the discrete equations, thereby eliminating the vertical eddy-diffusivity. As the purpose of the eddy-viscosity term is to model the introduction of momentum into the windmixed layer of depth h, we can substitute a model for predicting h (i.e. the mixed-layer model previously described) combined with a model for the distribution of momentum over the depth h (Imberger and Patterson, 1990): dU = dt u2 = i,j,k i,j,k i,j k Si,j,k h : -h < m=1 zm < (4.56) where is the free surface height in water column (i, j). Equation 4.56 is applied separately in the x and y directions. The mixing algorithm discussed above provides the depth of the wind-mixed layer in each water column so that the velocity in each cell can be appropriately incremented to account for the wind. Since the heat input at the surface plays an important role in the establishment of the mixed-layer depth, we use the thermodynamic equations to compute the temperature change for non-penetrative effects applied directly to the first layer of grid cells at the free surface. Short wave effects change the temperature at depth in accordance with the rate of light extinction. After the surface waters have been heated (or cooled) by the surface thermodynamics routines, the mixing routine is invoked to determine the ability of unstable density profiles and mixing energetics to mix down into the water column. Using this technique, we are ensured that the heat flux computed by the thermodynamics algorithm will be exactly introduced into the wind-mixed layer, regardless of the grid resolution. 4.10 Sidewall and bottom boundary conditions The momentum introduced into a lake by wind forcing is dissipated in the boundary layers of the lake and turbulent processes in the interior. In 2D depth-averaged models (e.g. Casulli, 1990) the Chezy-Manning bottom stress formulation has been applied to account for the dissipation in the entire water column. This is particularly useful in providing a method of calibrating depth-averaged coastal/estuarine models to reproduce a given set of tidal data. In a 3D model with stratification, detailed field data for calibrating a bottom friction coefficient is generally not available. Furthermore, without the transfer of basin-scale internal-wave energy to subgrid- Numerical Techniques in CWR-ELCOM, March 2000 31 scales waves, it is questionable as to whether any boundary condition model in the literature can capture actual boundary dynamics and predict the correct location and timing of dissipation and vertical fluxes. Sidewall boundary conditions (i.e. vertical solid boundaries) are often modeled as free-slip to effect simpler implementation in a numerical method (e.g. Casulli and Cheng, 1992). However, the lack of sidewall drag prevents vertical vorticity from being produced at boundaries. As lakes typically have regions of steep slopes and longshore currents (due to Kelvin waves), the neglect of sidewall drag may not be a suitable simplificition in modeling basin-scale motions. There remains much work to be done in developing appropriate bottom and sidewall boundary conditions in coarse grid models. ELCOM provides three basic forms of boundary conditions: (1) no-slip, (2) free-slip, and (3) specified stress5 . 5 In ELCOM version 1, the specified stress conditions are not available on sidewall boundaries Chapter 5 Future work The critical numerical modeling issues that remain unsolved are (1) the damping of the Kelvin wave, (2) numerical diffusion of the metalimnion, and (3) subgrid-scale modeling of internal wave effects. The first issue may be related to the use of no-slip lateral boundaries in the present work. Davey et al. (1983) and Hsieh et al. (1983) showed decay of model Kelvin waves to be sensitive to lateral viscosities when no-slip boundaries are used. The second issue arises due to the use of a fixed grid to discretize the regime. As the internal wave field evolves, isopycnals pass through the fixed grid and undergo numerical diffusion, resulting in vertical mixing (Griffies et al., 1999). This problem can be ameliorated with a moving isopycnal coordinate system, which introduces another set of non-trivial numerical problems (e.g. changes to grid topology due to isopycnal compression, gravity currents or unstable stratifications). To some extent, numerical diffusion of the density field can be minimized by careful choice of advective schemes and grid resolution; however, it is fair to say that numerical diffusion will, for the near future, remain the single most important issue in conducting long-term (seasonal or multi-year) prognostic or diagnostic models of stratified lakes at coarse resolution on a fixed grid. The third issue, subgrid-scale modeling of internal wave effects, is a fundamental flaw in all numerical models that do not resolve the full range of internal waves from basin-scale to the small-scale dissipative motions at the buoyancy frequency. The closure schemes adopted in the literature, whether simple eddy-viscosity closure, k - , algebraic stress, LES, or mixed-layer models, do not take into account the transfer of energy from resolved to subgrid-scale internal-wave motions. As basinscale waves degenerate into nonlinear solitary waves (Horn et al., 1998), energy is propagated through the domain and may be dissipated by wave-breaking on the boundaries (Michallet and Ivey, 1999). In a numerical model that does not resolve solitary waves, energy extracted from the basin-scale waves is dissipated locally by a turbulence model (or numerical dissipation) and generally results in local mixing (Horn et al., 1999), thus overestimating interior mixing at the expense of boundary mixing. This is particularly troublesome in the use of coupled biogeochemical/hydrodynamic models as the motivation for using 3D (vice 1D) techniques is the desire to capture the heterogeneity of the ecology. If the mixing energetics at the boundary are not captured, then the benthic boundary resuspension processes that play a key role in nutrient dynamics will be in error. 32 Numerical Techniques in CWR-ELCOM, March 2000 33 If the total energy extracted from the basin-scale waves (and dissipated locally) is a reasonable approximation of the net boundary and interior dissipation in the real system, errors in turbulence modeling should not significantly impact modeling of basin-scale motions. This is one reason the present work obtains adequate results with a relatively crude representation of boundary effects. Indeed, since basin-scale internal waves are the primary energy store for mixing, it can be inferred that correct modeling of basin-scale energetics allows field data to be applied to calibrated modeling of basin-wide mixing (a task not attempted in the present work); however, the placement and timing of mixing events cannot be calibrated. This is particularly troublesome in the use of coupled biogeochemical/hydrodynamic models as the motivation for using 3D (vice 1D) techniques is the desire to capture the heterogeneity of the ecology. If the mixing energetics at the boundary are modeled poorly, then the benthic boundary resuspension processes that play a key role in nutrient dynamics will be in error. Thus, calibration of mixing models is not a practical method for modeling 3D mass transport in a lake. Appendix Effect of hydrostatic approximation Neglecting the nonlinear and viscous terms, the z momentum equation for a non-hydrostatic flow can be written as 1 pd W = - t 0 z b a (5.1) where pd is the dynamic pressure. Integrating over a layer from a < z < b we can write pd W dz = - t 0 (5.2) Modelling the vertical velocity as piecewise linear, over a small region (z) we can say pd - 0 2 W t + k+ 1 2 W t z k- 1 2 (5.3) In the horizontal momentum equations, the neglected non-hydrostatic pressure term may be approximated as 1 pd W z W + (5.4) - 0 x 2 x t k+ 1 t k- 1 2 2 A measure of the "goodness" of the hydrostatic approximation is the relative size of the non-hydrostatic pressure gradient to the change in momentum. Let be defined as the ratio: 1 pd 0 x U t (5.5) In discrete form, the non-hydrostatic pressure gradient term in the x direction can be approximated as: 1 pd 0 x - z 2tx + W n - W n-1 W n - W n-1 i+1,j,k+ 1 2 i+ 1 ,j,k 2 - W n - W n-1 - W n - W n-1 i,j,k+ 1 2 i+1,j,k- 1 2 i,j,k- 1 2 (5.6) 34 Numerical Techniques in CWR-ELCOM, March 2000 35 Thus, 1 may be approximated by 1 2x {U n -z - U n-1 }i+ 1 ,j,k 2 W n - W n-1 + W n - W n-1 i+1,j,k+ 1 2 - W n - W n-1 - W n - W n-1 i,j,k+ 1 2 i+1,j,k- 1 2 i,j,k- 1 2 (5.7) In general, we can write that z x O (W ) O (U ) (5.8) Even in an extreme case where O (W ) O (U ) and the flow should be nonhydrostatic, equation 5.8 gives z/x. This implies that we must have aspect ratios such that z/x > O 10-2 for the nonhydrostatic term to make a significant contribution to the resolved flow field. Note that this does not imply that the nonhydrostatic terms are unimportant, rather it implies that we cannot resolve the nonhydrostatic term on a discrete grid of small aspect ratio. That is, the relative contribution of the nonhydrostatic pressure to the evolution of the horizontal flow field is a function of the aspect ratio of the discrete grid. Bibliography Al-Khafaji, A. W. and Tooley, J. R. (1986). Numerical Methods in Engineering Practice. Holt Rinehart and Winston. Amorocho, J. and DeVries, J. J. (1980). A new evaluation of the wind stress coefficient over water surfaces. Journal of Geophysical Research, 85:433442. Casulli, V. (1990). Semi-implicit finite difference methods for the two-dimensional shallow water equations. Journal of Computational Physics, 86:5674. Casulli, V. (1997). Numerical simulation of three-dimensional free surface flow in isopycnal co-ordinates. International Journal for Numerical Methods in Fluids, 25:645658. Casulli, V. and Cattani, E. (1994). Stability, accuracy and efficiency of a semi-implicit method for three-dimensional shallow water flow. Computers Math. Applic., 27:99112. Casulli, V. and Cheng, R. T. (1992). Semi-implicit finite difference methods for three-dimensional shallow water flow. International Journal for Numerical Methods in Fluids, 15:629648. Davey, M. K., Hsieh, W. W., and Wajsowicz, R. C. (1983). The free Kelvin wave with lateral and vertical viscosity. Journal of Physical Oceanography, 13:21822191. Glorioso, P. D. and Davies, A. M. (1995). The influence of eddy viscosity formulation, bottom topography and wind wave effects upon the circulation of a shallow bay. Journal of Physical Oceanography, 25:12431264. Griffies, S. M., Pacanowki, R. C., and Hallberg, R. W. (1999). Spurious diapycnal mixing associated with advection in a z-coordinate ocean model. Monthly Weather Review (in press). Gross, E. S., Casulli, V., Bonaventura, L., and Koseff, J. R. (1998). A semi-implicit method for vertical transport in multidimensional models. International Journal for Numerical Methods in Fluids, 28:157 186. Gross, E. S., Koseff, J. R., and Monismith, S. G. (1999). Evaluation of advective schemes for estuarine salinity simulations. Journal of Hydraulic Engineering, 125:3246. Hodges, B. R. (2000). Heat budget and thermodynamics at a free surface. Technical report, Centre for Water Research, The University of Western Australia. ED 1300 BH. Hodges, B. R., Imberger, J., Saggio, A., and Winters, K. B. (2000). Modeling basin-scale internal waves in a stratified lake. Limno. Oceangr (in press). Hodges, B. R. and Street, R. L. (1999). On simulation of turbulent nonlinear free-surface flows. Journal of Computational Physics, 151:425457. Horn, D. A., Imberger, J., and Ivey, G. N. (1998). The degeneration of basin-scale internal waves in lakes. Journal of Fluid Mechanics (submitted). Horn, D. A., Imberger, J., and Ivey, G. N. (1999). Internal solitary waves in lakes - a closure problem for hydrostatic models. In Proceedings of `Aha Halikoa Hawaiian Winter Workshop, January 19-22, 1999. University of Hawaii, Manoa. 36 Numerical Techniques in CWR-ELCOM, March 2000 37 Hsieh, W. W., Davey, M. K., and Wajsowicz, R. C. (1983). The free Kelvin wave in finite-difference numerical models. Journal of Physical Oceanography, 13:13831397. Imberger, J. and Patterson, J. C. (1981). A dynamic reservoir simulation model - DYRESM 5. In Fischer, H., editor, Transport Models for Inland and Coastal Waters, pages 310361. Academic Press. Imberger, J. and Patterson, J. C. (1990). Physical limnology. Advances in Applied Mechanics, 27:303 475. Iserles, A. (1996). A First Course in the Numerical Analysis of Differential Equations. Cambridge. Jacquet, J. (1983). Simulation of the thermal regime of rivers. In Orlob, G. T., editor, Mathematical Modeling of Water Quality: Streams, Lakes and Reservoirs, pages 150176. Wiley-Interscience. Kowalik, Z. and Murty, T. S. (1993). Numerical Modeling of Ocean Dynamics. World Scientific. Leonard, B. P. (1991). The ULTIMATE conservative difference scheme applied to unsteady onedimensional advection. Computer Methods in Applied Mechanics and Engineering, 88:1774. Lin, B. and Falconer, R. A. (1997). Tidal flow and transport modeling using ULTIMATE QUICKEST. Journal of Hydraulic Engineering, 123:303314. Martin, P. J. (1985). Simulation of the mixed layer at OWS November and Papa with several models. Journal of Geophysical Research, 90:903916. Michallet, H. and Ivey, G. N. (1999). Experiments on mixing due to internal solitary waves breaking on uniform slopes. Journal of Geophysical Research, 104:1347613477. Smolarkiewicz, P. K. and Pudykiewicz, J. A. (1992). A class of semi-Lagrangian approximations for fluids. Journal of the Atmospheric Sciences, 49:20822096. Spigel, R. H., Imberger, J., and Rayner, K. N. (1986). Modeling the diurnal mixed layer. Limnology and Oceanography, 31:533556. Staniforth, A. and C^t, J. (1991). Semi-Lagrangian integration schemes for atmospheric models a oe review. Monthly Weather Review, 119:22062223.

Find millions of documents on Course Hero - Study Guides, Lecture Notes, Reference Materials, Practice Exams and more. Course Hero has millions of course specific materials providing students with the best way to expand their education.

Below is a small sample set of documents:

University of Texas - SITE - 2006
Hydroinformatics 2000 Conference Iowa Institute of Hydraulic Research, 23-27 July 20001 of 14preprintModeling the Hydrodynamics of Stratified LakesB.R. HodgesThe University of Western Australia, AustraliaJ. ImbergerThe University of Wester
University of Texas - SITE - 2006
13th Australasian Fluid Mechanics Conference Monash University, Melbourne, Australia 13-18 December 1998Editors: M.C. Thompson and K. Hourigan Volume II, pp. 717 - 720WAVE-INDUCED ENSTROPHY AND DISSIPATION IN A SHEARED TURBULENT CURRENTBen R. HO
University of Texas - GRADHYDRO - 2007
CE 394K.2 Hydrology Homework Problem Set #2 Due Tues March 6 Problems in &quot;Applied Hydrology&quot; 3.5.1 Evaporation from a lake 3.6.3 Evapotranspiration computation 4.1.3 Soil water flux3.5.1.The computations are summarized In Table 3.5.1. For example,
University of Texas - GRADHYDRO - 2007
Hydrology of a Dynamic Earth A decadal research plan for hydrologic science DRAFT. Version 7.0. April 2, 2007 Released for community review and comment. Consortium of Universities for the Advancement of Hydrologic S
University of Texas - PHL - 358
University of Texas - COLA - 303
PhotographySeptember 12, 2008What do photographs do? Photographs exist. Theyre tangible; they have dimensionality. They demand our attention, our reaction, our interpretation. Its a human process, largely unconscious. They also persist. They dene
University of Texas - SESSION - 5
Getting Started in Programming Languages Implementation ResearchCristina Cifuentes Sun LabsOverview About Cristina Projects Lessons Learned Overall Advice 2007 Sun Microsystems, Inc.About Cristina Birth to PhD 2007 Sun Microsystems, Inc
University of Texas - SESSION - 1
Academic WritingWilliam Cook Univ. of Texas Austin1Encode a complex web of ideas2.as a linear stream of text. how?34about meWilliam CookHigh school drop-out PhD Brown 1989 HP Labs: Foundations of OOPLearn writing the hard wayIndustr
University of Texas - SESSION - 2
Hot Topics and Future Directions in Programming LanguagesGuy L. Steele Jr.Sun Microsystems Laboratories May 9, 2007Hot Topics and Future Directions in Programming LanguagesCopyright 2007 Sun Microsystems, Inc. (&quot;Sun&quot;). All rights are reserved
University of Texas - SESSION - 3
Getting Started in Programming Language Design ResearchGuy L. Steele Jr.Sun Microsystems Laboratories May 9, 2007Getting Started in Programming Language Design ResearchCopyright 2007 Sun Microsystems, Inc. (&quot;Sun&quot;). All rights are reserved by S
University of Texas - SESSION - 6
Career Paths: How to Get Started in Academia or IndustryCristina Cifuentes Sun LabsOverview Background Getting started in industrial research Case study research at Sun Labs Lessons learned 2007, Sun Microsystems, Inc.Background 2007, S
University of Texas - CS - 372
Disk Scheduling RevisitedMargo Seltzer, Peter Chen, John Ousterhout Computer Science Division Department of Electrical Engineering and Computer Science University of California Berkeley, CA 94720ABSTRACT Since the invention of the movable head dis
University of Texas - CS - 372
Operating Systems: Basic Concepts and History1What is an Operating System?A program and an interfaceAn abstract virtual machine A set of abstractions that simplify application designFiles instead of &quot;bytes on a disk&quot;For any OS area (CPU sch
University of Texas - CS - 372
Concurrent Programing: Why you should care, deeply1Student Questions 1. it is said that user-level threads are implemented by a library at the user-level. we have POSIX for starting user threads in C+. How do I start a kernel thread? 2. we all kn
University of Texas - CS - 372
Virtual Memory and Address Translation1ReviewProgram addresses are virtual addresses.Relative offset of program regions can not change during program execution. E.g., heap can not move further from code. Virtual addresses = physical address inc
University of Texas - CS - 372
Network and Distributed File Systems1From Local to Network File SystemSo far, we have assumed that files are stored on local disk . How can we generalize the design to access files stored on a remote server? Need to invoke file creation and mana
University of Texas - CS - 372
ProcessesLast Time Each process maintains its own state, that includes its text and data, procedure call stack, etc. The OS also stores process state for each process. This state is called the Process Control Block (PCB), and it includes the PC,
University of Texas - CS - 372
Mutual Exclusion: Primitives and Implementation Considerations1Too Much Milk: Lessons Software solution works, but it is unsatisfactorySolution is complicated; proving correctness is tricky even for the simple example While thread is waiting, it
University of Texas - CS - 372
Distributed Coordination1Topics Event Ordering Mutual Exclusion Atomicity of Transactions Two Phase Commit (2PC) DeadlocksAvoidance/Prevention DetectionThe King has died. Long live the King!2Event Ordering Coordination of requests (especia
University of Texas - CS - 372
Concurrent Programming Issues &amp; Readers/Writers1Summary of Our Discussions Developing and debugging concurrent programs is hard Safety: isolation and atomicity Scheduling: busy-waiting and blocking Synchronization constructsNon-deterministic int
University of Texas - CS - 372
Protection and SecurityHow to be a paranoid or just think like one12The ProblemTypes of misuseAccidental IntentionalProtection and security objectiveProtect against/prevent misuseThree key components:Authentication: Verify user ident
University of Texas - CS - 372
Memory Management Basics1Basic Memory Management ConceptsAddress spacessupported by the hardwarePhysical address space The address spaceStarting at address 0, going to address MAXsysMAXsysview of its own memoryLogical/virtual address
University of Texas - CS - 372
OS Structure, Processes &amp; Process Management1What is a Process?A process is a program during execution.Program = static file (image) Process = executing program = program + execution state.A process is the basic unit of execution in an operat
University of Texas - CS - 372
Introduction to I/O and Disk Management1Secondary Storage ManagementDisks just like memory, only differentWhy have disks?Memory is small. Disks are large.Memory is volatile. Disks are forever (?!)File storage.Short term storage for mem
University of Texas - CS - 372
Introduction to Operating SystemLast TimeAn operating system is the interface between the user and the architecture.User Applications Virtual Machine Interface Operating System Physical Machine Interface HardwareOS as juggler: providing the il
University of Texas - CS - 372
Semaphores and Monitors: High-level Synchronization Constructs1Synchronization Constructs SynchronizationCoordinating execution of multiple threads that share data structuresPast few lectures:Locks: provide mutual exclusion Condition variab
University of Texas - CS - 372
File System AbstractionLast TimeApplications Daemons Programmer Interface Device Indepedent Interface Device Interface Open() Close() Link() Sectors Seek()ServersShellRead() Write() Rename() Tracks WriteBlock()ReadBlock() Hardware DiskH
University of Texas - CS - 372
File Systems: Consistency Issues1File Systems: Consistency IssuesFile systems maintain many data structuresFree list/bit vector Directories File headers and inode structures Data blocksAll data structures are cached for better performanceWor
University of Texas - CS - 372
Experience with Processes and Monitors in Mesa1Butler W. Lampson Xerox Palo Alto Research Center David D. Redell Xerox Business SystemsAbstractThe use of monitors for describing concurrency has been much discussed in the literature. When monitors
University of Texas - CS - 372
Parallel ComputingBasics of Parallel Computers Shared Memory SMP / NUMA Architectures Message Passing Clusters1Why Parallel Computing No matter how effective ILP/Moore's Law, more is betterMost systems run multiple applications simultaneously
University of Texas - CS - 372
Condition Synchronization1Synchronization Now that you have seen locks, is that all there is? No, but what is the &quot;right&quot; way to build a parallel program.People are still trying to figure that out.Compromises:between making it easy to modify