The utility of static and adaptive mesh refinement (SMR and AMR, respectively) are examined for idealized tropical cyclone (TC) simulations in a two-dimensional spectral element f-plane shallow-water model. The SMR simulations have varying sizes of the statically refined meshes (geometry based) while the AMR simulations use a potential vorticity (PV) threshold to adaptively refine the mesh to the evolving TC. Numerical simulations are conducted for four cases: (i) TC-like vortex advecting in a uniform flow, (ii) binary vortex interaction, (iii) barotropic instability of a PV ring, and (iv) barotropic instability of a thin strip of PV. For each case, a uniform grid high-resolution “truth” simulation is compared to two different SMR simulations and three different AMR simulations for accuracy and efficiency. The multiple SMR and AMR simulations have variations in the number of fully refined elements in the vicinity of the TC. For these idealized cases, it is found that the SMR and AMR simulations are able to resolve the vortex dynamical processes (e.g., barotropic instability, Rossby wave breaking, and filamentation) as well as the truth simulations, with no significant loss in accuracy in the refined region in the vortex vicinity and with significant speedups (factors of 4–15, depending on the total number of refined elements). The overall accuracy is enhanced by a greater area of fully refined mesh in both the SMR and AMR simulations.
Atmospheric motions span a multitude of different spatial and temporal scales. Examples are planetary waves at spatial scales of 106 m that evolve over days to boundary layer turbulent eddies at scales of 101 m that evolve over minutes. With current computational resources, it is not possible to simulate the entire spectrum of atmospheric flows. One of the goals in the design of next-generation numerical weather prediction (NWP) models is that they be unified, or that one nonhydrostatic dynamical core has the capability of simulating a wide spectrum of atmospheric spatial and temporal scales of motion, from microscale to global, and weather to climate. Severe and high-impact weather can often take the form of localized weather systems, such as severe thunderstorms, tornadoes, fronts, and tropical cyclones. Often these weather systems have steep gradients in various fields such as moisture and vorticity. With limited computational resources, it would be preferred to perform local mesh refinement to resolve the details of these features, while resolving the large-scale features (e.g., synoptic-scale high pressure systems) at coarser resolution.
Currently, the primary method for tackling this scale discrepancy is by utilizing multiply nested NWP models (Kurihara et al. 1979; Hodur 1997; Kurihara et al. 1998; Skamarock et al. 2005; Doyle et al. 2014). However, a number of drawbacks exist with this method. First, there are multiple lateral boundaries, often with the existence of nonphysical blending zones. Second, there is inefficiency in performing the same forecast on each nest since the nests are embedded within each other. Third, because of the extra communication required between nests, it is not expected that these setups would scale that efficiently on multiple processors. An alternative method to embedded nests is static or adaptive mesh refinement. In the former method, a mesh could be statically refined over a region of interest (e.g., a city or coastline) providing more finescale details of the flow there. In the latter method, the mesh could adaptively refine and de-refine based on some feature of interest, such as a tropical cyclone (TC). The earlier work of Berger and Oliger (1984) and Skamarock and Klemp (1993) demonstrated the utility of adaptive mesh refinement (AMR) for hyperbolic equations. More recently, the application of AMR to the shallow-water equations has been examined (Berger et al. 2009; LeVeque et al. 2011; McCorquodale et al. 2015). A review of the current state of AMR for atmospheric modeling is described in Jablonowski (2004), Behrens (2006), and Weller et al. (2010).
There are many challenges that remain for the application of variable-resolution meshes for real TC prediction in NWP models. Foremost among these challenges are implementation of scale-aware physics and determination of suitable adaptation criteria. However, the use of statically refined variable-resolution meshes in climate simulations has yielded some promising results in comparison to uniformly refined grids (Zarzycki et al. 2014). Prior to real TC prediction with AMR in NWP models, many aspects of the application of AMR and static mesh refinement (SMR) to TC simulations must first be examined in an idealized framework. As such, the purpose of the present study is to examine the utility of both SMR and AMR for ideal TC simulations in a next-generation dynamical core. The model is a planar spectral element shallow-water model, with similar numerical methods used in the Nonhydrostatic Unified Model of the Atmosphere (NUMA; Giraldo and Restelli 2008; Kelly and Giraldo 2012; Giraldo et al. 2013). We examine the utility of SMR and AMR for four flows, representing idealizations of TC dynamics in the real atmosphere. First, we examine a TC advecting in a uniform flow, representing a TC tracking in the atmosphere in the environmental flow. Second, a binary vortex interaction is examined, representing the interaction of two TCs that are close together. Third, instabilities and mixing processes are examined in the hurricane inner core (eye and eyewall). Fourth, the instability of a simplified intertropical convergence zone (ITCZ), its breakdown, and formation of TC-like vortices are examined. In each case, we compare a series of SMR and AMR simulations with variable regions of refined mesh to a “truth” simulation with a uniform refined mesh in order to obtain an understanding of efficiency and accuracy trade-offs. While it remains to be seen whether or not AMR would be feasible for real hurricane prediction in future NWP models, this work is intended to provide a first step to this problem by examining its applicability in very idealized frameworks.
The remainder of this paper is organized as follows. In section 2, the continuous model equations and numerical method are given. The test cases and model setup are given in section 3. The initial conditions and results of each experiment are presented in section 4. Diagnostics of the results are given in section 5, including a quantitative analysis of the accuracy and efficiency. The main conclusions are presented in section 6.
2. Model equations and numerical method
a. Continuous equations
The numerical model is based upon the flux form of the f-plane shallow-water equations, which can be written in compact vector form as follows:
where , is the velocity vector (where the superscript T denotes the transpose), with h being the fluid height, and g is the gravitational constant. Other quantities requiring definition include , which is an identity matrix; is the vector pointing upward (along the z coordinate that coincides with the direction along which h is measured); denotes the tensor product operator; and denotes the divergence operator. An important property of the unforced, inviscid shallow-water equations (1) and (2) is the material conservation of potential vorticity,
where is the potential vorticity normalized by the mean and time-independent fluid depth , where is the relative vorticity, and is the material derivative.
b. Numerical method and mesh refinement algorithms
The shallow-water equations (1) and (2) are discretized using the continuous Galerkin (CG) numerical method. High-order CG methods for the shallow-water equations are given in Ma (1993) and Taylor et al. (1997), and the numerical model used in this study is based upon the specific methods discussed in Giraldo (2001) and Giraldo and Restelli (2008). The CG method has been applied to numerous idealized test cases, as well as more complicated idealized tests cases of atmospheric phenomena, such as two-dimensional (Gabersek et al. 2012) and three-dimensional (Marras and Giraldo 2015) moist experiments of a squall line.
A brief overview of the CG method is given here. Given a computational domain , the domain is first decomposed into a number of elements as
where is one element. Within each element the local contributions to the integrals appearing in the weak form equations are computed, and the solution is expanded as
where is a prognostic variable, since square elements are used, N is the polynomial order, and is a local basis function. In the CG method, neighboring elements share interface points. The solution is chosen to be obtained at the Legendre–Gauss–Lobatto (LGL) nodal points. LGL points allow for the use of collocated interpolation and integration [since LGL points have integration strength of O()]. Equally spaced points are not used as the interpolation points because the Lebesgue constant will grow and the interpolation will be poor.
The mesh refinement algorithm is based upon Kopera and Giraldo (2014, 2015), and uses a forest of quad trees (i.e., each element may have four children), similar to the approach used by St-Cyr et al. (2008). In the refinement procedure, the polynomial order N is held constant, while individual elements are refined. For the experiments here, a maximum of two levels of mesh refinement is used, so that a fully refined element is 4 times the horizontal resolution of a coarse element. Additionally, the refinement algorithm includes the functionality to generate an arbitrary number of layers of refined mesh cells extending away from the feature of interest. This is hereafter referred to as the “buffer” region. The AMR criterion for this study is the normalized potential vorticity . The normalized potential vorticity has the same units of s−1 as vorticity. Refinement and coarsening of the elements is accomplished based on different thresholds in P. An attractive feature of using PV for mesh adaptation is that linear inertia–gravity waves contain zero PV, eliminating the possibility of AMR tracking fast-mode inertia–gravity waves. This is of particular interest for this study since we are interested in the slow-manifold PV dynamics. The SMR criterion is based subjectively on different area sizes around the feature of interest. More details on the refinement methods, including how nonconforming elements, which are generated with refinement, are handled in the CG method, and the grid-to-grid interpolation method used during the AMR simulations, can be found in Kopera and Giraldo (2015).
3. Cases, model setup, and error norms
Four different test cases are examined and described in detail below. These cases are as follows: (i) TC vortex moving in a uniform flow; (ii) binary vortex interaction; (iii) dynamic instability of the hurricane eyewall, eye mesovortex formation, and mixing between the eyewall and eye; and (iv) formation of TC-like vortices from the barotropic instability of a shear zone. These cases are two-dimensional idealizations of TC processes occurring in the real atmosphere. To relate the simulations to the real processes, Fig. 1 shows an example real-case scenario that each test case is designed to represent. In Fig. 1a, Hurricane Andrew is shown moving west toward Florida, being advected by the easterly flow around the subtropical ridge to its north (case-1 idealization). In Fig. 1b, the binary vortex interaction of two storms, Typhoons Melor and Parma, are shown. The differential advection induced from each cyclone advects the other creating a net cyclonic motion (the Fujiwhara effect; Fujiwhara 1921). Binary vortex interactions may also occur at a much smaller scale inside a tropical cyclone, as two convectively generated PV anomalies interact with one another (Kuo et al. 2004) leading to the formation of concentric eyewalls (case-2 idealization). In Fig. 1c, the instability and breakdown of the eyewall of Hurricane Dolly (2008) is shown, leading to an asymmetric radar reflectivity pattern there (case-3 idealization). Finally, in Fig. 1d, the instability and breakdown of the ITCZ is shown over days. The deep convection along the ITCZ is observed to undulate and finally to break down into distinct tropical cyclones (case-4 idealization).
b. Discretization and model setup
Cases 1–3 use a 600 × 600 km2 domain and case 4 uses a 8000 × 8000 km2 domain due to the larger feature of interest (ITCZ). The fluid depth is set to m for each experiment, corresponding to a gravity wave phase speed of approximately 140 m s−1. For all simulations in section 4, fifth-order polynomials () are used in each element, and a fourth-order explicit Runge–Kutta scheme is used for the temporal integration. For each initial condition listed above, six numerical simulations are performed: (i) a high-resolution truth simulation (FINE), which has a uniform grid of fully refined elements; (ii) a large statically refined mesh around the TC processes (SMR2); (iii) a smaller statically refined mesh (SMR1); (iv) adaptive mesh refinement with a buffer of six fully refined elements (AMR3); (v) adaptive mesh refinement with a buffer of three fully refined elements (AMR2); and (vi) adaptive mesh refinement with no buffer (AMR1). The FINE numerical simulation is the truth, and is expected to simulate the phenomenon with the most accuracy. The FINE simulation has a uniform grid with the same resolution as the finest parts of the SMR and AMR simulations.
The varying SMR and AMR simulations are intended to help retain the accuracy of the solution in the vicinity of the TC while also saving on computational time. All simulations are run at the time step of the FINE simulation since our goal is understanding computational aspects only with regard to the spatial variation. The time steps for each simulation are set using a Courant number of approximately 0.25 using the gravity wave phase speed above. The simulation parameters and grid description for all experiments are given in Table 1. The first part of the table lists the initial number of elements for each simulation, followed by the element sizes and effective resolutions for the three possible mesh refinement levels: small, middle, and large. The large element is 2 times larger than the middle element, and 4 times larger than the small element. Since polynomials are used in each element, an approximate effective resolution can obtained by dividing the element size by a factor of 5. However, the actual minimum grid spacing is less than this number since the LGL points are unequally spaced, and closer together near the element boundary. All simulations are run with no diffusion (inviscid); however, in order to control aliasing instability, a modal filter is used. The modal filter is the Boyd–Vandeven filter, which transforms the solution on an element to modal space, trims the highest modes using a prescribed filter function (in our case we simply remove 5% of the highest modes), and transforms the solution back to nodal space. Detailed description of this filter applied to element-based Galerkin methods can be found in Giraldo et al. (2013). For the AMR simulations, the grid is checked for adaptation based on a threshold in PV every 10 model time steps. The numerical method herein uses the same numerical representation for height and velocity (equal order at the same LGL nodal points) and so will produce spurious computational modes. However, as LeRoux (2012) points out, many methods produce some form of spurious modes, which must be controlled by some dissipation mechanism. Here, we control the modes by the use of the Boyd–Vandeven filter.
c. Error norms
To quantify aspects of these results, normalized and errors were computed between the FINE, or truth, simulation, and the SMR/AMR simulations. To compute the and error norms, the solution at the AMR/SMR meshes is interpolated to the FINE mesh using a two-dimensional L2 projection as described in Kopera and Giraldo (2014). The errors are computed at the simulation end times for each case: (i) case 1: h, (ii) case 2: h, (iii) case 3: h, and (iv) case 4: h. The normalized and errors are defined as
where q is the predicted variable, N is the total number of points, are the nodal points, the superscript “” denotes the numerical simulations, and the superscript “” denotes the reference solution (or in this case the FINE solution). The normalized errors are computed for the magnitude of the velocity vector and the geopotential . Since we are interested in how well the SMR/AMR simulations resolve the local TC processes in comparison to the FINE simulations, the errors are computed in two regions: (i) the entire domain and (ii) in the localized region in which the AMR/SMR simulations are designed to resolve. For (ii), this region was defined as km for cases 1, 2, and 3, and km for case 4. These correspond to the regions of large PV in each case (or localized region of interest).
We also show contour plots of the local error, defined as
a. Case 1: TC vortex advecting in a uniform flow
The first test case is a TC-like vortex advecting in a uniform flow, which is an idealization of a TC moving with the environmental flow in the atmosphere. The initial vortex is constructed as a Rankine vortex in polar coordinates [, where and ϕ is the azimuthal angle in radians] according to
where is the tangential velocity, s−1, and km. With these parameters, the peak tangential velocity at km is 25 m s−1. A smooth radial decay function is multiplied by the tangential winds so that at , with a cutoff radius km, and beyond . The vortex Cartesian momentum components and are next specified and then the uniform zonal flow m s−1 is added to u. This experiment is conducted on an f plane with . The initial balanced fluid depth is determined by solving the nonlinear balance equation for the balanced departure of the height field:
using the CG method in conjunction with the generalized minimal residual method (GMRES), where , , and . The right-hand side of the nonlinear balance equation is computed directly using the horizontal velocity components, rather than the streamfunction (although the equation is listed in its general form above). The solver tolerance used for convergence was 10−8. The final balanced height field is then .
The equations are solved on a square domain of 600 km × 600 km. The setup is a zonal channel flow, with no-flux boundary conditions applied at the north and south lateral boundaries, and periodic boundary conditions applied at the west and east boundaries. The simulation is run for one revolution, so the final TC is located at the starting point, for comparison to the initial condition. For the AMR cases, adapting occurs when s−1.
The initial condition of case 1 is given in Fig. 2. Here the magnitude of the perturbation velocity vector for each simulation is shown in colored contours, with the elements overlaid. Only the element boundaries are overlaid and not the actual grid of nodal points inside each element. The uniform resolution mesh (FINE simulation) is shown in the top-left panel. In the top-right panel, the SMR2 simulation is shown, which has fully refined mesh for −200 y 200 km. The SMR1 initial condition is shown in the middle-right panel, which has fully refined mesh for −100 y 100 km. The initial conditions for the AMR3, AMR2, and AMR1 simulations are shown in the middle-right, bottom-left, and bottom-right panels, respectively. Here the mesh is adapting to the initial condition using the PV threshold providing a fully refined mesh around the hurricane eyewall (red region of stronger winds) and eye. The AMR3, AMR2, and AMR1 buffers are readily evident as layers of refined elements extending from the edge of the region of large PV.
In Fig. 3, a contour plot of the error is shown for the advection case at h, with the element boundaries at this time overlaid. For the SMR2 simulation, Fig. 3a shows that the largest errors exist at the static mesh northern and southern boundaries, in the intermediate-resolution mesh area. Similar results exist for the SMR1 simulation (Fig. 3b); however, the errors extend into the coarse mesh area. The region of larger errors corresponds to the location of where the vortex winds are forced to decay to zero by the blending function ( km). This blending is necessary to ensure the winds at the lateral boundaries are zero. Moving to Fig. 3c, the AMR3 simulation now has some more significant errors in the vortex inner-core region, as well as errors due to the coarser representation at the terminus of the radial decay function. In Figs. 3d and 3e, the AMR2 and AMR1 simulations show larger errors in the inner region, on the order of 10−2 and 10−1, respectively.
b. Case 2: Binary vortex interaction
The second case is a binary vortex interaction (Dritschel and Waugh 1992; Prieto et al. 2001, 2003). In this case, two TC-like vortices are offset by a certain distance and allowed to interact with one another. Depending on the offset distance, and the size and intensity of each vortex, different interactions can occur such as complete merger, complete straining out, partial straining out, and elastic interaction (Prieto et al. 2003). The case we choose here is a complete merger, where each vortex strains and elongates the other one, and they finally merge together into one central vortex.
The initial condition for the binary vortex interaction case consists of two offset Rankine vortices. Each vortex is constructed according to
where s−1 and km. The first vortex is positioned at km and the second vortex is positioned at km. The same as in case 1 is applied to each vortex to ensure the winds decay to zero prior to reaching the lateral boundary. No-flux boundary conditions are used at each lateral boundary. The nonlinear balance equation in (10) is solved using the initial condition of both vortices in order to obtain the corresponding field, only with s−1. For the AMR cases, adapting occurs when s−1.
The initial condition for the binary vortex interaction is given in Fig. 4. The two Rankine vortices are evident as the two PV maxima. All initial grids are able to resolve the initial PV anomalies well. Here the SMR2 simulation has a statically refined mesh between −200 (x, y) 200 km, and SMR1 simulation has a statically refined mesh between −100 (x, y) 100 km.
In Fig. 5, a close-up view of the evolution of the binary vortex interaction case is given. In Fig. 5a, the FINE simulation is shown in the inner region at the three times. In Figs. 5b–d, the AMR3, AMR2, and AMR1 simulations are shown, respectively. The close-up inspection confirms the previous results, showing the filaments are well resolved in the AMR3 and FINE simulations. It is worth noting that our PV adaptation criterion was s−1, which exceeds the filament magnitudes. The AMR1 simulation could have the PV filaments well resolved with a reduced threshold. However, in viewing these simulations as eventually being applicable to real-world TC prediction, adaptation would be done at a higher PV threshold as one would not want to adapt to any weak PV area on the globe, but rather stronger PV areas associated with tropical cyclones. Therefore, we feel that our method of using a larger PV threshold, with an adequate buffer region would be one way to adapt to evolving tropical cyclones (cf. the AMR3 simulation).
In Fig. 6, a contour plot of the error is given for the binary vortex interaction. In Figs. 6a and 6b, there are also errors at the terminus of the radial decay function for the SMR2 and SMR1 cases (more significant for SMR1 because of coarser representation of this area). Similar to the advection case, the inner-core errors begin to increase for the AMR3 to AMR1 simulations, with the AMR1 simulation having the largest inner region errors (Figs. 6c–e).
c. Case 3: Barotropic instability of the hurricane eyewall
The third case is the barotropic instability of the hurricane eyewall. The hurricane inner core can be described using a three-region model: (i) a low vorticity eye, (ii) a high vorticity eyewall, and (iii) a low vorticity environment. Observations of such vorticity structures in real hurricanes are given in Kossin and Eastin (2001) and Hendricks et al. (2012). An idealization of this structure can be constructed according to Schubert et al. (1999) using the tangential velocity profile:
which defines a discrete three region model of axisymmetric relative vorticity:
Here s−1, s−1, km, and km. The hurricane eye is defined as the region where , the eyewall is defined as the region between and , and the environment is defined as the region between and infinity. Similar to the previous experiments, a smooth radial decay function is multiplied by the tangential winds so that at . The cutoff radius km. No-flux boundary conditions are applied at each lateral boundary.
The nature of the instability is as follows. Each vorticity gradient of the ring supports a vortex Rossby wave (Montgomery and Kallenbach 1997). The inner vortex Rossby wave progrades relative to the mean flow, while the outer vortex Rossby wave retrogrades relative to the mean flow. Thus, it is possible for each of these waves to have the same angular velocity, or be phase locked, leading to the barotropic instability of the ring. A comprehensive linear stability analysis of this structure is provided by Schubert et al. (1999), and nonlinear simulations and discussions of aspects of this problem are given in Kossin and Schubert (2001) and Hendricks et al. (2009).
To initiate the instability process, a broadband perturbation was added to the basic-state vorticity equation (10) of the following form:
where s−1 is the amplitude and is the phase of azimuthal wavenumber m. For this set of experiments, the phase angles were chosen to be random numbers in the range . Each simulation was done using the same set of :
In real hurricanes, the impulse is expected to develop from a wide spectrum of background convective motions. The vorticity perturbation is inverted to produce a tangential velocity perturbation that is then added to the symmetric tangential velocity in (12). Similar to case 1, the nonlinear balance equation (10) is solved to obtain the corresponding field, with s−1. This ring has a thickness parameter and hollowness parameter (where is the average inner-core vorticity). According to the linear stability analysis of Schubert et al. (1999), this ring is most unstable to azimuthal wavenumber , with an e-folding growth time of 0.57 h. For the AMR cases, adapting occurs when s−1.
The initial condition of case 3 is given in Fig. 7. Here, a vortex with a very sharp gradient in tangential velocity is shown (ring of elevated PV). The initial condition of each simulation has a fully refined mesh over the ring of elevated PV. The AMR1 simulation (with no buffer) has a couple of coarse mesh cells at the very center as the initial mesh is adapting to the ring of large PV only.
A close-up view of the inner portion of the domain for this case is given in Fig. 8. In Fig. 8a, the FINE simulation is shown at different times moving downward. In Figs. 8b–d, the AMR3, AMR2, and AMR1 simulations are shown, respectively. The timing of the instability, mesovortex formation, mesovortex mergers, and final mixing into a monopole is identical in each simulation. Similar to the binary vortex interaction case, the weaker PV anomalies and filaments are better resolved in the FINE, AMR3, and AMR2 simulations than in the AMR1 simulation.
While the error results above are for the final time of each simulation, it is also insightful to examine the errors at intermediate times for case 3. In Fig. 9, the normalized and error norms are examined for h in the innermost 200 km × 200 km. In Figs. 9a and 9b, the errors in the velocity potential and geopotential are shown versus the total number of points. Again, in these panels, the AMR1 simulation is the leftmost circle, and the SMR2 simulation is the rightmost circle, with the other circles representing the other simulations. Broadly speaking, the results show that the error convergence is more rapid for early times than later times, and for each simulation the error grows with time. In Figs. 9c–e, contour plots of are shown for the AMR1, AMR2, and AMR3 simulations, respectively, for the magnitude of the velocity vector. The results show that error onset progression is more rapid for the AMR1 simulation than the AMR2 simulation, and the onset is more rapid for the AMR2 simulation than the AMR3 simulation.
Figure 10 shows the error for the instability case. Similar to the previous cases, there are increasing errors in the region of large PV moving from the SMR2 simulation (Fig. 10a) to the AMR1 simulation (Fig. 10e). This indicates that errors are reduced over the feature of interest by having fully refined elements in the vicinity.
d. Case 4: Barotropic instability of a shear zone
The fourth test is the examination of the formation of TC-like vortices through the barotropic instability of a region of large horizontal velocity shear. An example of such a process occurring in the real atmosphere is the instability of the ITCZ in the eastern North Pacific Ocean basin, causing it to undulate over days, and eventually break down forming pools of high PV (Fig. 2d). Provided favorable conditions, these PV pools can then form into TCs. An observational and modeling study of the genesis of TCs from this formation mechanism is given in Ferreira and Schubert (1997).
The initial condition is constructed as an idealization of the ITCZ, with easterly trade wind flow to the north and westerly trade wind flow to the south. A linear function is assumed to bridge the two regions, forming a thin strip of cyclonic vorticity. Mathematically, the initial condition is
where m s−1 and km. This specific width is used so that the zonal wavelength of the most unstable mode is an integer factor of the domain size in the zonal direction. Here, and is determined by solving the geostrophic balance equation analytically:
where s−1. The final height field is . While a β plane is more appropriate for this case, an f plane is used to obtain a first basic understanding of how the variable-resolution meshes simulate barotropic instability of this simple shear zone.
The strip of PV supports the existence of two counterpropagating Rossby waves: one on the northern PV gradient that propagates to the east and one on the southern PV gradient that propagates to the west (Heifetz et al. 1999). These two waves can phase lock and grow, leading to the breakdown of the strip into vortices. A linear stability analysis of this structure is given by Gill (1980), and here we provide a basic overview of the stability characteristics. Assuming separation of the meridional and zonal structure, Rossby wave solutions of the form are sought, where is the wave streamfunction and is an arbitrary meridional amplitude function. Here, is the complex phase velocity, and k is the zonal wavenumber. The linear stability analysis of this simple shear zone indicates the most unstable zonal wavenumber m−1, or 1600 km. The domain used here is 8000 km × 8000 km; therefore, the most unstable mode is zonal wavenumber 5. The growth rate s−1, corresponding to an e-folding time of approximately 3.5 h. To initiate the instability, a wavenumber-5 vorticity perturbation with amplitude of 1% of the basic-state vorticity is added for this case, and also inverted (in this case, integrated in the y direction) to produce a perturbation in the symmetric zonal velocity field in the region of constant background vorticity (shear zone). The lateral boundary conditions for this run are the same as case 1, no-flux conditions are applied at the north and south boundaries, and periodic conditions are applied at the west and east boundaries. For the AMR cases, adapting occurs when s−1.
The initial condition for case 4 is given for each simulation in Fig. 11. Decreasing amounts of refined elements are evident in moving from the FINE to the AMR1 case, with the AMR1 case only having refined elements over the PV strip itself. In Fig. 12, each simulation is shown at h. Each simulation produces the most unstable mode of wavenumber 5 (recall a weak wavenumber-5 perturbation was added to the initial conditions). Upshear tilt of each PV anomaly associated with the northern and southern counterpropagating Rossby waves is evident. The simulations are shown at h in Fig. 13. Here the breakdown of the PV strip has resulted in the five separate vortices. All simulations are able to reproduce the five vortices. In the AMR1 simulation, the vortices have begun to interact with one another, displacing them from the linear orientation. This interaction does not occur in the other simulations. This interaction happens because each vortex projects a larger-scale wind field, which is less well resolved in the AMR1 simulations than the others. As this simulation is executed for a long time (180 h), spatial and temporal discretization errors finally manifest in an earlier interaction of these five vortices. Therefore, the inclusion of a buffer region is very important in this case, particularly at the longer lead times.
Figure 14 shows a contour plot of the error for the ITCZ simulation. The results are similar to the other figures in that the errors increase in the region of fine variable-resolution mesh, particularly in the AMR1, AMR2, and AMR3 simulations (Figs. 14c–e). The AMR1 simulation has the largest errors due to the individual vortices being displaced from the linear orientation.
Here, we examine the efficiency of the variable-resolution mesh simulations in comparison to the FINE simulation. Each simulation was executed on a single central processing unit (CPU). In Fig. 15, the speedup of each simulation is shown versus the nodal point ratio. The speedup is defined as the ratio of the CPU time of the FINE simulation and the CPU time of each variable-resolution mesh simulation. The nodal point ratio is defined as the ratio of the number of nodal points in the FINE simulation and the fixed number of nodal points in the SMR simulations or the average number of nodal points in the AMR simulations. An ideal speedup versus point ratio would be the 1:1 line, which is shown in pink with dashes. Each test case is shown in a different color, and the AMR simulations have the largest point ratios and speedups, while the FINE simulation has a speedup to point ratio of 1 by definition.
All simulations were run with two compiler optimizations: O0 and O2. The compiler optimization O0 is the strictest level (also known as debug mode), while O2 optimizes various parts of the code for enhanced speedups. The results of O0 are shown in Fig. 15a. With optimization O0, the speedup is always less than the point ratio. For the SMR2 and SMR1 simulations, the speedups are closer to the point ratio. However, for the AMR cases, the speedup is significantly less than the point ratio. For example, the AMR1 simulation of the binary vortex case (red circle farthest to the right) has a speedup of approximately 7 for a point ratio of approximately 13. However, moving to Fig. 15b, with compiler optimization O2, the speedups obtained in this study are ideal, and approximately equivalent or slightly larger than the point ratio. These results indicate that the speedups are strongly sensitive to the compiler optimization, and that the superlinear speedups shown in Fig. 15b are a result of compiler optimization. All results in this manuscript were presented from simulations run with compiler optimization O2.
It is worth noting that these speedup results are obtained by comparing simulation results from the same code. It is not known what the speedups would be with regard to a fixed structured uniform code simulation at similar resolution, rather than the FINE simulation of this code with AMR switched off.
b. Error analysis
The normalized and errors are shown versus the average number of nodal points for each case and variable mesh simulation within each case in Fig. 16. In Figs. 16a and 16b, the errors in and are given, respectively, for all points in the domain. In Figs. 16c and 16d, the errors are shown for the localized region of the features that each variable mesh simulation attempts to resolve ( km for cases 1–3, and km for case 4). In each panel, the leftmost circle is the AMR1 simulation (fewest average number of nodal points), and the rightmost circle is the SMR2 simulation (most nodal points). The error is represented by the solid line and the error is dashed. Overall, each simulation shows the expected decrease in error with increasing nodal points. In Figs. 16a and 16b, the ITCZ error convergence is more rapid than the others. This may be due to lower numerical errors resulting from simulating a shear zone that is better aligned to the grid, instead of a vortex, or due to the fact that the balanced height field is specified for this case analytically. The AMR1 simulations have significantly larger errors than the other simulations, for both the geopotential and velocity vector magnitude. The geopotential error convergence curves are flatter than the velocity vector magnitude convergence curves, and they converge to higher errors. Examining the inner region where the active TC processes are occurring, there is a significant reduction in the errors from the AMR1 to the SMR1 simulation for velocity vector magnitude (Fig. 16c), while not much further reduction in errors moving to the SMR2 simulation. A similar progression is evident with regard to the geopotential (Fig. 16d); however, the AMR3 simulation has similar errors to the SMR2 simulation. Combining the results, AMR3 and SMR1 have significantly reduced errors over AMR1, and are a good trade-off between lower errors and computational expense (Fig. 15).
In Fig. 17, the same errors are shown versus the speedup. In each panel, the leftmost circle is now the SMR2 simulation (smallest speedup), and the rightmost circle is the AMR1 simulation (largest speedup). Each panel depicts the same result as in Fig. 16, only with regard to speedup instead of total points. In Fig. 17a, for cases 1–3, the slopes are positive indicating that the benefit of a faster simulation is hindered by higher errors in the velocity vector magnitude. The slopes are the error curves are more gradual for cases 1–3 than case 4. Examining Fig. 17b, the slopes are more gradual for the geopotential for cases 1–3, indicating that errors are not increased significantly by fewer refined elements near the feature of interest. Moving to the inner region of the domain (Fig. 17c), the slopes are larger for cases 1–3 in the errors for the velocity vector magnitude, and particularly large for case 4. The inner region geopotential error curves slopes are similar to the curves for the entire domain (Fig. 17d). For Figs. 16 and 17, the actual values for the average number of nodal points, point ratio, and speedup are given in Table 2.
In summary, the SMR2 and SMR1 simulations have the lowest errors overall, and low errors over most regions of the domain. The AMR1, AMR2, and AMR3 simulations have larger errors in the regions of interest. Overall the errors for velocity vector magnitude and geopotential were generally quite low and the AMR3 simulations appear to be a good trade-off between computational expense and reduced errors. This simulation significantly reduced errors in both velocity vector magnitude and geopotential in the region of interest over the AMR1 and AMR2 simulations. The error analysis confirms the utility of using a buffer region of fine elements that extends away from the region of high PV for these cases.
c. Energy conservation
The total energy conservation error for each case is given in Fig. 18. In each plot the domain-integrated value is shown at each time step over the full model integration. The total energy conservation error is defined as . Here, is the initial total energy and is the total energy at time t. Total energy is not conserved exactly due to the modal filter. For the advecting vortex (Fig. 18a), the final energy conservation error is approximately 5 10−7. The energy conservation error is lowest with the FINE simulation, and highest with the AMR1 simulation, indicating an inverse relationship with the total number of points. The binary vortex case (Fig. 18b) exhibits similar results; however, here the AMR1 case has a larger energy conservation error. For the unstable vortex case (Fig. 18c), the final energy conservation error for all cases is approximately 1 10−5. This case exhibited a cascade of energy to small scales during PV anomaly mergers and filamentation; therefore, it is expected that more energy would be dissipated in this simulation. The ITCZ case (Fig. 18d) exhibits a final energy conservation error similar to the unstable vortex. For this case a sharp spike in energy conservation error to lower values is seen at the onset of barotropic instability at h. In all cases, total energy is lost since is positive. Close inspection of the total energy conservation error was also performed in the early part of the simulation (first 10 and 100 time steps). For all simulations, the total energy conservation error was positive indicating that energy is lost initially.
d. Sensitivity of results to polynomial order
In the foregoing results, we have chosen a specific level of refinement based on the number of elements spanning the domain and the polynomial order. The specific level of refinement was chosen so that the feature of interest is well resolved. It is also insightful to examine the sensitivity of the results to the polynomial order N and the element size. There are three ways to examine more marginally resolved cases: (i) reduce N, (ii) increase the element size, or (iii) some combination of (i) and (ii). Here we elect to examine more marginally resolved cases using a reduction in the polynomial order from to .
The results of this sensitivity test are given in Fig. 19. Here we show the four cases we have examined here using the AMR2 simulations. In Fig. 19a, the advection case is shown at one-half revolution for and , respectively. There are no obvious differences in this case, signifying that the vortex advection is also reasonably well resolved at . In Fig. 19b, the binary vortex case is shown halfway through the simulation. Here, significant differences are seen, with the case having the wrong central vortex orientation and filaments that are much more diffuse. In Fig. 19c, the vortex instability case is shown at the two polynomial orders, respectively. In this case, running at causes the most unstable mode to be incorrect (azimuthal wavenumber 4). Finally, in Fig. 19d, the ITCZ case is shown at h. For this case the simulation produced the correct most unstable mode; however, at longer lead times, merging of the resultant vortices occurs, while in no merging occurs. Similar to the AMR1 versus AMR2 case at , this demonstrates that the larger numerical errors using are manifested in earlier merging of the vortices.
In summary, the polynomial order was reduced from to for each case for the AMR2 simulations. The reduction to leads to a more marginally resolved feature of interest. The strategy employed in this work was to design the refinement so that the local feature of interest is resolved well. These same cases were also examined at and the results were similar to indicated that the results are reasonably well converged at .
A two-dimensional f-plane shallow-water model based on the continuous Galerkin (spectral element) numerical method has been used to examine idealized tropical cyclone (TC) problems, with a focus on the applicability of static and adaptive mesh refinement (SMR and AMR, respectively). Four different idealizations of TC cases in the real atmosphere were simulated in this model, with varying degrees of SMR and AMR. The SMR and AMR simulations were compared to a high-resolution “truth” simulation (noted previously as the FINE run) with regard to solution accuracy and computational time. Three different AMR simulations were conducted with varying levels of buffer regions, or the number of extra layers of fine elements added to the finely resolved region. Two different SMR simulations were executed with varying amounts of refined elements. For the AMR simulations, a normalized potential vorticity threshold was used for refining and de-refining the elements.
With regard to solution accuracy of the local feature of interest, the SMR simulations overall had lower errors than the AMR simulations. However, the AMR3 simulation, with a buffer of six fine elements away from the PV anomalies, was shown to have only a slightly larger error than the SMR1 simulation. The AMR1 and AMR2 simulations (with a smaller buffer region) had larger errors. Significant speedups were obtained by using AMR. The AMR3 simulation was 4–8 times faster than the FINE simulation (at a minor error expense), while the AMR1 simulation had speedups of up to a factor of 15. Overall, these results indicate that AMR can be used at significantly less computational expense to resolve the TC feature nearly as well as the truth run, provided a sufficient buffer region is used. We have also demonstrated the utility of a potential vorticity criterion for mesh adaptation in conjunction with a buffer region of arbitrary layers of fully refined elements in order to track tropical cyclones, and individual features within tropical cyclones, in an idealized framework.
In closing, we wish to note that we have examined static and dynamic adaptive mesh refinement for TC applications in an idealized framework of a shallow-water fluid in constant rotation. In the real atmosphere, TCs are three-dimensional phenomena, with complex physics interactions (microphysics, boundary layer, vertical mixing, and radiation), as well as interactions with the environment (such as vertical wind shear and ocean surface fluxes). Additionally, one of the major challenges in the future with AMR is the development of scale-aware physical parameterizations that will seamlessly represent physical processes across scales. For the reasons above, the utility and applicability of these variable-resolution grids to real TC prediction is still an open question, and more research is clearly needed. For example, a PV threshold has been used here for refinement; however, using one arbitrary PV threshold would not work for a typical real-world situation where there are many different types and scales of features of interest where one may wish to do local refinement. However, these results demonstrate that from a purely dry dynamical modeling standpoint in two dimensions, AMR shows great promise for TC simulations. The results here could be readily extended to examine more complex ideal TC scenarios, such as dry two-dimensional and three-dimensional simulations (with SMR and AMR only in the horizontal) on the sphere. Many aspects of tropical cyclones (such as physical mechanisms governing intensity and structure change) are not well understood. These ideal simulations using SMR and AMR could provide very high resolution over the inner core (eye and eyewall) to resolve important aspects of intensity variability, while using coarser resolution in the outer wind field so that these simulations would be computationally feasible.
EAH, MSP, JDD, and QJ gratefully acknowledge support from the Chief of Naval Research PE-0601153N. EAH also acknowledges support from Naval Postgraduate School Research Initiation Program. The codes used in this paper were developed by MAK and FXG during the 4-month long program on Multiscale Numerics in the Atmosphere and Ocean at the Newton Institute, Cambridge University. MAK and FXG gratefully acknowledge the support of the Office of Naval Research through Program Element PE-0602435N, the National Science Foundation (Division of Mathematical Sciences) through Program Element 1216700, and the Air Force Office of Scientific Research through the Computational Mathematics program as well as the ONR Computational Mathematics Program. We thank Sasa Gabersek, Alex Reinecke, and Kevin Viner for helpful discussions. This manuscript was improved by the helpful comments of three anonymous reviewers and Dr. Hilary Weller.