Point by Point Response to the "White Paper" by Skamarock and Baldwin

Zavisa Janjic and Tom Black, NCEP

 

 

 

There are established mechanisms for independent evaluation of the WRF cores within the WRF project.  Thus, the effort volunteered by Skamarock and Baldwin was welcome but not necessary.

 

 

Interestingly, this one, and some of the subsequent claims by Skamarock and Baldwin were made on the basis of a single example of an NMM precipitation forecast.  We do not know the origin of this NMM forecast example.  However, in the upper panel of Fig. 1 we show an example in which contrary to the claims by Skamarock and Baldwin, the NMM precipitation forecast has more variance than the forecast by the Mass core (MC).  These examples demonstrate that it is not possible to make valid general conclusions on the basis of selected fragmentary evidence.

 

 

As will be explained later on, we disagree with this claim.

 

 

Note that all this analysis is done on a single example.  As already pointed out, we do not know the origin of this example.  The NMM was not yet operational at that time and has changed since then.

 

We have demonstrated in one of our previous responses that the rather smooth appearance of some of our precipitation forecasts is mostly a matter of tuning of precipitation physics, and not a result of numerical filtering.  The explanation of the alleged properties of the precipitation spectra by damping in the dynamics is unfounded.  Moreover, as can be seen from the random example from the Western Domain shown in Fig. 1, our precipitation forecasts are not necessarily smoother than the forecasts by the MC (upper two panels) as claimed by Skamarock and Baldwin.  On the contrary, the NMM forecast presented here shows more structure.  Needless to say, neither the precipitation example by Skamarock and Baldwin nor our example is sufficient for making general conclusions on the behavior of the dynamics of WRF cores as Skamarock and Baldwin do.

 

Use of the spectrum as a measure of success of precipitation forecasts is very unusual and its relevance is questionable.  The most important aspects of precipitation forecasts are ultimately how well models predict the location, timing, and amount of precipitation.  The success of precipitation forecasts concerning the location and the amount of predicted precipitation can be measured, respectively, e.g., by equitable threat score and by the bias score. 

 

The statistics on these two scores for the Western domain for November 2003 is shown in the lower panel of Fig. 1 for the Eta (red solid line), the GFS (dashed dark blue line), operational NMM (dashed green line), experimental WRF NMM (dashed light blue line) and MC (dashed purple line).  As can be seen from the plots, despite the fact that the MC generally shows less structure than the operational NMM and the experimental WRF NMM in this domain, the MC lags considerably behind other models in both the threat and the bias scores.  It is particularly disturbing that, as can be seen from the bias score, the MC tends to significantly over-predict the precipitation amounts.  This indicates that the model has problems with maintaining the hydrological balance, and more importantly, that the energy input by the phase changes of water is overestimated by a considerable amount.  This spurious energy input adversely affects both the model physics and dynamics.

 

Two other points can be made regarding Fig. 1.  The forecast from the MC predicted a large amount of widespread heavy precipitation over Missouri and northeast Oklahoma where virtually none actually fell.  It also greatly exaggerates the amount of rainfall in northwest Wisconsin while totally missing what did fall in east central Iowa.  The NMM and Eta forecasts did not produce these errors or did so to a far lesser degree.  The simple fact that the MC precipitation forecast generates structure certainly does not in itself translate into skill in forecasting the nature of convection, timing, location, or amount of rainfall.  Secondly we have demonstrated in one of our previous responses that the rather smooth appearance of some of the NMM precipitation forecasts is mostly a matter of tuning of precipitation physics and not a result of numerical filtering.  This is demonstrated again in Fig. 2 that was produced by modifying parameters within the BMJ convection scheme.  The result is that the detail and structure in the field have been greatly enhanced although some of the errors made by the MC forecast are also beginning to appear.  As development proceeds we want to provide as much detail to the forecasters as possible but not simply for the sake of detail itself.

 

 

The data that are supposed to support this claim are not available from the web page to which Skamarock and Baldwin refer.  We will come back to the question of dissipation later.

 

 

The claim that "direct measure of dissipation in the model's dynamics is its kinetic energy spectra" is untenable in the current context.  The kinetic energy spectrum certainly is not a direct measure of dissipation in the model's dynamics.  The spectrum produced by a mesoscale model in a short-range run in a limited area domain depends on the spectrum of initial data, on the synoptic situation, on the presence and characteristics of physical and spurious sources and sinks of energy, the size of the integration domain, the length of the integration, and, of course, on how well the model dynamics simulate the nonlinear energy cascade.

 

But let us first clarify what is the subject of the discussion.  Namely, Nastrom and Gage (1985) examined measurements made by commercial aircraft and found that one-dimensional kinetic energy spectra along their flight-paths in the lower stratosphere and in the upper troposphere follow the –5/3 slope in the range from several hundred kilometers to several kilometers.  As of now, there is no universally accepted explanation of this spectral shape.  Several possible explanations have been proposed (e.g. Gage, 1979; Lilly, 1983; Gage and Nastrom, 1986; Tung and Orlando, 2003).  They include the downscale nonlinear energy cascade and an inverse cascade from smaller to larger scales.

 

Using two-parameter quasi-geostrophic dynamics, Tung and Orlando (2003) demonstrated that given enough time to reach the statistical equilibrium, the –5/3 spectral range can be generated through downscale cascade of energy.  Similar statistical properties of the spectra were obtained in extended simulations using the GFDL SKYHI model with higher than usual yet still quite modest horizontal resolution for a climate model (e.g., Hamilton et al., 1999; Kosyik and Hamilton, 2001).

 

Concerning the properties of the motions in the –5/3 spectral range, as can be seen from the Nastrom and Gage (1985) spectrum, the spectral amplitudes of the waves with the largest energy are 103 to 106 times larger than the amplitudes in the –5/3 range.  So if the amplitudes in the large scale part of the spectrum correspond to a wind speed on the order of 10 m s-1, then the amplitudes in the –5/3 range correspond to wind speeds of a few cm s-1 to a few tens of cm s-1.  Apart from being representative for lower stratosphere and upper troposphere, this is obviously not enough to explain, e.g. an MCS and many other important mesoscale motions.  Thus, this part of the spectrum cannot be directly related to severe mesoscale phenomena as Skamarock and Baldwin repeatedly do.

 

Generally, it should come as no surprise if a mesoscale model develops the –5/3 spectrum considering that models with much simpler dynamics and coarser resolution have successfully done that (e.g., Hamilton et al., 1999; Kosyik and Hamilton 2001; Tung and Orlando, 2003).  However, the mesoscale runs typically are made on much shorter time scales than those considered in most previous studies.  Namely, the statistical properties of atmospheric spectra typically are investigated in extended integrations (tens or hundreds of days), and the spectra are averaged over long periods (tens or hundreds of days) in order to ensure that statistical equilibrium is reached.  The need for extended integrations and long averaging periods arise due to the time scale of the nonlinear cascade.

 

Despite the problem with the time scale, in short range runs the model spectrum perhaps still can be close to the statistical atmospheric spectrum if the initial data are well balanced and do not deviate too much from observed statistics.  However, consider the possibility that the initial data are not well balanced, and that they do deviate considerably from the observed atmospheric spectrum in the –5/3 range.  This could be due to the imperfections of the driving model or due to the set-up of the data assimilation system for example.  In addition, the size of the integration domain in mesoscale runs is typically smaller than the size of the large-scale atmospheric disturbances that feed the downscale nonlinear cascade.  Thus it appears that physical or spurious sources of energy other than the downscale nonlinear cascade from the large-scale motions are needed in order to generate and maintain the –5/3 spectra in mesoscale atmospheric models.

 

Possible sources that can generate and maintain the –5/3 spectrum in mesoscale models may be (a) physically justified mesoscale forcing, (b) early collapse of the spectrum due to spurious computational nonlinear cascade, (c) other small-scale computational errors such as the errors due to the representation of topography, etc.  In order to clarify the situation several tests have been performed using two versions of the NMM, namely the parallel version on the E grid (WRF-NMM) and the PC version on the B grid (NMM-B).  Both versions were designed applying the same principles in the discretization of the model dynamics and share the same physical package (Janjic et al., 2001; Janjic, 2003).

 

It should be noted that the WRF-NMM and the NMM-B are well qualified for investigating atmospheric spectra.  Their conserving of energy and enstrophy generally improves the accuracy of the nonlinear dynamics.  In particular, the energy and enstrophy conservation controls the nonlinear energy cascade and restricts an early spurious energy transfer toward smaller scales by nonlinear interactions.  The energy conservation improves the stability of the model and eliminates the need for excessive dissipation (either explicit or built into the finite-difference schemes) that could affect the spectra generated by the model.  In addition, the WRF-NMM and the NMM-B use a hybrid pressure-sigma vertical coordinate so that except for the errors propagating from below, in the upper troposphere and in the stratosphere they are relatively free of the sigma coordinate errors that are largest at higher altitudes.  Finally, explicit formulation of major dissipative processes allows precise “dosage” of dissipation.

 

The topography used in the tests is defined as grid-box means of the USGS 30’’ global data.  Except in ten rows along the lateral boundaries, no smoothing or filtering is applied.  The test with the operational set-up of the WRF-NMM over the Central Domain, but in the sigma coordinate, generated the –5/3 spectral slope on constant pressure surfaces in the upper troposphere and in the stratosphere.  An example of the spectrum (blue diamonds) obtained at the 300 hPa level, and averaged over forecast times from 24 to 36 hours with 3 hour intervals, is shown in Fig. 3.  As can be seen from the figure, the model develops the spectrum (blue diamonds) following the –3 (purple squares) and the –5/3 (yellow triangles) slopes that is in excellent agreement with observations.

 

However, there are problems with interpretation of this result.  Namely, when the operational set-up of the WRF-NMM is run in the hybrid coordinate, the spectrum at the small scales remained steeper than the –5/3 slope and approached the –3 slope (not shown).  Another reason for concern is that, as can be seen from Fig. 4, the spectrum of the square of unfiltered topography in the Central Domain (blue diamonds) follows the –5/3 law in the mesoscale range.  This result, together with the lack of the –5/3 range in the hybrid coordinate runs, challenges the interpretation of the result obtained in the sigma coordinate over mountainous areas as a legitimate Nastrom and Gage (1985) spectra.  Namely, the –5/3 spectrum in the upper troposphere and in the stratosphere could be generated by nonlinear cascade of spurious energy due to the sigma coordinate errors.  In this case, the computational noise would be mistaken for the Nastrom and Gage (1985) spectrum.  On the positive side one could argue that, although for a wrong reason, the nonlinear dynamics of the model still performed well by generating the –5/3 spectrum.  However, considering the shape of the mountain spectrum, the –5/3 spectrum observed in model integrations over mountainous areas could be simply a projection of the spectrum of topography and thus may have nothing to do with the nonlinear dynamics.

 

In order to investigate the possible effect of topography, the NMM-B was run using a resolution of 15 km in the horizontal and 32 levels in the vertical in a domain of the same size as the Central Domain, but over Atlantic Ocean.  In this way the possibility of mountains influencing the energy spectrum is eliminated.  In addition, as a major deviation from the operational set-up, the lateral diffusion was turned off in order to facilitate and hasten the accumulation of energy at the small scales.  However, weak divergence damping was still used.

 

The evolution of the NMM-B spectra at 300 hPa over Atlantic Ocean (blue lines) is shown in Fig. 5 at 6 hour intervals (6-24 hours, top to bottom in the left column and 30-48 hours, top to bottom in the right column).  The –3 (purple lines) and –5/3 (yellow lines) slopes are shown for comparison.  The time average over 36-48 hours of the spectra is shown in Fig. 6.  In the case considered, there was a small but vigorous extratropical cyclone in the northern part of the integration domain.  As can be seen from Fig. 6, by the end of the forecast, the model (blue diamonds) developed the –3 (purple squares) and the –5/3 (yellow triangles) spectral ranges that agree well with observations.  However, as can be seen from Fig. 5, it needed 24-36 hours to do so.  In the run with the same initial and boundary conditions, but with the physical package turned off, the –5/3 spectral range did not develop which indicates that in this case the physical processes provided the energy needed for spinning up the –5/3 spectrum.  It may be that the generation of the –5/3 spectrum was also aided by a more vigorous downscale nonlinear energy cascade since there was more energy than usual on the scales smaller than the size of the integration domain.

 

The WRF-NMM run in the Eastern Domain for the case of hurricane Isabel (initial data September 17, 18Z, from the Eta data) with the resolution of 8 km and 60 levels, and with no lateral diffusion, showed consistent behavior.  The evolution of the WRF-NMM spectra (blue lines) at 300 hPa in the Eastern Domain is shown in Fig. 7 at 6-hour intervals (6-24 hours, top to bottom in the left column, 30-48 hours, top to bottom in the right column).  The –3 (purple lines) and –5/3 (yellow lines) slopes are shown for comparison.  The time average over 36-48 hours of the spectra (blue diamonds) is shown in Fig. 8 together with the –3 (purple squares) and –5/3 (yellow triangles) slopes.  As can be seen from Fig. 8, by the end of the forecast the WRF-NMM also generated well developed –3 and–5/3 spectral ranges in agreement with the observations.  However, as can be seen from Fig. 7, the WRF-NMM needed some time to do so, but perhaps somewhat less than in the previous run of the NMM-B over the Atlantic Ocean.  Again, it may be that the generation of the –5/3 spectrum was also aided by a more vigorous downscale nonlinear energy cascade since there was more energy on the scales smaller than the size of the integration domain.

 

As expected on the basis of theoretical considerations, the presented results demonstrate that the nonlinear dynamics used in the NMM has been successful in reproducing the observed mesoscale atmospheric spectra, even at a rather modest resolution of 15 km, and even without the forcing by the mountains.  Whether the energy in the small-scale part of the spectrum comes from legitimate physical sources or from computational noise due to model imperfections is another issue that requires further investigation.  In other words, there is still no guarantee that the –5/3 spectrum of the model forecasts is generated by the same mechanisms as the Nastrom and Gage –5/3 spectrum observed in nature.  If it turns out that the model forecast spectra in short-range integrations are physically legitimate, and if the model simulates the nonlinear energy cascade well, the shape of the spectrum produced by the model could perhaps be used as guidance for more precise tuning of the model’s dissipation parameters.

 

The most important point here is that, contrary to the claims by Skamarock and Baldwin, the NMM dynamics apparently can and do spin-up the atmospheric spectrum virtually perfectly.  Interestingly, Skamarock and Baldwin have ignored these results.

 

 

 

It is not clear why Skamarock and Baldwin chose to show the results in their Figs. 4-6.  It is even more puzzling that they advertise these apparently problematic results as resounding success.  Namely, none of the spectra they show resembles the observed Nastrom and Gage (1985) spectrum, or the successfully simulated spectra by the NMM shown in our Figs. 3, 6 and 8.

 

For example, in our Fig. 9 we reproduced their Fig. 4, except that we added for convenience a few more lines in order to facilitate evaluation of the agreement of their spectrum with observed spectral slopes.  As can be seen from Fig. 9, the MC spectrum (red line) is significantly steeper than the –5/3 slope in the range where the –5/3 slope is observed (see the light blue line).  In addition, there is a sharp drop-off at the small-scale end of the spectrum that is not observed in nature.  But perhaps the most disturbing is the difficulty that the model apparently had in developing the –3 spectral range (see the dark blue line).  Namely, instead of following the –3 slope, and then gradually transitioning to the –5/3 slope, as in the case of the NMM spectra shown in Figs. 3, 6 and 8 and in the observations, their spectrum generally follows a relatively straight line (purple line in Fig. 9) that is both significantly shallower than –3 and steeper than –5/3.  This indicates that the MC has significant problems with computational nonlinear energy cascade.

 

It is also important that all spectra shown are computed over the central US.  We have pointed out to Skamarock and Baldwin that in this area the spectrum of topography squared follows the –5/3 slope (Fig. 4).  Thus, the deviation of the MC spectrum from the –3 slope may have been aided by the sigma coordinate errors (if the spectrum is indeed at the height of the Nastrom and Gage data), or even worse, generated in the postprocessing of the projection of topography on the sigma coordinate surfaces.

 

We also reproduced their Fig. 5 as our Fig. 10 except that we added for convenience a few more lines in order to facilitate comparison of their spectra with observed spectral slopes.  As can be seen from Fig. 10, the spectra again qualitatively follow straight lines throughout the entire spectral range (brown straight line on top of the red 22 km spectrum, dark green straight line on top of the green 10 km spectrum, and dark blue straight line on top of the blue 4 km spectrum) with slopes in between –5/3 and –3.  However, the slopes of the spectra become shallower with resolution. 

 

In comparison with the NMM, it should be noted that the MC cannot reproduce the –5/3 spectral range even with the 4 km resolution, while the NMM can do so with the 15 km resolution.  These results underscore the importance of energy and enstrophy conservation for proper simulation of atmospheric spectra.

 

Skamarock and Baldwin Fig. 6 does not add any new information.

 

 

Hamilton and Koshik showed climatological spectra, and not spectra spun up in short-range runs.  As already pointed out, the MC spectra are very problematic and therefore cannot serve as a standard for comparison.

 

 

Skamarock and Baldwin cannot know that the deviation from the –3 slope in their spectrum in the central domain is not due to the sigma coordinate errors unless they run their model using an efficient technique for controlling the sigma coordinate errors in the stratosphere and in the upper troposphere.  Such a technique is the application of the hybrid sigma-pressure coordinate used in the NMM.  The claimed similarity of the spectra over the ocean and over the land does not prove anything because the spectra in the two domains could have been spun up from different energy sources.  However, the Skamarock and Baldwin spectra over the land and over the sea are not quite similar.  It is encouraging to note that there is a slight deviation from the straight line in the MC ocean spectrum which supports the speculation about spurious forcing by topography in the runs over land.

 

 

 

 

Skamarock and Baldwin do not offer any scientifically valid evidence for their claim that “Neither of these arguments suggest that the spin-up time should be long.”  If this claim is true, why are the statistical spectra computed after spinning them up for hundreds of days and then averaging them over hundreds of days (e.g. Tung and Orlando 2003)?  The experimental evidence produced by the MC is not acceptable because Skamarock and Baldwin have not demonstrated that the MC can reproduce a realistic atmospheric spectrum (see our Figs. 9 and 10).

 

Concerning the eddy turnover time argument, take for example the scale of 1000 km.  From the Nastrom and Gage (1985) spectrum we find that the characteristic wind speed for the motions on this scale is of the order 1 m s-1.  Then, 1000 km/1 m s-1 = 1 000 000 s = 11.574 days.  For the scale of 100 km, the wind speed is of the order of 0.1 m s-1 so that the eddy turnover time is about the same as for 1000 km. 

 

Concerning the inverse cascade, there have been difficulties with simulating it in numerical experiments, and this hypothesis is still being debated (e.g. Tung and Orlando, 2003).

 

 

We have already pointed out that Skamarock and Baldwin erroneously associate the –5/3 spectral range with severe mesoscale weather phenomena.  Apparently there is not enough energy in the –5/3 range to explain such phenomena.  The bulk of the energy in the spectrum is too high up.

 

 

The “significant input of energy from external or physical mechanisms outside the dynamical equations” is precisely what was responsible for spinning up the –5/3 spectrum in our short-range runs, and that is why it is inappropriate to associate the –5/3 spectrum entirely with the dynamics in the present context.

 

 

Why would the frontal collapse affect the spectrum in the stratosphere?  Also, much of the forcing responsible for cyclogenesis belongs to the –3 range and not to the –5/3 range.

 

 

We also believe that mesoscale models should reproduce the observed –5/3 spectral slope.  Moreover, we have demonstrated that the NMM dynamics can and does produce spectra that agree perfectly well with the observations provided there is a sufficiently strong physical or spurious source of energy, and given enough time (see Figs 3, 6 and 8).  In contrast to that, and contrary to the claims by Skamarock and Baldwin, the evidence they presented does not show that the MC can spin up a spectrum that resembles the observed one at the resolutions at which the NMM can, or at any other resolution they showed (see Figs 9 and 10).  Before we can start a meaningful discussion on the present topic Skamarock and Baldwin need to demonstrate that the MC is capable of reproducing the atmospheric spectrum in mesoscale applications.

 

However, we disagree that the statistical properties of the atmosphere are necessarily reproduced in each realization in small limited area domains, and that mesoscale models should necessarily spin up the –5/3 spectral range in less than 6 hours if the slope in this wavenumber range is much steeper in the initial data.  We are unaware of the parts of turbulence theory and mesoscale dynamics that contradict our view.

 

As already pointed out, comparing the amplitudes of motions in the –5/3 spectral range with the amplitudes of the large scale motions, we arrive at the conclusion that the –5/3 spectral range is dynamically consistent with motions with wind speeds in the range from several cm s-1 to several tens of cm s-1.  Such motions are hardly visible with conventional methods for presenting meteorological information and certainly do not represent phenomena such as severe storms as Skamarock and Baldwin appear to believe.

 

 

The other two primary filtering mechanisms in the MC used at NCEP following NCAR recommendation are:

 

smdiv (default value is 0.)

         This parameter specifies the divergence damping (0.1 is typical – and used at NCEP).

         Used in RK schemes with time-splitting to selectively damp

         acoustic noise.

 

emdiv (default value is 0.)

         This parameter specifies the external mode damping (whatever that may be) in mass model.

         Use it only in real-data cases in mass model (0.01 is typical – and used at NCEP).

 

 

Our intention has been to represent the dissipation in the atmosphere in a physically justified way.  The physical models of dissipation we have chosen, and particularly the tuning parameters that control the intensity of dissipation, are open for discussion and will be addressed in due course in future parallel tests since the NMM has been frozen for a considerable time as the transition to WRF has taken place.

 

 

Takemi and Rotunno (2003) seemed to be largely unsuccessful in reproducing the –5/3  spectral slope for w2 with the MC in their 3D turbulence test.  In contrast to that, the NMM-B produced the –5/3 range in a similar experiment without difficulties.

 

 

The key issue in our misunderstandings appears to be the belief of Skamarock and Baldwin that we are removing the noise from our forecasts by heavy numerical filtering.  That impression is incorrect.  While conserving energy almost exactly in the spatial finite-differencing (except for weak built-in divergence damping), our model dynamics do not generate small scale noise because it was designed not to do so following physical principles that govern the nonlinear dynamics of the real atmosphere (energy and enstrophy conservation in case of rotational flow).  Namely, following Arakawa, our modeling principle has been to prevent generation of computational noise rather than to employ continuous removal of the noise by numerical filters after it is generated.  Thus, it is our contention that our sophisticated numerical schemes based on physical considerations, and not the numerical filtering, are the main reason why our forecasts are to a large extent noise free.  At the same time, despite the fact that some forecasts have a smoother general appearance, we have presented examples demonstrating that our model dynamics produces spectra that agree extremely well with observations in both the –3 and the –5/3 ranges, even with 15 km resolution.  Actually, improved accuracy of the nonlinear energy cascade is probably one of the reasons the NMM produces such spectra.

 

The effect of dissipation in the operational NMM can be clearly seen in comparison with the WRF-NMM forecasts shown at http://wwwt.emc.ncep.noaa.gov/mmb/mmbpll/nestpage/.  Namely, the WRF-NMM is run with very little lateral diffusion.  Still, the forecasts by the operational NMM and by the WRF-NMM are similar in appearance.  The differences in the forecasts between the two models are mainly due to the different specifications of initial and boundary conditions and to subtle differences in their respective physical packages.  In contrast to the two NMM’s, the MC forecast that is also available at the site shows a much higher level of computational noise.

 

Extrapolating experiences gained with the MC to our models is simply not justified.  NCAR has chosen a discretization approach based on higher-order formal accuracy and it is not surprising that the experiences of Skamarock and Baldwin are different.  Namely, this approach is not without problems.  Experience with fitting higher-order polynomials to noisy data indicates that the results of discretization based on higher order formal accuracy may be locally very inaccurate in case of large amplitudes of small scale motions such as those forced by the sigma coordinate errors over rugged topography, or by convection.  In addition, higher order finite-differences require additional computational boundary conditions which also generate computational noise.  Finally, higher order formal local accuracy does not give any guarantee about the accuracy of energy transport by nonlinear interactions and the only means for controlling the nonlinear cascade in the MC is computational filtering.

 

 

The amplification factors sent to Skamarock and Baldwin say otherwise (see Fig. 11).  So, the dissipation in the MC is not smaller than in the Eta and in the NMM, and in contrast to the dissipative MC advection schemes, the dissipation in the Eta and the NMM is not hardwired.  Therefore the dissipation in the Eta and the NMM can be easily reduced further while the dissipation in the MC can be reduced only by increasing resolution because of the dissipative advection schemes.

 

 

The spectrum is not evidence of damping.  Again, the spectrum produced by a mesoscale model in a short-range run in a limited area domain depends on the spectrum of initial data, on the synoptic situation, on physical and spurious sources and sinks of energy, the size of the integration domain, the length of the integration, and, of course, on how well the model dynamics simulate the nonlinear energy cascade.  And, again, the small-scale noise may not be there because it is not generated by the inherent flaws in numerical algorithms.

 

 

It is imperative to note that second order schemes are not all equal.  In particular, we have been using a compact energy and enstrophy conserving scheme (Janjic, 1984) which controls spurious nonlinear downscale cascade.  This scheme performs differently from the most straightforward second order scheme used in these tests and does not even reduce to such a scheme when linearized.  Since we are dealing with spectra, we should expect that the energy and enstrophy conservation matters.  And as can be seen from Skamarock and Baldwin Fig. 11, it does.  As in previous examples, the MC spectrum is somewhere in between the –3 and –5/3 in the entire spectral range except at the end where it falls off sharply.  This should be a reason for concern since even the large scale part of the MC spectrum is too shallow presumably due to a spurious nonlinear energy cascade.  It is well known that the dissipation is the primary means for controlling the spectrum in model formulations that disregard energy and enstrophy conservation.  The dissipation appears to be helpful in the case of MC since it at least improves the shape of the spectrum at the large scales where it should approach the –3 slope.

 

From the description of the experiments we understand that the tests with the second order centered advection were done as follows:

 

a.  “we use 2nd order centered advection (both vertical and horizontal, to remove the upwind filtering from the standard mass-core WRF configuration), we have turned off all other computational spatial filters,” (presumably the red line in the plot)

 

b.  as in a. but with “we use the lower-bound viscosity from NMM in a second order horizontal damping term” (presumably the green line),

 

c.  as in b. but with “we use horizontal divergence damping with the same coefficient used in NMM” (presumably the blue line).

 

Thus, in this experiment series only the sensitivity of the second order MC spectra to dissipation alone was examined.  This experiment bears no relevance to dissipation properties of the NMM.

 

 

Our point is that the NMM forecasts have generally smooth appearance primarily because its sophisticated nonlinear schemes produce little computational noise.

 

 

We view lateral diffusion as a legitimate physical process, not as a numerical filter, and thus model it accordingly.  It is well known that there is no justification based on physical considerations for using diffusion operators higher than second order. 

 

It is also important to note that our lateral diffusion is not a second order filter because it is not linear.  Whether the sixth order filter is more selective than the nonlinear lateral diffusion is not obvious, and in contrast to using nonlinear lateral diffusion, there is no physical justification for using hardwired numerical filters.

 

Besides, the lateral diffusion is also used in the MC in addition to numerical filters.  At NCEP we had to double the recommended MC diffusion coefficient in order to keep the integrations of the MC stable with our physical package.

 

The reason why we don’t use high order numerical filters is not because we haven’t heard of them.  The reason is that we believe that we are doing better without using them at all.

 

 

Historically, the lateral diffusion coefficient in the Eta was being adjusted to damp the two-grid-interval wave in a given time.  This has been a common practice in NWP models (see e.g. how it used to be done in the UK MetOffice models at http://www.cgam.nerc.ac.uk/um/doc/umug/html/node110.htm).  This practice implies that the lateral diffusion is actually little more than a numerical filter, and for this reason we consider it inappropriate.  However, the lateral diffusion coefficients obtained in this way were not excessive.

 

The NMM dissipation was tentatively set in such a way as to mimic the total dissipation in the Eta pending a more detailed study when parallel runs become available.  We do not claim that the dissipation in the current operational NMM may not be overestimated, but we want to approach this issue in a more rigorous way.

 

 

Skamarock and Baldwin again erroneously consider the –5/3 spectrum synonymous with significant mesoscale phenomena.  As we have indicated, there is not enough energy in the –5/3 range to explain such phenomena as MCS’s.

 

 

This is not true.  As can be seen from Fig. 3, consistent with theory, the NMM as operationally set up can reproduce the –5/3 spectrum provided a sufficiently strong physical or spurious energy source is present, and given enough time.  Furthermore, the example shown in Fig. 1 demonstrates that the operational setup of the NMM can produce more meso-scale structure than the MC, and, moreover, that this additional structure verifies as can be seen from NMM’s significantly better threat and bias scores in the Western Domain.

 

In contrast, the evidence presented indicates that the MC has a serious problem in reproducing the observed atmospheric spectrum.  A typical MC spectrum approximately follows a single straight line throughout the spectral domain, except at the short-wave end of the spectrum where it sharply drops off.  The slope of this straight line lies in between –3 and –5/3 and decreases with increasing resolution.

 

 

This is perhaps the most fantastic conclusion that has been reached by Skamarock and Baldwin in view of the numbers themselves and in attempting to consider perceived smoothness.  With the same spatial resolution, the NMM is about three times faster than the MC.  Considering that the NMM code can be further optimized, this factor may be increased to about four.  On the other hand, as can be seen from the sharp drop-off of the MC spectra, its nominal resolution is reduced by a factor of about two by its filters.  Smoothness is completely irrelevant to efficiency since the amount of structure in the NMM precipitation fields can be increased easily through the existing physics packages.

 

In other words, in the quasi-operational runs at NCEP with the same resolution and using NCEP’s standard skill scores, the MC produced forecasts that are inferior in important aspects but at three times the computational cost of the NMM forecasts.

 

 

 

Nobody at NCEP is advocating “running high resolution forecasts using numerical filters that systematically remove the resolution gained by the refined grid.” Our experience at NCEP has been that the forecasts have been systematically improving with increasing resolution, which would not have happened if we have been “systematically removing the resolution gained by the refined grid.” 

 

On the other hand, as a matter of principle, we do not see justification for producing and retaining in the forecasts computational noise that is difficult to distinguish from legitimate mesoscale processes.  This statement should not be interpreted as the advocating of numerical filtering.  Namely, our modeling principle has been to prevent generation of computational noise rather than to continually remove the noise by numerical filters after it is generated.

 

We have been trying to represent the dissipation in the atmosphere in a physically justified way.  The physical models of dissipation we have chosen, and particularly the tuning parameters that control the intensity of dissipation, are open for discussion and will be addressed in due course in future parallel tests since the NMM has been frozen for a long time.

 

 

Skamarock and Baldwin advocate changing the verification rules in such a way that the verification versus observed data be reduced in importance in favor of comparison with statistical properties such as the Nastrom and Gage (1985) spectrum despite the fact that comparison with observations remains one of the most fundamental means of verifying the forecasts.  Their approach is inapplicable in day to day verification, since it is inappropriate to use statistical properties of the atmosphere to verify a particular realization, as well as potentially misleading, since the –5/3 range of the spectrum can be spun-up by computational errors.  Thus, a model may be rewarded for computational noise it makes and penalized for forecasts that agree with observations.

 

Comparison of conventional scores for the forecasts run at NCEP shows that the MC scores are consistently lower than the scores of the Eta and the NMM.

 

The experience at NCEP indicates that the conventional skill scores are still being improved with increasing resolution at horizontal resolutions of 10 km or so.

 

We look forward to the use of new methods of measuring skill in high resolution forecasts.  Until such methods are available and widely accepted, the current ones continue to provide a valuable tool in gauging forecast skill.

 

 

 

The MC appeared to be largely unsuccessful in reproducing the –5/3 spectral slope in the Takemi and Rotunno (2003) tests on cloud scales.  In contrast, the NMM had no difficulties with the –5/3 spectral range in a similar test.

 

 

Generally, this is not true.  The NMM dynamics has been tested on the cloud resolving scales in many runs performed primarily with the NMM-B.  The model passed successfully all the tests to which it has been subjected.

 

            Still, this claim by Skamarock and Baldwin is perhaps partly understandable considering that we have not published much about our results.  Due to resource constraints, the efforts at NCEP have been focused on the resolutions that will be most relevant for NCEP’s NWP applications in the near future.  Nevertheless, we are confident that the NMM will perform well at the cloud-resolving scale should researchers want to test this model application.

 

 

This paragraph does not make sense.  No one would run a model with the same dissipation parameters at 10 km resolution and at 1 km resolution.  At 1 km resolution even the parameterization of dissipation should be different (3-dimensional).

 

 

Although our very limited resources have not allowed us to publish much, particularly about the results on the cloud resolving scales, some relevant results have been presented to the scientific community.  For example, the results of explicit simulations of convection using different microphysics schemes were presented in the discussion at the Convection Workshop held at NASA in December 2001.  These runs were made with the horizontal resolution of 1 km, and 32 layers in the vertical.  The –5/3 spectrum spun up in a test analogous to that of Takemi and Rotunno (2003) (that the MC had difficulties with) was shown at the SRNWP meeting on nonhydrostatic modeling in Bad Orb this year.  Finally, the NMM is running daily with 2 km and 4 km resolutions in Switzerland, so that presumably there is more experience with the NMM forecasts run at 2 km and 4 km resolutions than with those produced by the MC.

 

 

Why before?

 

 

The same applies to the MC considering its notable tendency to overpredict hurricanes.

 

 

As demonstrated many times in this text, the claim that the NMM WRF is heavily damped, and in particular, that it is more damped than the MC WRF, is arbitrary and unfounded.  Even if there is too much damping in the current operational set-up of the NMM WRF, this damping is not hardwired as in the MC WRF finite-difference algorithm and can be easily reduced if needed.  In contrast, the dissipation in the MC due to its high order filters can be reduced only by increasing resolution.  The simple enhancement of detail in such fields as precipitation in the NMM can be done easily through the current physics packages.

 

REFERENCES

 

Gage, K.S., 1979:  Evidence far a k−5/3 Law Inertial Range in Mesoscale Two-Dimensional Turbulence. J. Atmos. Sci.  36, 1950–1954.

Gage, K.S., Nastrom, G.D., 1986: Theoretical Interpretation of Atmospheric Wavenumber Spectra of Wind and Temperature Observed by Commercial Aircraft During GASP. J. Atmos. Sci, 43, 729–740.

Hamilton, Kevin, Wilson, R. John, Hemler, Richard S., 1999: Middle Atmosphere Simulated with High Vertical and Horizontal Resolution Versions of a GCM: Improvements in the Cold Pole Bias and Generation of a QBO-like Oscillation in the Tropics.  J. Atmos. Sci, 56, 3829–3846.

Janjic, Z. I., 2003:  A Nonhydrostatic Model Based on a New Approach.  Meteorol. Atmos. Phys., 82, 271-285.

Janjic, Z. I., J. P. Gerrity, Jr. and S. Nickovic, 2001: An Alternative Approach to Nonhydrostatic Modeling.  Mon. Wea. Rev., 129, 1164-1178.

Koshyk, John N., Hamilton, Kevin, 2001: The Horizontal Kinetic Energy Spectrum and Spectral Budget Simulated by a High-Resolution Troposphere–Stratosphere–Mesosphere GCM. J. Atmos. Sci, 58, 329–348.

Lilly, D.K., 1983: Stratified Turbulence and the Mesoscale Variability of the Atmosphere. J. Atmos. Sci, 40, 749–761.

Nastrom, G.D., Gage, K.S., 1985:  A Climatology of Atmospheric Wavenumber Spectra of Wind and Temperature Observed by Commercial Aircraft. J. Atmos. Sci, 42, 950–960.

Takemi, T., Rotunno, R., 2003:  The Effects of Subgrid Model Mixing and Numerical Filtering in Simulations of Mesoscale Cloud Systems.  Mon. Wea. Rev, 131, 2085-2101.

Tung, Ka Kit, Orlando, Wendell Welch, 2003: The k−3 and k−5/3 Energy Spectrum of Atmospheric Turbulence: Quasigeostrophic Two-Level Model Simulation.  J. Atmos. Sci, 60, 824–835.

 

 

 

 

 

Fig. 1.  Examples of the NMM (upper lef panel) and the MC (upper right panel) forecasts and threat (upper part of lower panel) and bias (lower part of lower panel) scores for November 2003 for the Eta (red solid line), the GFS (dashed dark blue line), operational NMM (dashed green line), experimental WRF NMM (dashed light blue line) and MC (dashed purple line).

 

Fig. 2

 

 

Fig. 3.  The WRF-NMM spectrum (blue diamonds) at 300 hPa averaged over forecast times from 24 to 36 hours with 3 hour intervals produced with operational set-up of the NMM in Central Domain, but in the sigma coordinate.  The –3 (purple squares) and –5/3 (yellow triangles) slopes are shown for comparison.

 

Fig. 4.  The spectrum (blue diamonds) of the square of unfiltered topography height in the NMM-B Central Domain with 15 km resolution.  The –3 (purple squares) and –5/3 (yellow triangles) slopes are shown for comparison.

 

 

Fig. 5.  Time evolution of the NMM-B spectra (blue lines) at 300 hPa over Atlantic Ocean.  .  The –3 (purple lines) and –5/3 (yellow lines) slopes are shown for comparison Run starting from 12 UTC, 09/07/2003, GFS data, 15 km, 32 levels resolution.  No lateral diffusion, weak mass divergence damping.  Plots every 6 hours, top to bottom, 6-24 left column, 30-48 right column.

 

 

Fig. 6.  Time average over 36-48 hours of the NMM-B spectra (blue diamonds) at 300 hPa over Atlantic Ocean.  The –3 (purple squares) and –5/3 (yellow triangles) slopes are shown for comparison.  Run starting from 12 UTC, 09/07/2003, GFS data, 15 km, 32 levels resolution.  No lateral diffusion, weak mass divergence damping.

Fig. 7.  Time evolution of the WRF-NMM spectra (blue lines) at 300 hPa in the Eastern Domain.  Run starting from 18 UTC, 09/17/2003 (Isabel), Eta data, 8 km, 60 levels resolution.  The –3 (purple lines) and –5/3 (yellow lines) slopes are shown for comparison.  No lateral diffusion, weak mass divergence damping.  Plots every 6 hours, top to bottom, 6-24 left column, 30-48 right column.

 

Fig. 8.  Time average over 36-48 hours of the WRF-NMM spectra (blue diamonds) at 300 hPa in the Eastern Domain.  The –3 (purple squares) and –5/3 (yellow triangles) slopes are shown for comparison.  Run starting from 18 UTC, 09/17/2003 (Isabel), Eta data, 8 km, 60 levels resolution.  No lateral diffusion, weak mass divergence damping.

 

 

 

Fig. 9.  Skamarock and Baldwin Fig. 4, except that light blue (-5/3 slope), dark blue (-3 slope) and purple (fit to the Skamarock and Baldwin spectrum) lines were added in order to facilitate comparison of the Skamarock and Baldwin spectrum (red line) with observed spectral slopes.  Instead of following the –3 slope, and then gradually transitioning to the –5/3 slope, the Skamarock and Baldwin spectrum generally follows a straight line (purple) that is both significantly shallower than –3 and steeper than –5/3.

 

 

Fig.  10. Skamarock and Baldwin Fig. 5 except that brown (fit to their 22 km spectrum), dark green (fit to their 10 km spectrum) and dark blue (fit to their 4 km spectrum) straight lines were added for convenience.  The spectra qualitatively follow straight lines throughout the entire spectral range with slopes in between –5/3 and –3, and generally become shallower with resolution. 

 

 

Fig. 11.  Amplification factors of the Mass Core RK vertical advection scheme (upper left panel) and the Mass Core horizontal advection scheme (lower left panel).  The Courant number for flow along the coordinate axes is 0.5, and it is 0.707 along the diagonal.  For convenience, the missing contour labels are added in red and marked by z.j. with the aid of another plot from the same source (a presentation by Wicker).