Numerical diffusion in stratified lake models

in proceedings of the Fifth International Symposium on Stratified Flows, Vancouver, Canada, July 7-11, 2000. Pg 343-348 Numerical Diffusion in Stratified Lake Models Bernard Laval, Ben R. Hodges, Jörg Imberger Centre for Water Research, University of Western Australia Nedlands, Perth,Western Australia, 6907, Australia. blaval@cwr.uwa.edu.au, hodges@cwr.uwa.edu.au, jimberger@cwr.uwa.edu.au 1. Abstract Three-dimensional numerical models that solve the shallow water equations on a coarse fixed-grid, are limited in their ability to predict internal wave evolution, by the cumulative effects of numerical diffusion. Numerical diffusion of the metalimnion alters the internal response of the basin, and precludes simulations on seasonal to annual time scales. The discrete solution of the equations of motion leads to unintentional smoothing of advected gradients, which irreversibly increases the thickness of the metalimnion. The hydrostatic approximation exacerbates numerical diffusion by steepening internal basin-scale waves due to nonlinear effects until limited by numerical diffusion and dissipation. This steepening enhances horizontal gradients at the wave front which increases numerical diffusion locally, leaving a diffused metalimnion in the wake of the wave front. A method for correcting the numerical diffusion of mass based on conservation of background potential energy is proposed. 2. Introduction As wind blows over the surface of a lake, it sets up a basin-scale tilt in the pycnocline which evolves into long internal waves. These waves steepen due to non-linear effects until the rate of steepening is balanced by dispersion, leading to solitary waves (Horn et al., 1999). Upon reaching a shoaling boundary, these solitary waves will break losing a considerable amount of their energy to local irreversible mixing (Michallet and Ivey, 1999). The mixing induced by the shoaling and breaking of these waves is an important component of the net vertical transport in a stratified lake (Imberger, 1998). The discrete solution of the Navier-Stokes equations on a fixed grid leads to unintentional smoothing of gradients known as artificial or numerical diffusion. Numerical diffusion affects any fluid properties that are advected, including momentum and any active (ie. salt, heat, sediment load) or passive scalars (ie. nutrients, dissolved oxygen, plankton, color, dye) as shown in Figure 1. Considerable research has gone into advection methods which minimize numerical diffusion. However, these still require a relatively fine grid mesh to be effective (ie.: Leonard, 1991). Numerical diffusion of momentum (numerical dissipation), leads to enhanced decay of kinetic energy, while numerical diffusion of active scalars leads to enhanced diapycnal flux of mass. The former is of lesser importance when forcing time scales are shorter than that associated with the decay of motion. The latter has a cumulative effect on the density stratification which changes the timescales of the internal response of the water body. As the relative phase of wind forcing and internal waves changes, so too does the coupling between the wind forcing and pycnocline setup, altering the internal response to a given wind (Antenucci et al., 2000). Computer power limits the total number of grid points used in a simulation. This limitation combined with the large horizontal to vertical aspect ratios of most natural water bodies, leads to large aspect ratios of the numerical grid mesh if a reasonable vertical resolution is sought. A coarse grid mesh limits the effectiveness of including the unsteady term in the vertical momentum equation, and so it is reasonably neglected, using the shallow water, or hydrostatic approximation. This Figure 1: Schematic diagram of numerical diffusion. Consider two regions having property concentration c1 and c2 , separated by a tilted interface that is sharply defined along the cell edges in a). Horizontal advection in b) introduces c2 water in c1 cells requiring spatial averaging due to the discrete grid, which results in water with property concentration c3 in c). approximation combined with numerical dissipation and numerical diffusion substantially alters the dynamics of the simulation results by removing the dispersive effects from long wave dynamics. This allows non-linear steepening of basin-scale waves to evolve until limited by numerical diffusion of mass and momentum. This paper describes three-dimensional, coarse-grid numerical simulations of a wind-forced saltstratified basin to demonstrate the effects of numerical diffusion. A method is then proposed to correct for the numerical diffusion of mass and forestall the cumulative effects that lead to numerical thickening of the metalimnion over seasonal time scales. 3. Simulations The simulations described in this paper were performed using the numerical code ELCOM (Hodges, 2000; Hodges et al., 2000). ELCOM solves the three-dimensional, Reynolds-averaged, Navier-Stokes equations using the hydrostatic and Boussinesq approximations. The Reynolds stress terms are approximated in the horizontal as eddy diffusivities, and a mixed-layer model based on mixing energy budgets is used in the vertical dimension. To minimize the effects of numerical diffusion, ELCOM uses the ULTIMATE QUICKEST method for the solution of the scalar transport equations (Leonard, 1991). All simulations were performed in a rectangular domain with seiche motions parallel to the sidewalls. The domain was one cell wide, 24 km long and 25 m deep, with a 1 m layer thickness. The domain had a constant temperature of 20◦ C, while surface waters had a salinity of 4 parts per thousand (ppt), deep waters had a salinity of 8 ppt, with a smooth gradient region between 10 an 15 m. For the purposes of demonstrating the effects of numerical diffusion, all horizontal diffusivities and mixing were set to zero. Free-slip boundary conditions were used, so that all dissipation was numerical. An initially tilted pycnocline was released at time zero and subsequent motions allowed to develop freely. This simulation was preformed at two different horizontal resolutions to demonstrate that the amount of non-linear steepening of basin-scale waves is highly dependent on grid resolution (Figure 2). The wave front at 250 m horizontal resolution is much steeper than that simulated at 1000 m horizontal resolution. Due to the grid-mesh aspect ratio, horizontal advection of steep fronts leads to large amounts of diapycnal numerical diffusion. While the wave front simulated at high resolution appears similar in form to a packet of solitary waves, the diapycnal flux occurs at the wave front rather than at the domain boundary as is the case for true solitary waves. a) b) Figure 2: Contour plot of salinity from the numerical simulation of a free, internal, basin-scale seiche at grid aspect ratios of a) 250x1m, and b) 2000x1m. Higher horizontal grid resolution leads to steeper wave fronts, which can lead to greater numerical diffusion. The free surface is at 0.25 m depth and the bottom at 25 m depth As a second set of simulations to demonstrate the long term effects of numerical diffusion a uniform wind was applied across the surface at 5 ms−1 for one half of the fundamental seiche period (Ti /2) based on the initial stratification, and subsequently no wind was applied until the wind-induced motions had died down (∼30 days), whereupon a 5 ms−1 was reapplied for Ti /2. This continued for one year of simulation time. The mixing model coefficients altered to minimize mixing while still allowing the transfer wind momentum into the surface layer to induce seiching motions. Horizontal diffusivities remained zero. Figure 3a shows the initial salinity field, and figure 3b shows the salinity field after 365 days of simulation. The pycnocline thickness has nearly doubled over the duration of the simulation, while the surface salinity has increased by 0.6 ppt and the bottom salinity has decreased by 0.2 ppt. The asymmetry in the salinity changes at the surface and bottom are likely due to a small amount of wind stirring by the mixed layer model. However, the dramatic increase in the pycnocline thickness is due to numerical diffusion. 4. Correction Technique The background potential energy, Eb , is the potential energy of a system brought to rest adiabatically (e.g. Winters et al., 1995). Eb is dependent on the adiabatically relaxed density distribution which can only change through irreversible mixing, or boundary fluxes of heat or mass. Pure advection cannot occur across diapycnals, hence does not alter the adiabatically relaxed density distribution and therefore does not change Eb . ELCOM calculates mixing and advection separately, so the change in Eb during scalar transport computations is due to numerical diffusion. The change in Eb during scalar transport is then a a) b) c) Figure 3: Contour plots of salinity from numerical simulations of a periodically wind forced internal basin-scale seiche with a grid size of 1000x1m. a) shows the initial salinity field. b) and c) show the salinity field at the end of one year for unfiltered and filtered simulations (see §4). Wind forcing was 5 ms−1 for Ti /2 once a month and 0 ms−1 otherwise. The free surface is at 0.25 m depth and the bottom at 25 m depth. measure of the domain integrated numerical diffusion in a time step. The goal of our present work has been to apply our knowledge of the integrated numerical diffusion to compute a local correction to the density field that counters the diffusion over time. Our technique applies a sharpening filter to each water column such that: |Eb− − {Eb+ }| < threshold (1) is satisfied; where Eb− and Eb+ are the background potential energy before and after the computation of the advection respectively, the curly braces represent the application of the sharpening filter, and threshold is a predetermined quantity. The sharpening filter is a three point recursive filter based on the inverse Z-transform of a linear response to a heavy-side step function (Fozdar et al., 1985). The filter was originally developed and is presently used to enhance field instrument response. In order to avoid biasing, and to simulate the bi-directionality of diffusion in 1D, the filter is applied both backwards and forwards in the vertical dimension in each water column. Application of the filter to the wind-driven simulation described in the previous section allows for the effects of numerical diffusion to be “undone”. Figure 3c shows the salinity field after 365 days of simulation with the filter applied. Surface salinity has increased by ∼0.2 ppt, bottom salinity remains unaltered, while the pycnocline is about 1 m deeper and remains approximately the same thickness as it did initially. The increase in salinity of the surface layer and deepening of the pycnocline are likely due to the wind stirring induced by the mixing model. 5. Conclusion Three-dimensional numerical modeling of stratified lakes is limited by its inability to resolve the energy cascade from basin-scale internal waves to solitary waves, break at the boundaries. These limitations arise from a combination of the hydrostatic assumption and a coarse grid mesh, which results in large amounts of numerical dissipation of momentum and numerical diffusion of mass. The hydrostatic assumption leads to the non-linear steepening of basin-scale internal waves into bores which leave a trail of numerically diffused properties. Hence, mixing and dissipation are occurring in the interior of the domain rather than at the boundaries. Also, diffusion of the pycnocline reduces the period of basin-scale seiches, which changes the relative phase of internal motion and wind forcing, altering their coupling. A method to correct for the numerical diffusion of mass is presented which is effective in removing the effects of numerical diffusion over seasonal time scales. This method improves the ability of numerical models to simulate the internal wave field over similar time scales. However, the long term net transport requires . . . Acknowledgments This research was supported by the Centre for Environmental Fluid Dynamics at the University of Western Australia, and by an International Postgraduate Scholarship. 6. References Antenucci, J. P., Imberger, J., and Saggio, A. (2000). Seasonal evolution of the basin-scale internal wave field in a large stratified lake. Limnology and Oceanography. in press. Fozdar, F. M., Parker, G. J., and Imberger, J. (1985). Matching temperature and conductivity sensor response characteristics. Journal of Physical Oceanography, 15:1557–1569. Hodges, B. R. (2000). Numerical rechniques in CWR-ELCOM. Technical report, Centre for Water Research, University of Western Australia, Nedlands, Western Australia, 6907. Hodges, B. R., Imberger, J., Saggio, A., and Winters, K. (2000). Modeling basin-scale internal waves in a stratified lake. Limnology and Oceanography. in press. Horn, D., Imberger, J., and Ivey, G. (1999). Internal solitary waves in lakes — a closure problem for hydrostatic models. In Proceedings of the 11th ’Aha Haliko’a Hawaiian Winter Workshop, University of Hawaii, Manoa. Imberger, J. (1998). Flux paths in a stratified lake: A review. In Imberger, J., editor, Physical Processes in Lakes and Oceans, pages 1–17. American Geophysical Union. Leonard, B. P. (1991). The ULTIMATE conservative difference scheme applied to unsteady onedimensional advection. Computer Methods in Applied Mechanics and Engineering, 88:17–74. Michallet, H. and Ivey, G. N. (1999). Experiments on mixing due to internal solitary wavesbreaking on a uniform surface. Journal of Geophysical Research, 104(C6):13467–13478. Winters, K. B., Lombard, P. N., Riley, J. J., and D’Asaro, E. A. (1995). Available potential energy and mixing in density-stratified fluids. Journal of Fluid Mechanics, 289:115–128.