The Effects of Numerical Dissipation on Hurricane Rapid Intensification With Observational Heating

The computational fluid dynamics of hurricane rapid intensification (RI) is examined through idealized simulations using two codes: a community‐based, finite‐difference/split‐explicit model (WRF) and a spectral‐element/semi‐implicit model (NUMA). The focus of the analysis is on the effects of implicit numerical dissipation (IND) in the energetics of the vortex response to heating, which embodies the fundamental dynamics in the hurricane RI process. The heating considered here is derived from observations: four‐dimensional, fully nonlinear, latent heating/cooling rates calculated from airborne Doppler radar measurements collected in a hurricane undergoing RI. The results continue to show significant IND in WRF relative to NUMA with a reduction in various intensity metrics: (a) time‐integrated, mean kinetic energy values in WRF are ∼20% lower than NUMA and (b) peak, localized wind speeds in WRF are ∼12 m/s lower than NUMA. Values of the eddy diffusivity in WRF need to be reduced by ∼50% from those in NUMA to produce a similar intensity time series. Kinetic energy budgets demonstrate that the pressure contribution is the main factor in the model differences with WRF producing smaller energy input to the vortex by ∼23%, on average. The low‐order spatial discretization of the pressure gradient in WRF is implicated in the IND. In addition, the eddy transport term is found to have a largely positive impact on the vortex intensification with a mean contribution of ∼20%. Overall, these results have important implications for the research and operational forecasting communities that use WRF and WRF‐like numerical models.

stagnant (Marks & Shay, 1998;Rappaport et al., 2009), with large forecast errors still present today (DeMaria et al., 2014). As described in DeMaria et al. (2014), none of the deterministic models had RI forecasting capability from 1991 to around 2015. There has been some ability to forecast RI with both dynamical and statistical models since 2015, however, a significant under-prediction or low bias in RI cases is still present today (DeMaria et al., 2021).
The increase of kinetic energy in hurricane intensification, as well as RI, is driven by the vortex response to heating in convective clouds (e.g., Shapiro & Willoughby, 1982), where the source of moist enthalpy flux comes from the thermodynamic disequilibrium between the ocean and atmosphere, (e.g., Emanuel, 1986). Dissipation of energy occurs most prominently in the boundary layer through surface friction and a hierarchy of turbulent eddies of various scales. However, new research has shown that the hurricane boundary layer is not purely dissipative and contains coherent turbulent structures that can "backscatter" energy to larger scales (Sroka & Guimond, 2021). In numerical models, dissipation of energy can also occur implicitly through the algorithms used to solve the fluid-flow equations ("implicit numerical dissipation [IND]") or explicitly through the addition of filters. Implicit numerical dissipation can occur from the use of low-order discretizations of the governing equations in both space and time. For example, the use of second-order or upwind-biased spatial discretizations of the advective terms can result in substantial numerical dissipation error when compared to high-order (e.g., fifth) or centered schemes (Hoffman & Frankel, 2001;Skamarock & Klemp, 2008).
In general, minimal numerical dissipation is desired in highly nonlinear computational fluid dynamics problems, such as hurricanes, because errors incurred from excessive dissipation can quickly propagate through the system. Kravchenko and Moin (1997) examined numerical errors in spectral and finite difference codes as well as the effects of sub-grid scale models in turbulent channel flow. They demonstrated that the high wavenumber portion of the energy spectrum is severely damped by truncation errors in low-order (e.g., second) finite-difference schemes, and the contribution of the sub-grid scale model is small in this context. By increasing the order of the finite-difference approximations, the results of their large eddy simulations and the performance of the sub-grid scale model were enhanced. Larsson et al. (2007) found that maintaining low numerical dissipation was important for simulating shock/turbulence interactions, especially for coarse resolution simulations where the fields are under-resolved and a sub-grid model is required. This setup is also the case in mesoscale atmospheric modeling, which is our focus in this paper. In the under-resolved simulations of (Larsson et al., 2007), the numerical dissipation was large enough to dampen or erase the smaller-scale motions on the grid and from the sub-grid model. However, not all numerical dissipation is detrimental to the success of a simulation as some amount is required for stability and it has been shown that for large eddy simulations, IND can mimic the effects of a sub-grid turbulence model (Implicit Large Eddy Simulation: Computing Turbulent Fluid Dynamics, 2007). Furthermore, the use of some numerical dissipation (as well as artificial viscosity) can help control the effects of numerical dispersion errors, which occur in the high wavenumber portions of the flow field where the dissipation is active.
Continuous Galerkin (CG) and discontinuous Galerkin (DG) numerical methods have several unique properties which distinguish them from low-order (i.e., the order of accuracy equal to or smaller than two) methods and other high-order methods, such as finite volume/difference schemes. These include: (a) possessing low dissipation and dispersion errors for turbulent flows with highly disparate spatial and time scales; (b) achieving arbitrary high-order discretization for all spatial derivatives; and (c) highly scalable, and efficient on massively parallel supercomputers, such as those accelerated by Graphics Processing Units (GPUs) (Abdi et al., 2017(Abdi et al., , 2019. These superior numerical properties make high-order CG and DG methods attractive for hurricane research. The advantages of high-order numerical methods over their low-order counterparts for low Reynolds number flow problems has been demonstrated through a workshop series, The International Workshop on High-Order CFD Methods (Wang et al., 2013). However, when simulating under-resolved problems that require a sub-grid turbulence model and problems with discontinuities, high-order numerical schemes can have issues with excessive grid scale noise, aliasing and stability (Gassner & Beck, 2013;Honein & Moin, 2004;Moura et al., 2017). These issues can lead to errors in the simulated flow or a failure of the simulation. To address these potential problems, de-aliasing techniques (Blaisdell et al., 1996;Fischer & Mullen, 2001;Gassner, 2013;Gassner & Beck, 2013;Karamanos & Karniadakis, 2000), localized artificial viscosity (Persson & Peraire, 2006;Yu et al., 2015) and limiters (Cockburn & Shu, 1998;Qiu & Shu, 2005;Zhang & Shu, 2010) can be used to stabilize flow simulations. In this study, we rely on a combined approach that applies artificial viscosity based on output from a turbulent kinetic energy (TKE) sub-grid model for turbulent diffusion; see details in Section 2.3. Takemi and Rotunno (2003) studied the effects of sub-grid mixing and numerical filtering in squall line simulations using the Weather Research and Forecasting (WRF) model, which is a finite difference based code that can provide high order discretization for the advective (or flux divergence) terms only (Skamarock & Klemp, 2008). The authors found that using a fifth-order, upwind-biased advection scheme along with a standard TKE sub-grid model resulted in many noisy, grid-scale convective cells. Rather than applying an explicit numerical filter, which could damage the physical modes, the authors tuned the TKE sub-grid model coefficient to produce reasonable convective structures and energy scaling with wavelength. They also tested the inclusion of an explicit numerical filter and found that it had a much larger effect on the solutions than the sub-grid TKE model, which highlights the importance of analyzing numerical dissipation. It is clear that in order to ensure both high accuracy and stability of a simulated flow, a careful balance between signal and noise must be achieved.
In a theoretical hurricane study, Guimond et al. (2016) showed that the vortex response to simple, impulsive, asymmetric thermal anomalies can produce significant differences in system intensity across models due to the amount of IND. The community atmospheric model (WRF) was shown to have anomalously large IND when compared to research atmospheric codes (the High-Gradient [HIGRAD] model and the Non-hydrostatic Unified Model of the Atmosphere [NUMA]), which resulted in a muted intensity response from asymmetric thermal anomalies. We note that a community atmospheric model is created and used by a large and diverse group of scientists/engineers using well-established methods. The model is designed to be versatile and user-friendly for a wide range of users. A research atmospheric model is a model created and used by a small and specialized group of scientists/engineers using cutting-edge methods. This type of model necessitates extensive knowledge of computational methods.
The HIGRAD and NUMA models produced a more energetic response due to much less numerical dissipation. Spectral kinetic energy budgets showed that the pressure gradient term was the dominant source of the anomalous dissipation in WRF with the flux of inertia-gravity wave energy describing most of the variance in the pressure term. Acoustic and inertia-gravity waves are considered fast modes in the equation set, which are split off from the slow modes in WRF. This understanding lead to the recommendation that the time integration scheme was the main culprit for the numerical dissipation in WRF. Evidence for this hypothesis was shown through sensitivity tests with different time integration schemes in NUMA, which showed significant differences in the amount of energy and role of the pressure term.
In this work, we study the response of a tropical storm-like vortex to time-dependent, 3-D observational heating calculated from airborne Doppler radar measurements in the RI of Hurricane Guillermo (1997). The remainder of the paper is organized as follows. In Section 2, a detailed description of the numerical models and simulation setup is presented. Therein, we introduce the WRF and NUMA models, vortex initialization and heating strategies, and eddy viscosity and diffusivity setup. Comparison of the wind field features, for example, maximum and azimuthal mean wind speed as well as kinetic energy, from WRF and NUMA is discussed in Section 3. In this section, we also present kinetic energy budget analyses to explain the wind field disparity between WRF and NUMA. Similar features for longer time period simulations are studied on Section 4. Important implications of this work in the hurricane research and operational fields are given in Section 5. Future work is also discussed in this section.

Description and Setup of Numerical Models
A comprehensive introduction of the governing equations and numerical methods used in the WRF-ARW (hereafter WRF) and NUMA models has been given in Guimond et al. (2016). For completeness, we review them below and provide a quick reference of the numerical methods, boundary conditions and sub-grid scale diffusion settings used in each model in Table 1.

The WRF Model
The WRF model solves the compressible, non-hydrostatic Euler equations in a conservative form with a η mass vertical coordinate (Laprise, 1992;Skamarock & Klemp, 2008). For comparisons with NUMA, all variables are interpolated to regular height levels in post-processing. Note that the differences between the η levels and height are very small for these idealized simulations. The simplified model equations for a dry atmosphere can be expressed as follows, using a Laplacian operator for explicit diffusion and η as the vertical coordinate: 10.1029/2021MS002897 4 of 25 Here u, v and w are the velocities in three dimensions, m = m(x, y) is the mass per unit area within a column, θ is the potential temperature, ρ is the dry air density, is the perturbation pressure, f is the Coriolis parameter, g is gravity, μ is the eddy viscosity, κ is the thermal diffusivity, S is the heating rate source term, and ∇ is the three-dimensional gradient operator. Variables with a hat denote perturbations from the hydrostatically balanced reference state.
A combined finite-difference/finite-volume spatial discretization of the governing equations is employed in WRF. In the horizontal and vertical directions, a spatially staggered Arakawa C grid is utilized where velocities are defined on the cell faces and scalars at the cell centers. For the nonlinear advective terms, a fifth-order, upwind-biased discretization in the horizontal and a third-order scheme in the vertical are typically used. We have utilized these settings here, but also tested the impacts of the less diffusive, even-ordered schemes (sixth-order and fourth-order in the horizontal and vertical dimensions, respectively). The differences between the even-ordered and odd-ordered schemes were small (maximum values of ±0.5 m/s in the eyewall) and therefore, we utilize the odd-ordered formulations in all presented results. The WRF model relies on a split-explicit time integration process, where acoustic and gravity wave modes are calculated using a small time step and advection is computed on a larger time step (Klemp & Wilhelmson, 1978;Skamarock & Klemp, 2008;Wicker & Skamarock, 2002). Horizontal modes are solved explicitly within the small time stage, while vertical modes are implicitly solved. The implicit solve is done with backward Euler. A third-order Runge-Kutta scheme is used to perform the overall time integration, including both the small-and large-time step equations. The small time step results are applied as a correction to the large time step calculations during the Runge-Kutta time integration. More details on WRF can be found in (Skamarock & Klemp, 2008). Finally, we seek to produce minimally dissipative WRF solutions and therefore, we have turned off all filtering/damping options: explicit sixth-order numerical filtering, vertical velocity damping, divergence damping and external mode damping. Artificial viscosity is applied at the model top and through the Laplacian operators in the above equations are discussed further in Section 2.3.

The NUMA Model
The NUMA model is capable of using various forms of the Euler equations (for example, Giraldo and Restelli, 2008;Giraldo et al., 2010). However, for this study, we use the non-conservative form using potential temperature as the thermodynamic variable (Kelly and Giraldo, 2012;Giraldo et al., 2013) to be consistent with Guimond et al. (2016). The choice of conservative or non-conservative equation set is not expected to make a significant difference because the error resulting from the non-conservative set is much lower than the temporal error (Giraldo & Restelli, 2008). Instead of the η mass vertical coordinate, physical height z is used in NUMA.
The governing equations are expressed as: The spatial discretization of Equations 6-10 is carried out using the CG spectral element method (CG-SEM) (Giraldo & Restelli, 2008;Giraldo et al., 2013). Specifically, the physical domain is decomposed into a set of non-overlapping hexahedral elements and inside each element, the state variables are represented by polynomial expansion using Lagrange basis functions of a chosen order. The continuous spatial derivatives are constructed in discrete form by analytically taking derivatives of the polynomials that approximate the solutions. The state variables in each element are collocated with each other and placed at unequally spaced Legendre-Gauss-Lobatto points. In this study, we utilize fifth-order polynomial basis functions in all three spatial dimensions, which also provides fifth-order accuracy for all spatial derivatives and is identical to that presented in Guimond et al. (2016). Note that the stencil for all polynomial orders in NUMA is symmetric about the element centroid, so upwind-biased diffusion for fifth-order polynomials is not present. For time integration, the three-dimensional semi-implicit methodology (Giraldo et al., 2013) is used along with a second-order leapfrog scheme (LF2) with backward Euler off-centering. A first-order Robert-Asselin time filter is applied to stabilize the LF2 scheme with a coefficient of 0.2, which effectively removes the computational mode. The above description of NUMA comprises our control simulations. Several other time integration formulations are available in NUMA and we will note where sensitivity tests have been conducted. Interested readers are referred to Giraldo and Restelli (2008); Giraldo et al. (2013) for more details of the NUMA model.

Details of Simulation Setup
Careful analysis has been undertaken to setup WRF and NUMA nearly exactly the same to isolate any differences in the model solutions to the numerical schemes that comprise the dynamic core. In addition, as described below, we have imposed a number of idealizations in the physics in order to isolate the response of the vortex to heating, which is the fundamental dynamics controlling hurricane intensification. Our approach is to study hurricane intensification using a model hierarchy whereby idealized simulations with different numerical methods are used to build a base of knowledge and then layers of physics are added to carry this knowledge forward to the "full physics" regime. Therefore, we do not intend to reproduce the details of a specific case.
The computational domain is a box extending 800 km in the horizontal directions and 20 km in the vertical direction. In WRF, 2 km grid spacing in the horizontal is chosen to match that of the radar observations used as forcing and to be consistent with Guimond et al. (2016). The first model level is found at 167 m with constant vertical spacing of 333 m up to the model top (60 levels). To match the horizontal and vertical grid spacing in WRF, we have used 80 elements in each horizontal direction and 12 elements in the vertical direction along with fifth-order polynomials in all dimensions for the NUMA grid, as described in Section 2.2. These settings yield an element-averaged grid spacing in NUMA of ∼2 km in the horizontal and ∼333 m in the vertical. A time step of 2 s is used in each model.
Periodic boundary conditions are used in both horizontal directions in each model. A gravity wave absorbing zone (sine-squared function) is imposed at the top of the computational domain with the WRF zone occupying the top 4 km along with a small coefficient (0.00833) and the NUMA zone representing the top 1 km with a large coefficient (1.0). The differences in the absorbing zones are due to stability issues and sensitivity tests show the results are not sensitive to these choices. The free-slip boundary condition is applied at the bottom of the computational domain in each model, which disables fluxes of quantities (such as heat) from the surface and prevents a frictional boundary layer from developing. These idealizations enable the focus to be on the vortex dynamic response to the imposed heating. The simulations are run without moisture, but instead, four-dimensional latent heating/cooling rates derived from airborne Doppler radar observations are used to force the model as described below.
In post-processing, both WRF and NUMA fields are interpolated to a uniform, collocated grid at the horizontal/vertical grid spacings listed above. Linear interpolation is used to post-process the WRF results. To post-process NUMA results, a high-order interpolation based on Lagrange polynomials is applied, which is facilitated by the spectral element method (since by construction the solution exists everywhere in the element). To ease interpolation, any hexahedral element in the physical space (x, y, z) can be transformed into a standardized space Thus, for any Nth-order standard hexahedral element, there are N + 1 Legendre-Gauss-Lobatto points, namely, (α i , β i , γ i ), i = 1, …, N + 1, in each direction α, β and γ. The Lagrange polynomial basis L IJK (α, β, γ), I, J, K = 1, …, N + 1, can be constructed using the tensor product as Then, any value of the flow variables V(α, β, γ), such as the wind speed, inside a hexahedral element can be interpolated from the solutions V IJK , I, J, K = 1, …, N + 1, on the Legendre-Gauss-Lobatto points as The initial conditions are identical in each model. The hydrostatic and gradient-wind balanced, tropical storm-like vortex with axisymmetric tangential winds described in Equation 10 of Guimond et al. (2016) is utilized here, which is similar to the study of Nolan and Grasso (2003). The tangential velocity field in the radius-height plane for this initial vortex is shown in Figure 1. On top of this vortex, four-dimensional latent heating/cooling rates derived from airborne Doppler radar measurements in rapidly intensifying Hurricane Guillermo (1997) described in Guimond et al. (2011), are added as a source term in the potential temperature equation in both WRF and NUMA. These heating fields are computed on a grid covering the inner-core of the system out to a radius of ∼60 km with a grid spacing of 2 and 0.5 km in the horizontal and vertical dimensions, respectively. There are 10 heating snapshots covering ∼5.7 hr with a time step of ∼34 min. The peak heating is located at a radius of 25-30 km, which is well inside the radius of maximum wind (RMW) of the initial vortex (∼50 km). This heating and vortex configuration represents the RI process well, where convective bursts are the main driving force, for example, (Guimond et al., 2010;Montgomery et al., 2006;Reasor et al., 2009;Rogers et al., 2013). Guimond et al. (2011) conducted an extensive uncertainty analysis of the latent heat retrievals and found they were reasonably accurate, especially for convective bursts with randomly distributed errors in the heating magnitudes of ∼16% for updrafts greater than 5 m/s. In addition, Guimond and Reisner (2012) inserted the heating retrievals into realistic forecasts of Guillermo and found very good agreement in the predicted wind fields relative to observations.
Starting from the initial conditions, the first heating snapshot is introduced into the model over a 30 min period using a hyperbolic tangent function. Then, the remaining heating snapshots are linearly interpolated to the next  Figure 2 shows the three-dimensional structure of the heating for three snapshots and the time evolution function used to control the forcing in the models. Note that we have also added an exponentially decaying function at the upper-edge of the observational domain (10 km) to smoothly transition the data into the model grid, which helps maintain numerical stability.
Both models are also supplied the same explicit diffusion settings. While we can utilize the same sub-grid scale turbulence scheme in WRF and NUMA for our comparisons, the differences in the dynamic core of each model and associated dissipation characteristics will likely produce different eddy viscosity values during the course of the simulations. To isolate the model differences to the numerical formulation, we developed a simple height-dependent eddy viscosity model based on output from the WRF 3D TKE sub-grid turbulence scheme. Initially, we conducted a WRF simulation with the 3D TKE scheme to get an idea of the eddy viscosity and diffusivity values produced from the vortex and heating. Following a parcel, the sources and sinks of TKE in this scheme depend on the shear, buoyancy and dissipation. Details describing the implementation of this scheme in WRF, including the parameterization for dissipation, can found in Skamarock et al. (2021). The observational heat forcing will generate TKE from both the buoyancy and shear terms, but we only focus on the output eddy viscosities and diffusivities, which are calculated as where e is the TKE, C k is a constant of 0.15, and l is a length scale, which is around 2,000 m in the horizontal and 375 m in the vertical. Figure 3 shows plots of the horizontal eddy viscosity from the 3D TKE scheme at 0.50 and 9.80 km height at 6 hr. The eyewall is visible in the figures with viscosity values of ∼240 m 2 s −1 or larger in a thin ring at 0.5 km height and a broader region of 500-750 m 2 s −1 values at 9.80 km height.
Localized regions of higher viscosity values near 400 m 2 s −1 and 1,500 m 2 s −1 at lower and upper levels, respectively, are connected to large, vertically coherent heating pulses from convective bursts. Note that we have set the turbulent Prandtl number in WRF, which has a default value of 1/3, to 1 which enables the same eddy viscosity/ diffusivity values for momentum and scalars. There are areas of the WRF software where the default value is hard coded and we have taken careful steps to maintain values of 1 throughout the code. The turbulent Prandtl number is set to one in NUMA as well. Figure 4 shows the maximum horizontal and vertical eddy viscosity values produced over the 6 hr WRF simulation as a function of height. Both curves have relatively small values at lower levels, but they increase sharply at middle levels with some additional oscillations up to ∼16 km height. The large values found at middle to upper levels are from the strong heating pulses associated with convective bursts as seen in Figure 3. The values drop off sharply at ∼16 km height because that is where the gravity wave sponge is introduced into the model. The maximum horizontal viscosity values are about five times larger than the vertical values. Overlaid on top of the maximum viscosity curves are high-order polynomial fits that approximate the general structure and values of the eddy viscosity data. These fits take the following form where z is the geometric height in meters and the coefficients of the fifth-order polynomial are found in Table 2.
To stabilize the NUMA code, an extra constant was added to the ̃ coefficient for both the horizontal and vertical polynomial fits. The ̃ coefficient is shown in Table 2 and the curves in Figure 4 include this offset. Finally, these height-dependent eddy viscosity values are used in both WRF and NUMA for momentum and scalar diffusion in the comparison simulations. This simple explicit diffusion model is intended to both stabilize each numerical model and also represent, to some degree, realistic sub-grid scale turbulent diffusion from the TKE scheme.

Time Series and Windspeed Structure
In this section, we compare the solutions from the WRF and NUMA control simulations in terms of time series of horizontal wind speed and kinetic energy as well as the structure of windspeed perturbations. Here, perturbation is defined as the total wind speed at a particular time minus the wind speed of the initial condition, which helps identify the wind structures produced from the observational heating. Figure 5 shows the maximum windspeed output every 30 min for the control WRF and NUMA simulations with solid red and green lines, respectively. The maximum winds increased by about 45 m/s in 6 hr, which is a very large RI rate. The reason for this high rate, besides the idealized setup, is the much larger (and weaker) initial vortex compared to that which occurred in nature, which drives a large inward movement of angular momentum and associated increase in winds. Guimond and Reisner (2012) considered the same observational heating as the present study, but used an initial vortex based on radar observations of Hurricane Guillermo (1997) that had a much smaller RMW (∼30 km) compared to the vortex used here (∼50 km). Guimond and Reisner (2012) found that the minimum pressure dropped by about 12-15 hPa in 6 hr with the realistic initial vortex, which compared more favorably with observations than the present vortex. Nevertheless, the goal of this study is to analyze the idealized vortex response to heating pulses derived from observations in a RI system and examine the effects of IND in this process. We do not intend to simulate and reproduce the Guillermo case study. Thus, the current initial vortex and model setup are sufficient for the goals of this work.
Figure 5 also shows that after ∼2 hr, the NUMA winds begin to increase relative to WRF and during the last couple of hours of the simulations, the maximum windspeed is 2-7 m/s or 4%-12% higher in NUMA compared to WRF. A time series of minimum pressure perturbation (not shown) was also analyzed and the evolution is similar to Figure 5 with NUMA producing minimum pressures 5.5-6.5 hPa lower than WRF over the last several hours of the simulation.
In an attempt to more closely match the time series of WRF and NUMA, three sensitivity tests were conducted with WRF where the eddy viscosity values were set to 25%, 50% and 75% of the default values. The 75% tests still show significantly reduced maximum winds relative to NUMA, while the 25% tests generally seem too high, especially before 3.5 hr. In general, the 50% tests show a much closer match to NUMA, especially up to and including 3.5 hr, despite the anomalously high value at 2.5 hr. There is some larger variability between 4 and 6 hr, but smoothing through that variability indicates a reasonable match to NUMA. Therefore, these results indicate that in order to produce a similar intensity time series to NUMA, the explicit diffusion in WRF must be turned down significantly, with a reduction in eddy viscosity values of ∼50% relative to those in NUMA.
Additional time series diagnostics for azimuthal mean quantities were also calculated. While the environment surrounding the vortex has no mean flow, the observational heat forcing has an azimuthal wavenumber-one structure as can be seen in Figure 2, which produces a wavenumber-one flow asymmetry that slowly moves the vortex to the southeast. The storm center is computed through a simple iterative method that finds the position which maximizes the azimuthal mean windspeed. The data are interpolated onto a cylindrical grid with this storm center and the azimuthal mean quantities are calculated. Figures 6 and 7 show the times series of maximum azimuthal mean windspeed and mean kinetic energy, respectively. These figures show a similar qualitative pattern as the maximum windspeed with NUMA producing larger azimuthal mean windspeeds and kinetic energy values relative to WRF, with those differences growing over time. In the last couple of hours, the maximum azimuthal mean windspeed is about 4 m/s or 8% higher in NUMA compared to WRF. The mean kinetic energy follows a very similar pattern to  the windspeed and will be used as a reference for a dynamical budget analysis, presented in Section 3.2, to explain the reasons behind the model differences. Figure 6 also shows the maximum azimuthal mean windspeed for WRF simulations at 1 and 0.5 km horizontal resolution. For these simulations, everything is the same as the WRF control run with the exception of the horizontal grid spacing, which allows one to examine how the IND will change with the addition of a finer mesh. There is a noticeable increase in the WRF mean windspeed for the 1 and 0.5 km runs relative to the 2 km control case between 3 and 6 hr with peak differences of ∼1 m/s. However, there is very little, if any, differences between the WRF 1 and 0.5 km simulations for the maximum azimuthal mean  windspeed metric. These results indicate that some of the IND in WRF is reduced through increasing the horizontal resolution by at least a factor of two, which creates a slightly stronger vortex that is more like NUMA. Further grid refinements, especially in the vertical (not examined), could produce additional increases in intensity. The details on why WRF is producing a lower intensity response relative to NUMA will be described in Section 3.2.
While the differences between WRF and NUMA shown in Figures 6 and 7 are not large for these mean quantities, there is more substantial variability on local space and time scales, which is demonstrated next. In addition, it is important to keep in mind the short time period of these simulations (dictated by the available observations) and the idealized nature of the setup, both of which will limit the variability in the models.
At 5 and 6 hr, the vortex has reached peak intensity with perturbation windspeeds of ∼60 m/s found on the Northeast and Northern side of the storm as shown in Figure 9. At these time periods, the RMW of the vortex is 15-20 km with NUMA on the lower side and WRF on the higher side of that interval. The RMW of the initial vortex was ∼50 km and this large, rapid contraction rate is consistent with the rapid increase in winds from conservation of angular momentum.
At 5 hr, NUMA shows significantly stronger winds than WRF by 7-8 m/s within large portions of the northern eyewall including peak differences of up to +12 m/s. Similar to the previous time period, the low-wind eye of NUMA is a bit wider than WRF, which produces large negative differences in the eye and a larger radial wind gradient. Some thin bands of higher wind differences can also be seen to the North and Northeast of the storm, which may be related to vortex Rossby wave dynamics. At 6 hr, the local wind differences are smaller, but still significant in that much of the eyewall has positive differences of ∼5 m/s with larger values on the Eastern side of the vortex. Figure 10 shows the windspeed differences (NUMA-WRF) at 0.19 km height at 4 hr, 5 and 6 hr as a function of WRF grid spacing. The differences described above at 2 km WRF spacing are shown again in Figure 10 in the first row for ease of comparison. It is clear from Figure 10 that going from 2 to 1 km spacing in WRF reduces the differences with NUMA and creates a stronger vortex in the eyewall region by ∼2-3 m/s for all time intervals. The largest changes tend to occur in localized regions, such as the Northern portions of the eyewall in Figures 10b and 10e, as well as the eye in these same figures. Transitioning from 1 to 0.5 km WRF simulations does not produce significant changes with only ∼0.5 m/s reductions in windspeed differences across all time intervals. These results are consistent with the maximum azimuthal mean windspeed intensity metric shown in Figure 6.

Budget Analyses
The previous section established clear differences between the two numerical models with NUMA producing a larger intensity response. How can this be given that each model was set up exactly the same with the same initial conditions and same heat forcing? The answer to this question lies in the design of the numerical schemes that make up the dynamic core of each model and in this section, we analyze how the intensity differences are produced and highlight parts of the numerical scheme that are driving this effect.
The horizontal kinetic energy for the azimuthal mean vortex in cylindrical coordinates (r, θ, z) is expressed as where u is the radial windspeed, v is the tangential windspeed and the overbar indicates an azimuthal mean quantity. After azimuthally averaging, these variables and those below are functions of radius (r) and height (z) unless noted otherwise. After multiplying the radial and tangential equations of motion by their corresponding velocity, summing the two equations and applying Reynolds decomposition in the azimuthal direction (the over bar and prime notations below indicate azimuthal mean and eddy variables, respectively), we arrive at the transport equation for azimuthal mean kinetic energy,̄= where, The full derivation of this equation can be found in the Appendix A. In Equation 16, M defines the mean kinetic energy transport terms, E defines the eddy transport terms which represent the Reynolds stress contributions in the azimuthal dimension, P defines the pressure gradient term and D defines the total explicit diffusion term. Figure 11 shows a times series of azimuthal mean kinetic energy budget tendencies from both WRF and NUMA after averaging the fields over the eyewall (∼10-50 km radius) and height (∼0.19-1.5 km). The largest term is the pressure gradient, which contributes positively to the increase in mean kinetic energy of the vortex shown in Figure 7. This large positive contribution is from the input heating, which leads to significant integrated warming in the storm core and an associated increase in the radial pressure gradient between the undisturbed outer regions and lowered pressures in the core region. This larger radial pressure gradient drives strong inflow at low levels, which transports high angular momentum from the outer regions of the storm toward the center, increasing the tangential velocity of the vortex. Large differences between WRF and NUMA are present in the pressure gradient term, especially beyond 2 hr with NUMA larger than WRF by ∼25% on average with a maximum of ∼40%.
The second largest contribution is the mean transport term, which shows largely negative values that increase with time as the mean flow intensifies. The vertical mean transport dominates over the horizontal transport, which is dictated by the heating profile that is maximized near middle levels. Therefore, there is a significant positive flux of kinetic energy out of the lower levels of the vortex, which results in a net sink of energy. However, the differences between WRF and NUMA are much smaller for the mean transport term, relative to the pressure gradient, by a factor of ∼4.5. The eddy transport term oscillates around zero tendency up until ∼3.5 hr after which a clear positive contribution to the mean kinetic energy is visible. After summing all the budget terms, the eddy transport contributes up to 15%-40% to the increase in mean KE of the vortex over the 6 hr period. When integrating the budget terms over time, the eddy transport contributes 18% to the mean kinetic energy with no notable differences between the models. At these later time periods, the vertical divergence of the vertical tangential momentum flux (the third term in the E contribution in Equation 16) has the largest values and provides a positive tendency to the mean kinetic energy in our analysis domain.
For the total explicit diffusion term, the values are very similar before ∼1 hr, but after that time the values from NUMA start to slowly increase relative to WRF with consequential differences at later times into the simulation (NUMA larger than WRF by ∼36% when averaged from 2 to 6 hr). The reason for this is that while the eddy viscosity values are fixed, the velocity gradients in NUMA are larger than WRF (described in the previous section), which produces a larger magnitude in the Laplacian operator. This is not a truly fair comparison of the IND in each code and a simple diagnostic calculation that accounts for this effect is outlined below. Figure 12 shows a time series of the summation of all terms on the right hand side of Equation 16 for WRF and NUMA (dashed lines). This figure represents the slope of the curves displayed in the time series of mean kinetic energy in Figure 7. Clearly, there is intensification throughout the simulation with the exception of the last output time at 6 hr, but note that this weakening is not shown in Figure 7 because the last output time is right at 6 hr. The intensification rate hovers between 0.03 and 0.05 m 2 /s 3 over most of the simulation with the exception of a large spike at 4.5 hr. After integrating these time series curves with the trapezoidal rule over 6 hr with the 30-min output interval, we find that NUMA has a larger mean (azimuthal mean and averaged over r = 10-50 km and z = 0.19-1.5 km) kinetic energy than WRF by ∼8%. However, this difference does not account for the larger explicit diffusion in NUMA mentioned above, which obscures the ability to isolate the effects of IND (the main goal of this study). To correct for this, we re-calculated the integrated mean kinetic energy with the explicit diffusion term (D) removed and we find that NUMA has larger values than WRF by ∼18%. We view the explicit diffusion term as a stabilizer of the numerical schemes and to some degree part of the sub-grid physics (turbulence), which is crude in this context despite the efforts to fit curves to the TKE output. Future work will examine the role of sub-grid physics in a much more detailed manner. Of course, the offline integration without explicit diffusion is not a perfect method for isolating the effects of IND. There is a coupled, nonlinear evolution of the fields whereby differences in the explicit diffusion affect the other terms in the mean kinetic energy during the simulations. This is difficult to control, and we do not address this issue here.
In summary, the main differences in the mean vortex intensity between WRF and NUMA is due to the radial pressure gradient contribution to the mean kinetic energy with smaller effects from the transport terms and explicit diffusion. However, even the transport and explicit diffusion terms are controlled, for the most part, by the pressure gradient term because the differences in each model's response to the input heating, via the pressure gradient, results in different velocities even very early (e.g., 1 hr) into the simulations (see Figure 11). Thus, differences in the calculation of the nonlinear advective terms in each model is not a significant source of diffusion and this result is consistent with Guimond et al. (2016). In addition, this result follows from the fact that both WRF and NUMA utilize high-order discretization of the advective terms and WRF showed no tangible differences when switching from the fifth order to sixth order stencil. Guimond et al. (2016) also identified the pressure gradient term as the controlling factor in dynamic core comparisons between three different models (including WRF and NUMA) for the same vortex analyzed here. However, Guimond et al. (2016) only considered idealized potential temperature perturbations to the initial state of the model as opposed to the time-dependent, 3-D observational heating used here. Guimond et al. (2016) conducted sensitivity tests with different order time integration schemes and found significant differences in the solutions, which led to the conclusion that the diffusion in WRF was due to the temporal discretization. Similar sensitivity tests were conducted in this work by comparing the control NUMA run (essentially a first-order in time method with the filter active) to a second-order in time Runge-Kutta method, which is very similar to WRF. These sensitivity tests in NUMA revealed small differences with peak absolute values of ∼1.5 m/s (not shown) indicating that the temporal discretization is not significantly affecting the solutions in either NUMA or WRF. For the spatial discretization of the pressure gradient term, WRF uses second-order finite differences while NUMA is utilizing fifth-order polynomials, which also provides fifth-order accuracy for the pressure gradient evaluation. Given the kinetic energy budgets and sensitivity tests described above, we assess that the low-order spatial approximation of the pressure gradient term in WRF is the source of the significant diffusion in the vortex intensity identified in this paper. The difference between Guimond et al. (2016) and the present paper is due to the nature of the forcing. In Guimond et al. (2016), the forcing was an initial condition (sharp in the time domain, smooth in space) while in the present study the forcing is from observations (smooth evolution in time, but sharp in space). This is the fundamental difference between these studies.
Given the large model differences in the pressure gradient contribution to the mean kinetic energy demonstrated in Figure 11, it is imperative to examine the full structure of this term to identify any potential localized signals. The components of the horizontal pressure gradient contribution to the kinetic energy in cylindrical coordinates are given as − and − where u and v are the radial and tangential wind speed, respectively. Figures 13 and 14 show horizontal cross sections of these terms, averaged over low-levels (∼0.19-1.5 km), at 4 and 5 hr into the simulations, respectively. In addition, select height-averaged (over the full column) heating inputs to the models leading up to these time periods are also shown in Figures 13 and 14. Figures 13a and 13b show the heating inputs at 3.33, 4 hr, which represent the heating snapshots leading up to the 4 hr mark in the simulations. The 4 hr heating has the larger weight in the model results, but there is still some "memory" of the heating from earlier times. The radial component of the pressure term in WRF ( Figure 13c) and NUMA ( Figure 13e) shows an azimuthal wavenumber-2 structure in the eyewall region, which is connected to the input heating structure most closely at 4 hr. The heating at 4 hr shows localized regions of large positive and negative heating rates (see, e.g., the feature to the West of the storm center in Figure 13b), which are correlated with the positive/negative couplet in the radial pressure term in a similar region. Note that due to the vortex drift to the South-East over time, the heating input and radial pressure term are not exactly aligned. Comparing the radial pressure term from WRF and NUMA shows that NUMA has larger values than WRF, especially in the localized positive regions. This is the reason why the azimuthal average of these fields (see Figure 11) shows NUMA with much larger values than WRF. This result indicates that strong, localized heating regions associated with convective bursts are producing a concomitant, enhanced pressure gradient response in NUMA that is driving the differences in the intensification of the vortex. The reduced, localized pressure gradient values in WRF are due to diffusion from the low-order spatial discretization of this term, as described above.
The azimuthal component of the pressure term for this time period in both WRF ( Figure 13d) and NUMA ( Figure 13f) shows a clear azimuthal wave structure with an average wavelength of ∼20 km and ∼wavenumber-five structure. These waves are very likely convectively forced vortex Rossby waves (Montgomery & Enagonio, 1998), but are not discussed in detail. The anticipated vortex Rossby waves in NUMA have a larger amplitude than those in WRF and this can be seen most clearly to the south of the vortex center. Note that small-scale oscillations in the azimuthal component of the pressure term are visible in both the WRF and NUMA fields. These oscillations are likely gravity waves as opposed to grid-scale noise in the models, but more analysis is needed to verify this hypothesis further. Figure 14 shows the same fields as Figure 13, only at 5 hr into the simulations. Similar results to those at 4 hr are observed, including larger magnitude, localized positive anomalies in the radial pressure term in NUMA (Figure 13e) compared to WRF (Figure 13c) that are connected to the heating snapshots at this time (Figures 14a  and 14b). The azimuthal pressure term also continues to show evidence of vortex Rossby waves with clearly larger amplitude features to the North and East of the vortex center. The heating input at 5 hr ( Figure 14b) shows that the majority of the heating is on the eastern side of the vortex with large, localized regions to the north-east of the center, which is consistent with the larger amplitude waves in NUMA.

Examining Longer Time Period Simulations
In the previous section, we studied the vortex response to 4-D observational heating applied over a short, 6 hr period in both WRF and NUMA. The 6 hr time period of these simulations was dictated by the specialized observations processed to compute latent heating rates. In this section, we extend these results to 18 hr simulations by introducing idealized heating for the first 12 hr and then adding in the 6 hr observational heating.  The idealized heating is constructed by summing potential temperature perturbations from several azimuthal wavenumbers (0, 2, 3, 4, and 5) with amplitudes of (2.0, 1.5, 1.0, 1.0 and 1.0 K), respectively. A wavenumber-one perturbation is not included because this asymmetry would cause a cross-vortex flow that would move the system out of the domain center. This combined potential temperature perturbation is converted to a heating rate by scaling the values by the maximum and multiplying by 100 to achieve a peak heating rate of 100 K/hr. The idealized heating is centered at 5 km height and 40 km radius, which is just inside the RMW of the initial vortex (∼50 km). Figure 15 shows a horizontal and vertical cross section of the idealized heating, which depicts an isolated convective cell in the Northern eyewall (Figure 15a) along with a ring-like structure in other regions of the eyewall with embedded asymmetries. This heating structure and its vertical characteristics (Figure 15b) represents a reasonable approximation to heating in the early stages of the RI process (Guimond et al., 2011;Park & Elsberry, 2013). The idealized heating is increased to maximum capacity in both models over a 30 min period and then held constant up to 12 hr. After that time, the 6 hr observational heating is applied in the same manner as the previous experiments to achieve an 18 hr simulation. Note that the control setup for both WRF and NUMA are utilized for this test. Figure 16 shows a times series of the maximum azimuthal mean windspeed in WRF and NUMA for the 18 hr simulation with output every 0.5 hr. After about 2 hr, the windspeed in NUMA becomes larger than WRF by ∼1.5 m/s, which is due to the implicit numerical diffusion effect identified in Section 3.2. However, between 8.0 and 14 hr, the windspeed in NUMA becomes significantly lower than WRF by up to ∼3.25 m/s at about 10 hr. What is going on in the models during this time period?
Figures 17a, 17d, and 17g show horizontal cross sections of windspeed perturbations for WRF, NUMA and the differences (NUMA-WRF), respectively, at 10.5 hr and 0.19 km height. The WRF windspeed shows very little azimuthal variability and a structure that is also relatively invariant in time when the fields are animated (not shown). The NUMA fields are in stark contrast at this time with significant azimuthal variability due to the mixing of momentum between the eye and eyewall by heating-induced mesovortices. The difference plot shows a strong signal of this momentum mixing as NUMA has higher windspeed values in the eye and weaker values in the eyewall, which produces a weaker maximum azimuthal mean windspeed as shown in Figure 16.
The idealized heating perturbation has a wavenumber-four and -five anomaly in the total field, which over time, generates a strong enough wavenumber-four and -five circulation to enable mixing between the eye and eyewall of the vortex. WRF does not produce this mixing phenomenon because the circulation asymmetries are damped out from the IND in the pressure gradient discretization. Thus, sometimes even qualitative differences can appear in the simulations in addition to significant quantitative differences. We believe that the mixing produced by NUMA is a more accurate reflection of the dynamics for this setup, but further analysis of this effect is left for future research. Returning to Figure 16, after 14 hr the differences between WRF and NUMA are small. Part of this is due to the fact that the mean intensity of the NUMA vortex has been reduced from the mixing process (recall that the observational heating is plugged into the simulations at 12 hr). In addition, the maximum azimuthal mean windspeed statistic is not always a good indicator of the total vortex intensity. Figure 17 shows the WRF and NUMA windspeed perturbations and associated difference fields for 15 hr (Figures 17b, 17e and 17h, respectively) and 17 hr (Figures 17c, 17f and 17i, respectively). The NUMA panels show much stronger windspeeds than WRF in the Northern and Western portions of the eyewall at both 15 and 17 hr with a large swath of ∼20 m/s differences at 15 hr and ∼10 m/s differences at 17 hr. Clearly, these large differences (under-predicted intensity from WRF) would be a problem for operational forecasters, but the intense regions of the eyewall in NUMA are not reflected in Figure 16. It is also worthwhile to note the moderate reductions in NUMA windspeed, relative to WRF, radially outward of the eyewall (∼−5 to −8 m/s) at both 15 and 17 hr, which is due to a more compact windfield in NUMA.  In addition to the control 18 hr simulations shown in Figure 16, simulations with WRF and NUMA that include noise in the initial conditions were also conducted. Random perturbations, within the range of ±0.5 m/s, in a ring centered on the RMW and below 2 km height were added to the initial zonal and meridional wind fields in each model. The motivation for these noise simulations is to represent uncertainty in the initial conditions that will stimulate nonlinear interactions and potential large scale random error at later times ("chaos"). Both the WRF and NUMA noise simulations show little to no differences in the maximum azimuthal mean windspeed intensity metric at all times. This indicates that at least for time periods before 13 hr, the implicit numerical diffusion signal is significantly larger than the noise. Beyond 13 hr, additional analysis is required to examine these effects given the substantial asymmetric windspeed differences demonstrated in Figure 17. Figure 18 shows a time series of the maximum azimuthal windspeed differences for the 18 hr simulations. To clarify, we calculated the azimuthal mean windspeed at each output time, a function of radius and height, for each model run (NUMA, WRF, NUMAnoise, WRFnoise) and subtracted the same quantity for a different member of this collection as denoted in the key of Figure 18. After the difference, the maximum value is reported. This procedure effectively characterizes the difference fields shown in Figures 17g-17i as a time series plot.
The control simulation (NUMA-WRF; green line) shows a general increase in the model differences with time up to peak values greater than 10 m/s at 14.5 hr, which reflects the larger overall vortex intensity response from NUMA at these later time periods where the observational heating is active. The control simulations with noise (NUMAnoise-WRFnoise; red line) shows how this numerical diffusion signal looks after considering chaos, revealing random perturbations to the control but no systematic changes, as expected. This noise-perturbed curve is one realization of an envelope of uncertainty that will exist about the control curve.
It is beyond the scope of the present study to conduct a large set of ensemble simulations to identify a more robust uncertainty estimate. However, we can identify the magnitude of the chaos signal present in each specific model by comparing differences of a noise-perturbed run with a noise-free run. The purple and blue lines in Figure 18 show this model specific chaos signal for NUMA and WRF, respectively. The magnitude of these curves are well below the control simulation (green line), especially for WRF in the time period beyond ∼12 hr, which is a critical period where the observational heating is applied. These results indicate that the reduction in overall vortex intensity identified in WRF from IND (the "signal") is larger than the chaos introduced by perturbing the initial conditions with small amplitude uncertainties (the "noise"). Furthermore, the magnitude of the simulation "noise" is larger in NUMA compared to WRF because of the extensive numerical dissipation in WRF.

Summary and Conclusions
In this paper, we have studied the computational fluid dynamics of the hurricane RI process by considering idealized simulations of the vortex response to time-dependent, 3D latent heating estimates derived from airborne radar measurements collected in the RI of Hurricane Guillermo (1997). Idealized heating perturbations that contain a summation of low-order (0, 2, 3, 4, 5) azimuthal wavenumbers were also analyzed. Two types of numerical models were considered: a community-based, finite difference and split-explicit model called WRF and an advanced, spectral element and semi-implicit model called NUMA. The models are carefully setup and analyzed to ensure the differences can be isolated to the numerical schemes that comprise the dynamic core. This includes explicit diffusion settings, which are parameterized based on output from a 3D TKE subgrid model experiment.
Prior studies used simple thermal perturbations to the initial conditions to represent the effects of convective heating and found that the WRF model had significant IND when compared to advanced research codes, including NUMA (Guimond et al., 2016). The current study also finds significant IND in WRF with a reduction in several intensity metrics over a 6 hr period: (a) maximum wind speeds in WRF are ∼12% lower than NUMA when matching the eddy diffusivity values, (b) time-integrated, mean kinetic energy values in WRF are ∼20% lower than NUMA when accounting for differences in the Laplacian diffusion operator and (c) peak, localized wind speed differ- ences in WRF are ∼12 m/s lower than NUMA. Reed and Jablonowski (2012) compared global simulations of hurricanes with several dynamic cores and found that the spectral element core produced the strongest storms, suggesting that this numerical method was the least diffusive. However, their results are clouded by the fact that seemingly different explicit diffusion mechanisms and coefficients are used in some dynamic cores, which does not allow isolation of the solution differences to the numerical methods.
Sensitivity studies show, that in order to achieve a similar intensity time series to NUMA, the explicit diffusion in WRF must be reduced drastically, with eddy viscosity values set to 50% of those in NUMA. In the control simulations, the NUMA windspeeds are visibly larger than WRF by roughly 5 m/s when averaged over the eyewall with local regions exceeding 10 m/s. In addition, NUMA's low wind region in the eye is slightly wider than WRF's, resulting in larger velocity gradients (and larger Laplacian diffusion) when accounting for the enhanced values in the eyewall.
Results from an 18 hr simulation, which incorporated both idealized and observational heating, continued to show large overall vortex intensity reductions in WRF with large swaths of the low-level eyewall 10-20 m/s lower than NUMA for several hours. Other intensity metrics, such as the maximum azimuthal mean windspeed, showed small differences between the models during the later portions of the simulations, which highlights the need to appropriately characterize the total vortex wind field. Of particular interest, during the middle portion of the 18 hr simulation, NUMA produced turbulent mixing between the eye and eyewall of the vortex from heating-induced mesovortices, which was not present in WRF due to the effects of IND.
To understand the nature of the differences between the models, the azimuthal mean kinetic energy budget was examined. At all time periods in the 6-hr simulation, the pressure gradient force contribution to the kinetic energy is significantly higher in NUMA compared to WRF by ∼23% in the mean and ∼40% in the maximum. Examination of the horizontal components of the pressure term reveal that NUMA produces localized pressure gradient anomalies that are larger in magnitude when compared to WRF. These localized regions are tied into the observational heating inputs that contain the presence of convective bursts, which are prevalent during the RI process. In addition, the presence of azimuthal waves in the pressure gradient term are visible in the simulations, likely vortex Rossby waves, with larger amplitudes in NUMA. While the axisymmetric transport of kinetic energy, especially by vertical fluxes, is substantially larger than the asymmetric transport, we find that these eddy processes contribute 15%-40% at 30-min output intervals over the 6 hr period and ∼18% when integrating the terms over time.
Sensitivity tests with different time integration schemes in NUMA were conducted to identify the root numerical causes of the model differences. However, employing a second-order in time scheme, compared to an essentially first-order in time method for the control run, did not produce any notable differences in the NUMA solutions. This is in contrast to the results in Guimond et al. (2016), where higher order time integration schemes produced even more energetic solutions. The reason for the discrepancy is likely due to the nature of the problem: Guimond et al. (2016) analyzed a freely evolving vortex initialized with a perturbation while the present study considered strong, 4-D forcing. Therefore, the significant numerical dissipation in WRF is controlled by a spatial discretization error and this is consistent with the fact that WRF relies on a diffusive, second-order approximation to the pressure gradient force while NUMA utilized a fifth-order accurate approximation.
A natural question to ask is: which model solution is correct? The best method for answering this question for the idealized simulations presented in this paper is to conduct a numerical convergence study. WRF simulations at 1 and 0.5 km horizontal grid spacing showed an increase in the vortex intensity (relative to the 2 km control), toward that of NUMA, through various metrics. These results demonstrate that some of the IND errors in WRF are being reduced or shifted to smaller scales by adding more grid points. However, a reduced vortex intensity (relative to NUMA) still remains in the 0.5 km WRF simulations. Considering that the overall accuracy of WRF is about second-order in space due to the pressure gradient discretization and NUMA is fifth-order in space, WRF would require at the bare minimum a factor of two increase in resolution in the x, y, and z dimensions (relative to the 2 km control) to approach the numerical error in the 2 km NUMA simulations. These estimates are based on the convergence curves for the 1-D advection equation for a second-order and fifth-order numerical method described in Giraldo (2020). Thus, further grid refinements, especially in the vertical (not analyzed), in WRF may bump up the intensity even closer to NUMA, which should be a more "correct" solution.
Excessive numerical dissipation is not a desired aspect of a modeling system because it reduces the effective resolution of the simulations and can damage the effects of physics-based sub-grid models and observations used to initialize the model in data assimilation practices. Returning to the major issue of the significant low bias in RI cases in operational models described in the introduction, this paper indicates that part of that problem could be due to IND. The chaos predictability horizon is always a looming factor in nonlinear dynamical systems, but this uncertainty should appear as a random error as opposed to the bias currently observed for RI cases.
Much work remains to examine these computational physics effects further. The simulations in this paper were restricted to dry dynamics and thus, the positive feedback loop between windspeed, surface fluxes of enthalpy, microphysical heating and pressure anomalies is shut off. We expect additional divergence in the models, with potentially further intensity reductions in WRF, when considering moist dynamics and other physical processes. Nevertheless, this paper makes an important step forward in an attempt to develop a holistic, thorough investigation of the computational fluid dynamics of the hurricane RI process.

Appendix A: Azimuthal Mean Kinetic Energy Equation Derivation
Beginning with the radial and tangential equations of motion on an f-plane and considering the anelastic approximation =̄( ) : Here u and v are the radial and tangential windspeeds, respectively, and all other variables are defined as in the main text. The total derivative in cylindrical coordinates is: Breaking Equations A1 and A2 into mean and turbulent parts, azimuthally averaging and transforming the resulting equations into flux form: for the radial equation, for example, we find:

Data Availability Statement
The WRF code is available at https://github.com/wrf-model/WRF. An overview of the WRF modeling system, along with information regarding downloads, user support, documentation, publications, and additional resources can be accessed through the WRF Model Users' Web Site: https://www2.mmm.ucar.edu/wrf/users/. Permission to access the NUMA source code can be requested from F. X. Giraldo (fxgirald@nps.edu). NUMA is a research code funded by ONR that shares similarities with the U.S. Navy's NEPTUNE model that is scheduled to go into operations at the end of 2023. We have stated available upon "request" because ONR had asked us to seek agreements with collaborators wishing to work with the code due to the connection of the base code to potentially sensitive Navy operations. Uploading the code to Zenodo would preclude those agreements from taking place. The WRF and NUMA output data are available at the following link: https://figshare.com/projects/The_Effects_ of_Numerical_Dissipation_on_Hurricane_Rapid_Intensification_with_Observational_Heating/126469.