[August 24, 2011] A Scale-Invariant Formulation of the Anticipated Potential Vorticity Method [Monthly Weather Review] Tweet (Monthly Weather Review Via Acquire Media NewsEdge) ABSTRACT The long-term success of climate models that operate on multiresolution grids depends on access to subgrid parameterizations that act appropriately across a wide range of spatial and temporal scales. As the first step in a series of efforts to obtain such scale-aware subgrid parameterizations, the authors focus on the anticipated potential vorticity method (APVM) on a sequence of quasi-uniform grids with varying resolutions. Through a scale analysis technique and phenomenological theories for two-dimensional turbulent flows, they derive a new formulation of the APVM that depends on a single parameter that is formally independent of the time-step size, the grid resolution, and the flow itself. Results of numerical experiments with this new formulation demonstrate that the optimal parameter of the new APVM formulation is invariant with respect to the time-step size, is insensitive to the flows, and is only weakly dependent on the grid resolution. (ProQuest: ... denotes formulae omitted.) 1. Introduction The pressing need to have detailed knowledge of local features of climate systems and the high costs of obtaining such knowledge has motivated research on regional climatemodeling. Anew global-to-regionalmultiresolution approach has been proposed for which only one grid, with variable resolutions and with smooth transitions between fine and coarse grid regions, is involved (Ringler et al. 2011). This approach gives rise to a new set of challenges. One challenge is how to produce suitable meshes with local resolutions (see, e.g., Ringler et al. 2008). Another challenge is related to the implementation of parameterizations on a variable-resolution grid. This work is part of a series of efforts to address the second issue. Developing scale-aware parameterizations for the atmosphere and ocean has been a difficult and largely unmet challenge. While the closures for clouds in the atmosphere (Arakawa and Schubert 1974) and eddies in the ocean (Gent and McWilliams 1990) have clearly been successful, neither has been generalized across spatial and/or temporal scales. The long-term success of models that operate on meshes with multiple resolutions will depend on access to closure parameterizations that act appropriately across a wide range of spatial and temporal scales with little or no ad hoc tuning. Such closures, assuming that they exist, would demonstrate an understanding of the modeled physical process(es) far beyond our current capability. Two of the most important components of the climate system, the ocean and the atmosphere, can both be considered as thin layers of fluids surrounding the earth. In other words, the large-scale motions of the ocean and atmosphere are nearly two-dimensional flows and, thus, the fundamental principles of two-dimensional turbulent flows provide guidance into the dynamics of the ocean and atmosphere. Themost striking features of such turbulent flows are the inverse cascade of the kinetic energy and the cascade of the enstrophy (Batchelor 1969; Kraichnan 1967; Lilly 1971, 1969; Maltrud and Vallis 1991). For two-dimensional turbulent flows, the enstrophy plays the role that the kinetic energy does for threedimensional turbulent flows. In adiabatic systems, such as the shallow-water equations or the primitive equations cast in isopycnal coordinates, this role is played by the potential enstrophy. The presence of the upscale (inverse) energy cascade in two-dimensional turbulence acts tomove kinetic energy away fromthe grid scale; thus requiring very little, if any, grid-scale dissipation of kinetic energy in numerical simulations. At the same time, the downscale cascade of potential enstrophy necessitates the removal of grid-scale variance in potential vorticity in order to maintain robust simulations. The fundamentally different behavior of energy and enstrophy in twodimensional turbulence has led to the development of energy-conserving, enstrophy-dissipating numerical schemes (see, e.g., Arakawa and Hsu 1990; Sadourny and Basdevant 1985). The anticipated potential vorticity method (APVM) of Sadourny and Basdevant (1985) is a subgrid eddy closure scheme that is perfectly energy conserving and, at the same time, enstrophy dissipating. In the context of a shallow-water flow, the APVM can be described by ... (1) with u and gh denoting the velocity and geopotential fields, respectively; ... is the potential vorticity; D is a correction to the potential vorticity q; and k is the unit vector in the vertical direction. The modified system remains energy conserving, as the term k × u is perpendicular to the horizontal velocity component. As discussed in detail in section 2 and the appendix, the termD is chosen so that the modified system dissipates the potential enstrophy (i.e., the variance of the potential vorticity). In its simplest form, we have ... (2) It was proposed in Sadourny and Basdevant (1985) that g1 should be a time-scale-selective parameter and thus should take the following form: ... (3) with dt denoting the time step and s is an independent parameter. After substituting (3) into (2), we see that, with dt × u being the displacement vector and ?q being the change rate of the potential vorticity over space, D represents the variation of the potential vorticity along the fluid path. As a result, the quantity (q - D) represents the value of the potential vorticity at an upstream point, hence, the name anticipated potential vorticity method. As we have mentioned, we intend to study parameterizations on variable-resolution grids. Thus, we are pursuing approaches for the development of scaleaware parameterizations. As the first step, here we focus on the APVM (1)-(2) on a sequence of quasi-uniform grids with varying resolution. We ask the following questions: what is the optimal form-value of the coefficient g and how does the optimal parameter change in response to changes in the grid resolution and time step? Clearly, the primary challenge in the development of any scale-aware parameterization will be the identification of parameter(s) that are largely insensitive to the spatial and temporal resolution of the numerical model. So whereas the form of ? in (3) is appealing because of its clear physical interpretation, it is clearly sensitive to the choice of model time step. There is no analysis available to support the invariance of the parameter s when using the APVM; in fact, our experience shows that the optimal value of s is indeed influenced by the time-step size and the grid resolution of each particular simulation. In our opinion, this puts a severe constraint on the applicability of the form in (3) for ? because, for each particular simulation, without extensive fine-tuning and comparisons, it is not clear what the optimal value of s one should use. If we can cast theAPVMin terms of a scale-insensitive parameter, then we can use basic parameter optimization techniques to find the appropriate value of this parameter through comparisons to high-resolution reference solutions. Parameter optimization is really only practical for scale-insensitive parameters. To overcome the difficulty associated with the form in (3) for g, we endeavor to find a form of the APVM parameterization that is a function of a single, largely scale-invariant parameter, which is denoted as a in our analysis below. Our main tools are a scale analysis and the phenomenological theories of two-dimensional turbulent flows. The nonlinear advective term plays an important role in the phenomenological theories of turbulence (see, e.g., Frisch 1995). However, the nonlinear advective term usually does not participate in the scale analysis carried out in designing subgrid eddy closure schemes (Berselli et al. 2006). The dissipation term always has a key place in scale analysis because the dissipation rate is assumed to be representative of the eddy fluxes throughout the inertial range. A nonlinear dissipation term, such as the one withAPVM[see (17) as well as (A10) and (A11)] will lead to nonlinear convoluted terms in the contributions to the (potential) enstrophy dissipation at a given scale [see (A15) for the two-dimensional incompressible case]. We deal with the difficulty associated with these terms by analyzing the scale interactions and singling out the dominant term, which is then taken as an approximation to the (potential) enstrophy dissipation rate at that scale. Details are presented in section 2 for shallow-water flows and in the appendix for two-dimensional incompressible flows. In section 3, we present the results of numerical experiments conducted to determine the optimal value of a that appears in the new form of ?. One role of those experiments is to show that the scale analysis is valid, despite the assumptions made. In fact, it is seen from the experiments that a is invariant to the size of the time step and has a very weak dependence on the grid resolution. In section 4 we provide the concluding remarks. 2. Subgrid eddy closure based on the APVM In this section, we perform a scale analysis on the APVM applied to the shallow-water equations, and derive a new formulation that depends on a single parameter that is formally independent of the time-step size and the grid resolution and therefore is suitable for parameter optimization. The scale analysis is an extension of the analysis given in the appendix for twodimensional incompressible turbulent flows for which the phenomenological (Batchelor 1969; Kraichnan 1967; Leith 1968) and mathematical (Bardos 1972; Constantin 2007) theories are mature. This extension is possible if we assume that the shallow-water flow is predominantly two-dimensional and variations in the fluid thickness variable h are negligible compared to the mean fluid thickness. For many realistic cases, this is a reasonable assumption. The APVM method for the system of the one-layer shallow-water equations is given by ... (4) where h denotes the fluid thickness and ... is the kinetic energy. For shallow-water flows the potential vorticity is defined as the ratio between the absolute vorticity and the fluid thickness: ... (5) with the absolute vorticity h given by ... (6) and the relative vorticity ? given by ... (7) The term D is the APVM correction to the potential vorticity. Taking the curl of the second equation of (4), we obtain ... (8) Because the Coriolis force can be considered time independent, we then write (8) as ... (9) Replacing ? by ?q in (9) and using the first equation of (4) for the height h, we infer that ... (10) Multiplying (10) by q, we obtain ... (11) Multiplying the first equation of (4) by q2/2, we obtain ... (12) Combining (11) and (12), we have ... (13) LetMdenote the spatial domain. IfMis bounded, then the nonpenetration boundary condition is assumed on the boundary (i.e., u * n = 0, n being the outward normal vector on the boundary ?M). IfMis the whole sphere, then no boundary conditions are needed. Integrating (13) overM, we obtain the potential enstrophy balance equation: ... (14) It is easy to see that, if the correction term D is taken as ... (15) then the right-hand side of (14) is negative. Indeed, with this choice of D, (14) becomes ... (16) Because the flow is predominantly two dimensional and the variations in h are small compared to h itself, we assume that h on both sides of (16) is constant. The implications of this assumption are discussed in section 4. Then, that equation reduces to ... (17) Let l denote the grid scale and ?l denote the enstrophy dissipation rate at the scale '. Noticing the relations ..., we derive from the left-hand side of (17) that ... (18) where ql and ?l are the potential vorticity and streamfunction, respectively, at the scale l and the symbol "~" means ''equal within an order one constant.'' From the right-hand side of (17), we infer that ... (19) In the above we have replaced u by u0, the magnitude of the zeroth mode of the velocity field, based on the conjecture that the term ... dominates in all the contributions to the potential enstrophy dissipation at the scale l. This conjecture can be supported by arguments similar to thosemade in the appendix for two-dimensional Euler flows. Comparing (18) and (19), we obtain ... (20) Using again the relation ..., we can rewrite (20) as ... (21) We remark that (21) has the dimension of time and this coincides with the traditional choice of sdt as the form of ?; see Ringler et al. (2011). The expression (21) is not of use in practice because ?l is not readily available. We shall apply the same trick as used in the appendix to equate ?l to the overall model potential enstrophy dissipation rate ?. To do so, we first determine the relation between ql and ?l According to (5)-(7), the potential vorticity q is defined as ... It still holds true that, as in the two-dimensionalEuler case (see the appendix), variations in the Coriolis parameter f are unimportant at the grid scale '. On the other hand, we continue to use the assumption that the variation in h is small compared to h itself and replace h in the above equation by its mean value h. Hence, we have ... (22) With (22), the potential enstrophy dissipation rate ?l in (18) takes the form ... (23) Using this new form for ?l, we can infer from (20) that ... (24) Note that the overall model enstrophy dissipation rate is given by ... (25) Replacing ?l in (24) by x and using (25), we have ... (26) The full form for ? is then given by ... (27) where a is a constant independent of the scale and state of the flow. We remark that the average kinetic energy ke, the advection of the potential vorticity u * ?q, and the average height h can all be computed from the corresponding dynamical variables of the model. The grid scale l can be represented by the distance between two neighboring grid points. Thus, a is the only parameter in the expression in (27) for ? and is the value of this parameter that we wish to optimize. 3. Numerical experiments a. The spectra and the basic optimization technique Our purpose here is to evaluate the effectiveness of the scale-aware APVM closure developed above. In particular, several reasonable sounding assumptions were made in carrying out the scale analysis, so that numerical testing of the scale invariance of a is useful in verifying that result of the analysis. As with the evaluation of most closures, the measure of effectiveness is based on the ability of a low-resolution simulation with the closure to reproduce certain important aspects of a high-resolution, reference simulation. In the simulations below, we evaluate the closure based in its ability to reproduce the spectrumof the enstrophy obtained in the high-resolution simulation. It has been conjectured (see, e.g., Batchelor 1969; Kraichnan 1967) and approximately verified in the literature (see, e.g., Lilly 1969, 1971) that for two-dimensional incompressible turbulent flows, the spectrum of the enstrophy satisfies a21 power law. Equivalently, the spectrum of the kinetic energy satisfies a23 power law, in contrast to the famous25/ 3 power law for three-dimensional turbulent flows. Our reference solution is computed on a very fine grid, with the least damping that is still capable of producing a potential enstrophy spectrum that approximates the 21 power law well. The damping mechanism used to produce the reference solution is provided by conventional hyperviscosity that has been widely used in geophysical fluid dynamics. To study the spectra of the fluids on the whole sphere, the natural choice of basis functions are the spherical harmonic functions (Adams and Swarztrauber 1997). We briefly recall that the spherical harmonics are the eigenfunctions of the Laplace operator on the sphere. In spherical coordinates, spherical harmonics are conventionally represented as ... (28) where m denotes the zonal wavenumber and Pmn the nthis the associated Legendre polynomial of mth order. Here Ymn (?, f) satisfies the eigenvalue problem: ... (29) Because of (29),...is often taken as the scalar (spherical) wavenumber of the spherical harmonic Ymn . It can be shown that {Ymn} forms a complete orthogonal basis for functions on the sphere (Courant and Hilbert 1953). For this reason, the complete set of the spherical harmonics makes an ideal candidate for the basis functions used for studying the spectra of various dynamical variables on the sphere. Because we have chosen to define the effectiveness of the closure based on its ability to reproduce the potential enstrophy spectra of the reference solution, our definition of ''error'' represents the difference, on a log scale, between the spectra over a subrange of wavenumbers. Hence, for a spectral density function P and the reference spectral density function Pr, we define the difference as ... (30) where the spherical wavenumber k is defined as [see (29)] ... and k0 and k1 are the starting and ending wavenumbers for the subrange, respectively.Werefer to d as the spectral distance. b. The optimal parameters The numerical experiments are conducted with the one-layer shallow-water equations on the whole sphere. The numerical model is based on the numerical scheme presented in Thuburn et al. (2009) and Ringler et al. (2011). This approach is being used to develop global atmosphere and ocean models that are a part of the Model for Prediction Across Scale (MPAS) project. The MPAS modeling approach is particularly attractive for this study since it solves the vector-invariant form of the momentum equation from which the APVM is derived in (4), conserves potential vorticity in analogy to the continuous system in (9), and guarantees that the APVM correction term does not create or destroy kinetic energy. While the MPAS system can be used with a wide class of meshes, for our experiments we use spherical centroidal Voronoi tessellations (SCVTs) because of their global uniformity and isotropy (Du et al. 1999). The remapping from SCVTs to regular Gaussian grids on the sphere (as input to the spherical harmonics package; Adams and Swarztrauber 1997) is done with the Spherical Coordinate Remapping and Integration Package (SCRIP; Jones 1999). The first test case is the standard shallow-water test case 5 (Williamson et al. 1992) that involves a zonal flow impinging on a mountain. The initial zonal velocity u of the flow and the surface height h are given by ... respectively. The mountain, which is part of the lower boundary, has the form ... The fluid thickness h is then given by ... In the above, ? represents the latitude as usual, and the other physical parameters are set following Williamson et al. (1992), namely, ..., with l being the longitude, and .... The midlatitude Rossby radius in these simulations is approximately 2000 km. The coarsest simulation conducted has a nominal resolution of approximately 450 km. Thus, in all simulations the Rossby radius is resolved. We remark that the topography is not a source or sink of potential vorticity, but instead acts as an inhomogeneity that leads to the breakdown of the nonlinearly balanced geostrophic zonal jet. There is no other external forces added, and therefore this is a slowly decaying turbulent flow. The reference solution is computed on an SCVT mesh with 655 362 cells (dx [asymptotically =] 30 km), with the traditional ?4 dissipation of 109 m4 s21. The shallow-water test case 5 evolves into turbulence in about 20-25 days. A snapshot of the potential vorticity field on day 50 is shown in Fig. 1. The axis of the globe is slightly tilted in order to better display the structure of the potential vorticity field on the Northern Hemisphere. We plot the spectrum of the potential enstrophy of the reference solution on day 150 in Fig. 2. An inertial range of width approximately 1 decade appears between wavenumbers 20 and 120, which approximately verifies the 21 power law for potential enstrophy spectra. Before we present comparison and optimization results, we remark on the choice of the starting and ending wavenumbers of the subrange used in the comparisons. The subrange for comparison should be part of the inertial ranges of the reference and approximate solutions. The reason for this is that the high noise level in the scales larger than the inertial range will render the comparison results largely stochastic, and hence noninformative, whereas the spectra at scales smaller than the inertial range heavily depend on the damping mechanism used in each particular simulation, and hence cannot give reliable results. Because the approximate solutions are computed using the APVM on a coarser grid to measure the effectiveness of the closure, these simulation display a spectral range that is narrower than that of the reference solution. Their inertial range will also be narrower than that of the reference solution. This fact should also be taken into consideration when choosing the starting and ending wavenumbers of the subrange used for comparisons. Approximate solutions are computed on a 10 242-cell (approximate resolution of 240 km) SCVT grid, with the APVMparameter a taking values from 0 to 0.01, with an increment of 0.0001. The APVM in (15) is the only closure used in these low-resolution simulations (i.e., the ?4 dissipation is turned off). Hence, there are 101 simulations. The spectrum of the potential enstrophy for each of these simulations is compared to that of the reference solution and a distance between the spectra is calculated using (30). The approximate solutions on the 10 242-cell grid have an inertial range between wavenumbers 20 and 80. Hence, we take k0 = 20 and k1 = 80 in (30), per the discussion in the preceding paragraph. Then, we plot the spectral distance against the APVM parameter a in Fig. 3. It is seen that the minimum distance is obtained at a = 0.0013. Just to show what happens if a is too small or too large, we plot in Fig. 4 the potential enstrophy spectra of the reference solution and of the simulations with a = 0, 0.0013, 0.0080. With a = 0, the potential enstrophy is transferred down scales by the nonlinearity and then, because of the lack of dissipation, it piles up at high wavenumbers. On the other hand, with a = 0.008, the overdamping is obvious as it pulls down the tail of the spectrum significantly below that of the reference solution leading to an inertial range with a slope significantly steeper than the -1 slope. The APVM, as it was formulated in Sadourny and Basdevant (1985) and used in the literature, is not suitable for parameter optimization because the way it is formulated, the parameter will depend on the time-step size as well as the grid resolution. The advantage of our formulation is that the parameter is invariant with regard to the time-step size, up to time truncation error. It is also formally independent of the grid resolution, though we suspect that the parameter may have a weak dependence on the grid resolution. In what follows, we demonstrate the time-step invariance aspect of our formulation and we explore whether and how the parameter depends on the spatial grid resolution. To show the invariance of the parameter with respect to the time-step size dt, we redo the comparison study, but with dt varying over f691.2 s, 345.6 s, 172.8 s, 86.4 s, 43.2 s, 21.6 sg. For each of these step sizes, a curve is plotted in Fig. 5 showing how the distance in the spectra varies with respect to the APVM parameter a. These curves for different time-step sizes agree very well with each other. In particular, they all point to the same minimizing parameter a = 0.0013. This is not surprising, as our modified subgrid eddy closure model consists of partial differential equations that are independent of the time-step size, and hence the solutions of the model are independent of the time-step size as well. To explore whether and how the parameter a depends on the grid resolution, we fix the time-step size at dt = 172.8 s and perform the optimization on a on 2 additional grids, one coarser with 2562 cells, and the other finer with 40 962 cells. For each of these grids (nCells = 2562, 10 242, and 40 962) we plot, in Fig. 6, a curve depicting how the spectral distance changes with respect to the APVM parameter a. The curve for the coarser grid (nCells = 2562) attains its minimum at a = 0.0010 whereas the curve for nCells = 10 242 attains its minimum at a = 0.0013, as we have already seen. The curve for the finer grid (nCells = 40 962) attains its minimum at a = 0.0018. Through this experiment, we see that the minimizing parameter a tends to increase as the grid is refined. The dependence of the minimizing parameter a is weak because, in this experiment, the grid resolution changes by a factor of 2, but the change in the minimizing a is much slower. As a matter of fact, the minimizing a is approximately proportional to nCells0.212, with the exponent 0.212 being much smaller than the exponent 3 over l in (27). There are a few other noteworthy features about Fig. 6. We see that for the same parameter a, finer grids usually produce better results (smaller spectral distance when compared to the reference solution). We also note that as the grid gets finer, the spectral distance curve becomes less sensitive to the increment in the parameter a. In fact, even though the curve for nCells = 40 962 attains its minimum at a = 0.0018, its value at a = 0.0013 is very close to the minimum value. To corroborate on this point, we compare, in Fig. 7, the spectra of the results on the finer grid (nCells = 40 962) with a = 0.0010, 0.0013, 0.0018, and 0.0050. The spectrum of the reference solution is also plotted for comparison. We see that the spectra for a = 0.0010, 0.0013, and 0.0018 stick together and are close to the spectrum of the reference solution in the inertial range. On the other hand, the spectrum for a = 0.0050 is pulled down below the reference solution due to overdamping. These results indicate that, although the minimizing value of the parameter a is weakly dependent on the grid resolution, there is some robustness present: using optimal values for the parameter determined using one grid resolution can seemingly be safely used for other grid resolutions. Universality is a key concept in Kolmogorov's K41 theory of three-dimensional turbulence (Kolmogorov 1941a,b). It is also one that has been the subject of intense debate (see, e.g., Frisch 1995; Landau and Lifshitz 1987). It is beyond the scope of the current article to discuss the implications that this concept bears on twodimensional turbulence modeling. Instead, we examine the sensitivity of the optimal APVM parameter with respect to the flow itself by performing another set of experiments with arbitrarily imposed initial data. In the results that we present, the initial velocity is derived from the streamfunction: ... where ? denotes the latitude and [straight phi] the longitude. This seemingly simple function is highly oscillatory in the spherical harmonics space, leading to the presence of a wide range of scales. The height h is initially taken as a constant 5000 m. There is no external forcing applied, and we let the flow evolve freely. Therefore, this is also a slowly decaying turbulent flow. The highresolution solution is again computed on the 655 362- cell SCVT grid, with a ?4 diffusion of 109 m4 s21. On day 230, we observe an inertial range spanning approximately from wavenumber 30 to 300 (see Fig. 8). The APVM in the approximate solutions in (4) are computed on the 10 242-cell SCVT grid. The errors in the inertial range of potential enstrophy spectra of these approximate solutions are plotted against the APVM parameter a in Fig. 9. The trend in the spectral distance with respect to the APVM parameter a resembles that depicted in Fig. 3. The minimizing parameter found for this test case, a = 0.0020, is close to the optimal value found for the standard shallow-water test case 5. 4. Conclusions The original form of the coefficient ? for the anticipated potential vorticitymethod, as suggested in Sadourny and Basdevant (1985), is not suitable for parameter optimization because it involves a parameter that is not invariant with respect to the time-step size and the spatial grid resolution. Using a scale analysis technique and the phenomenological theories of two-dimensional turbulence, we propose a newformfor theAPVMcoefficient ? such that the new parameter involved is formally invariant with respect to the time-step size and spatial grid resolution. Numerical experiments are conducted on the whole sphere with different grids, each of which is a quasiuniform spherical centroidal Voronoi tessellation of the sphere. Two test cases have been used. One is the standard shallow-water test case 5 involving a mountain topography, and the other is a flow evolving from arbitrarily imposed initial data. Our basic optimization technique is to compare the potential enstrophy spectra of the APVM solutions to that of a reference solution that is calculated on a very fine grid with ?4 hyperviscosity. The numerical results demonstrate the time-step invariance aspect of our formulation for ?. Over a sequence of grids having the number of cells varying from 2562 to 40 962, the optimal APVM parameter is found to be within the range 0.001 to 0.002. Because the measure function is relatively flat near the optimal value of a, this factor of 2 change in the optimal parameter is not significant. One issue that we have not touched upon in this article is the performance of the APVM compared to the traditional hyperviscosity (iterated D) method. This issue was discussed in Sadourny and Basdevant (1985) for the APVM in its original formulation; it was shown that, for the same spatial grid resolution, the APVM produces more realistic results than the traditional hyperviscosity method. A comparison between the APVM in the new scale-invariant formulation and the traditional hyperviscosity method is then in order. With the standard shallow-water test case 5, the reference solution is computed on a 655 362-cell (approximate resolution of 30 km) grid with a ?4 dissipation of 1.0 × 109 m4 s-1, which is the result of extensive fine tuning. Basic scaling arguments say that the ?4 parameter goes like l4. Thus, the optimal ?4 parameter on a 10 242-cell grid should be 4.0 × 1012 m4 s-1. For comparison, we plot, in Fig. 10, the spectra between the spherical wavenumbers 10 and 100, of the high-resolution reference solution, the 4.0 × 1012 m4 s-1 ?4 solution, and the APVM solution with the optimal parameter a = 0.0013. It is seen that the solution of theAPVMand the solution of the ?4 method match very well at low wavenumbers; at high wavenumbers the APVM solution has a more extensive inertial range that matches well with that of the reference solution. The spectrum curve of the ?4 solution also appears slightly steeper than that of the reference solution in the inertial range. This test shows that the APVM in the new formulation developed in this article is at least as good as the traditional ?4 method at producing solutions that conform to the phenomenological theories of two-dimensional turbulence. The full three-dimensional primitive system will require some type of ?4 closure to dissipate the downscale cascade of kinetic energy. Even in this situation, we expect that the APVM is a valuable component of the overall model closure by allowing the ?4 parameters to be smaller, and therefore less dissipative, than they would be otherwise. The closure presented here includes the assumption of small deviations in fluid thickness. Namely, before starting the scaling analysis in (17), we assume that the fluid thickness h that is included on both sides of (16) cancels out. For the simulations conducted above and, likely, for the broad class of shallow-water systems where this closure might be utilized, this will, in general, be a valid assumption. As we apply the closure to more realistic systems, such as the isopycnal model of the ocean circulation, this assumption will have to be carefully reevaluated. This work is the first in a series of efforts to address the issue of subgrid eddy parameterizations on variableresolution grids in global models. The experiments in this work are conducted on quasi-uniform grids and we address the question about how the resolution of the quasi-uniform grid affects the optimizing parameter for the APVM. Naturally, our next step is to address the issue on a variable-resolution grid. There, the question becomes: how the resolution of each region affects the optimal parameter a. The form of (27) for ? involves the grid resolution l and therefore the regional resolution has been accounted for in the coefficient ?. Hence, we believe the form of (27) can be the starting point for parameter optimizations on variableresolution grids. However, we also expect new challenges (e.g., the interactions between the flows from different regions). Our first step in doing parameter optimizations is to identify a parameter that is invariant with respect to other configurations of the model (e.g., the time-step size and the spatial grid resolution). We believe that this approach applies to and can be taken toward other parameterizations or subgrid eddy closure schemes. It is our intention to apply this approach to certain other important parameterizations (e.g., the Gent-McWilliams parameterization of turbulent transport on variableresolution grids; Gent and McWilliams 1990; Gent et al. 1995; Ringler and Gent 2011). Acknowledgments. The authors owe thanks to Mat Maltrud for a careful reading of a draft of this paper and for helpful comments. The authors also thank the anonymous referees, whose comments and suggestions helped to improve the manuscript. Q. Chen and M. Gunzburger were supported by theU.S.Department of Energy Grant DE-SC0002624 as part of the Climate Modeling: Simulating Climate at Regional Scale program. T. Ringler was supported by the DOE Office of Science's Climate Change Prediction Program Grant DOE 07SCPF152. 1 In the foregoing reference, ? instead of ? is used for the APVM coefficient. We have switched to ? to avoid confusion with the spherical coordinate u introduced later. REFERENCES Adams, J. C., and P. N. Swarztrauber, 1997: Spherepack 2.0: A model development facility. NCAR Tech. Note NCAR/TN- 436-STR, 59 pp. Arakawa, A., and W. Schubert, 1974: Interaction of a cumulus cloud ensemble with the large-scale environment, Part I. J. Atmos. Sci., 31, 674-701. _____, and Y.-J. G. Hsu, 1990: Energy conserving and potentialenstrophy dissipating schemes for the shallow water equations. Mon. Wea. Rev., 118, 1960-1969. Bardos, C., 1972: Existence et unicité de la solution de l'équation d'Euler en dimension deux (Existence and uniqueness of the solution of the two dimensional Euler equations). J. Math. Anal. Appl., 40, 769-790. Batchelor, G. K., 1969: Computation of the energy spectrum in homogeneous two-dimensional turbulence. Phys. Fluids, 12 (12), 233-239. Berselli, L. C., T. Iliescu, and W. J. Layton, 2006: Mathematics of Large Eddy Simulation of Turbulent Flows. Springer, 366 pp. Constantin, P., 2007: On the Euler equations of incompressible fluids. Bull. Amer. Math. Soc., 44, 603-621. Courant, R., and D. Hilbert, 1953: Methods of Mathematical Physics. Vol. I. Interscience Publishers, xv 1 561 pp. Du, Q., V. Faber, and M. Gunzburger, 1999: Centroidal Voronoi tessellations: Applications and algorithms. SIAM Rev., 41 (4), 637-676. Frisch, U., 1995: Turbulence: The Legacy of A.N. Kolmogorov. Cambridge University Press, 312 pp. Gent, P. R., and J. C. McWilliams, 1990: Isopycnal mixing in ocean circulation models. J. Phys. Oceanogr., 20, 150-155. _____, J. Willebrand, T. J. McDougall, and J. C. McWilliams, 1995: Parameterizing eddy-induced tracer transports in ocean circulation models. J. Phys. Oceanogr., 25, 463-474. Jones, P. W., 1999: First- and second-order conservative remapping schemes for grids in spherical coordinates. Mon. Wea. Rev., 127, 2204-2210. Kolmogorov, A. N., 1941a: Dissipation of energy in locally isotropic turbulence. Proc. U.S.S.R. Acad. Sci. (Atmos. Ocean. Phys.), 32, 16-18. _____, 1941b: The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. Proc. U.S.S.R. Acad. Sci. (Atmos. Oceanic Phys.), 30, 299-303. Kraichnan, R. H., 1967: Inertial ranges in two-dimensional turbulence. Phys. Fluids, 10 (7), 1417-1423. Landau, L. D., and E. M. Lifshitz, 1987: Course of Theoretical Physics. Vol. 6, Fluid Mechanics, 2nd ed. Pergamon Press, xiv 1 539 pp. Leith, C. E., 1968: Diffusion approximation for two-dimensional turbulence. Phys. Fluids, 11 (3), 671-672. Lilly, D. K., 1969: Numerical simulation of two-dimensional turbulence. Phys. Fluids, 12 (12), 240-249. _____, 1971: Numerical simulation of developing and decaying two-dimensional turbulence. J. Fluid Mech., 45 (2), 395-415. Maltrud, M. E., and G. K. Vallis, 1991: Energy spectra and coherent structures in forced two-dimensional and beta-plane turbulence. J. Fluid Mech., 228, 321-342. McWilliams, J. C., 1989: Statistical properties of decaying geostrophic turbulence. J. Fluid Mech., 198, 199-230. Pedlosky, J., 1987: Geophysical Fluid Dynamics. 2nd ed. Springer, 728 pp. Rhines, P. B., 1975: Waves and turbulence on a beta-plane. J. Fluid Mech., 69, 417-443. Ringler, T., and P. Gent, 2011: An eddy closure for potential vorticity. Ocean Modell., in press. _____, L. Ju, and M. Gunzburger, 2008: A multiresolution method for climate system modeling: Application of spherical centroidal Voronoi tessellations. Ocean Dyn., 58, 475-498. _____, J. Thuburn, J. B. Klemp, and W. C. Skamarock, 2010: A unified approach to energy conservation and potential vorticity dynamics for arbitrarily-structured c-grids. J. Comput. Phys., 229 (9), 3065-3090. _____,D. Jacobsen,M.Gunzburger, L. Ju,M.Duda, and W. Skamarock, 2011: Exploring a multiresolution modeling approach within the shallow-water equations. Mon. Wea. Rev., in press. Sadourny, R., and C. Basdevant, 1985: Parameterization of subgridscale barotropic and baroclinic eddies in quasi-geostrophic models: Anticipated potential vorticity method. J. Atmos. Sci., 42, 1353-1363. Smagorinsky, J., 1963: General circulation experiments with the primitive equations. Mon. Wea. Rev., 91, 99-164. Thuburn, J., T. Ringler, W. Skamarock, and J. Klemp, 2009: Numerical representation of geostrophic modes on arbitrarily structured c-grids. J. Comput. Phys., 228, 8321-8335. Williamson, D. L., J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber, 1992: A standard test set for numerical approximations to the shallow water equations in spherical geometry. J. Comput. Phys., 102, 211-224. QINGSHAN CHEN AND MAX GUNZBURGER Department of Scientific Computing, The Florida State University, Tallahassee, Florida TODD RINGLER Los Alamos National Laboratory, Los Alamos, New Mexico (Manuscript received 1 October 2010, in final form 18 January 2011) Corresponding author address: Qingshan Chen, Department of Scientific Computing, The Florida State University, Tallahassee, FL 32306. E-mail: [email protected] APPENDIX Scale Analysis for the 2D Incompressible Flows In this appendix,we give details about the scale analysis for two-dimensional inviscid incompressible flows, which leads to a scale-aware formulation for the anticipated potential vorticity method (APVM). Such an analysis serves as a base for the extension to shallow-water flows given in section 2. We also hope that it will provide guidelines for the analysis of other geophysical flows that are predominantly two-dimensional, such as the primitive equations cast in isopycnal coordinates. Two-dimensional inviscid incompressible flows obey the Euler equations: ... (A1) The term fk × u represents the Coriolis force, as the equations are cast in a rotating reference frame. Before we present the APVM, we need to rewrite (A1) in the vector-invariant form: ... (A2) where ... (A3) denotes the absolute vorticity. The APVM adds a correction D to ? in (A2) so that it becomes ... (A4) Wefirst remark that this correction term has no effect on the energy balance of the system because ... In other words, the kinetic energy is still conserved for the modified system in (A4). Taking the curl of the velocity equation and using the incompressibility condition shown in (A4), we obtain ... (A5) We shall determine a form of D so that the term on the right-hand side of (A5) will effectively dissipate enstrophy. We denote the spatial domain by M. If M is bounded, then the no penetration boundary condition is assumed (i.e., u * n = 0, n denoting the outward normal vector on the boundary ?M). IfMis the whole sphere, then no boundary conditions are needed. We multiply (A5) by ?, integrate by parts over the spatial domainM, and again making use of the incompressibility condition, we have ...(A6) The right-hand side of (A6) is negative if ... (A7) for some positive definite operator A, where A being positive definite means, for each scalar function J: ... (A8) A natural choice for the operator A is the iterated Laplacian operator: ...(A9) with ... being a positive parameter yet to be determined. Wecan verify thatAis positive definite according to (A8) via integration by parts, assuming suitable boundary conditions on ? and u to annihilate the boundary terms (again, on the whole sphere, no boundary conditions are needed). In this appendix we analyze the general form (A9) ofA. The analysis will be extended to a special form (m = 0) of A for the shallow-water model in section 2. The numerical experiments reported on in section 3 are conducted with this special form as well. We note that with A taking the form (A9), (A5), and (A6) become ... (A10) and ... (A11) respectively. In what follows we employ the scale-analysis technique to determine a suitable form for ... Before doing so, it will be helpful to recall that, given the streamfunction ?, the velocity field u and the absolute vorticity ? can be obtained via ... (A12) We denote by l the grid scale, and by a subscript ? the mode of certain variables at the scale ? (e.g., ?l for the mode of the streamfunction ? at the scale l and xl for the mode of the enstrophy dissipation rate at the scale l). At the grid scale l, the flow is assumed to be in the turbulent regime; isotropy is assumed for motions at this scale. Therefore, it is justified to replace the overall enstrophy dissipation rate by the per-unit-volume enstrophy dissipation rate (Berselli et al. 2006). The absolute vorticity h and the streamfunction are related by (A12). We first note that the variations in the Coriolis parameter f are unimportant at the grid scale l. A rigorous justification of this assertion would requires the use of Fourier expansion or spherical harmonics. Formally speaking, the Coriolis parameter f is only present in the first few modes (large scales), and is absent from higher modes (smaller scales). Thus at l, we have ... (A13) with ?l denoting the magnitude of the mode of the absolute vorticity at the scale l. Noticing that ..., we derive from the left-hand side of (A11) that ... (A14) The enstrophy dissipation rate xl can also be calculated using the right-hand side of (A11). However, there is a pitfall as the right-hand side involves the integral of the square of a nonlinear term. We have briefly discussed the difficulty associated with this term in the introduction. This fourth-order nonlinear term comes from the nonlinear term on the right-hand side of (A10) multiplied by h. Generally speaking, in the presence of nonlinearity, the scale separation principle no longer holds. However, it can be argued that ... is the dominant term in the enstrophy dissipation rate at the scale ', with u0 denoting the magnitude of the zeroth mode of the velocity field in its representation in terms of, for example, the spherical harmonics. Indeed, we go back to (A10), multiply the right-hand side by ?l, and we see that the convoluted terms that contribute to the enstrophy dissipation at the scale l has the following form: ... (A15) with l1, l2, '3 =l satisfying the following relation: ... (A16) In the above, for i = 1, 2, 3, the ul and ?li are the component of the velocity and the absolute vorticity, respectively, at the scale li. By (A12), these two variables are related to ?li the component of the streamfunction ? at the scale li in the following way: ... Substituting uli and ?li in (A15) by the corresponding expressions listed above, we obtain ... (A17) Phenomenological theories (Batchelor 1969; Kraichnan 1967; Leith 1968) of two-dimensional turbulence implies that the enstrophy cascade rate ? is a constant in the inertial range, which includes the grid scale ?. Therefore, letting L denote the largest scale in the inertial range (or, the origin of the enstrophy tunnel), and assuming ..., we infer from (A14) that ... (A18) Now that each ?li can be written in terms of ?L, li, and L, we replace the ?li s in (A17) by their corresponding expressions so that, also using (A16), we obtain ... (A19) It is clear that this term is large when both l1 and l2 are large (or equivalently when l3 is small) and it is small when both l1 and l2 are small (or equivalently when l3 is large). When l1 and/or l2 fall outside the inertial range, the relation (A18) does not hold and the preceding analysis is no longer valid. However, it is reasonable to conjecture that when l1 and l2 are the largest scales and l3 coincides with l, the term in (A19) dominates in all terms that contribute to the enstrophy dissipation at the scale l. Such a conjecture is supported by the postulation of inverse energy cascade in two-dimensional turbulent flows (see e.g., McWilliams 1989; Pedlosky 1987; Rhines 1975), as u0, the magnitude of the largest scale of the velocity field, tends to grow and dominate as the flow evolves. Well-behaved numerical results presented in section 3, albeit for the shallow-water flows, also support this conjecture. Thus, we replace ul1 and ul2 in (A15) by u0, and ?;3 by ?l, and take the resulting expression as an estimate for ?l. Therefore, we have ... (A20) In practice, u0 can be calculated by taking mean of the kinetic energy ke over the spatial domain. Thus (A20) is equivalent to ... (A21) Comparing (A21) with (A14), we find that ... (A22) Using (A14) again, we express ... in terms of ?l and the scale l: ... (A23) The expression in (A23) is not ready for practical use because, formany realistic turbulent flows, the enstrophy dissipation rate ?l is not known a priori. To overcome this difficulty, we apply a trick that is often used in constructing subgrid eddy closure schemes (e.g., the classic Smargorinsky closure; Smagorinsky 1963). The trick is to equate ?l with the overall model enstrophy dissipation rate ?. The rationale for doing this is that if the closure method is to faithfully model the flow, the dissipation mechanism at the scale l must be able to replicate the overall model enstrophy dissipation rate. To apply the trick to our case, we first derive from (A11) that ... (A24) Replacing ?l in (A23) by x and using (A24), we find ... The complete form of ... is therefore ... (A25) where ... is a constant independent of the grid resolution and of the state of the flow. (c) 2011 American Meteorological Society [ Back To TMCnet.com's Homepage ]