Preprint, International Workshop on Limited-area and Variable
Resolution Models. Beijing, China, October, 1995.

The Advanced Regional Prediction System and Real-time Storm-scale Weather Prediction
Ming Xue1, Kelvin K. Droegemeier1,2 and Vince Wong1
1 Center for Analysis and Prediction of Storms
2 School of Meteorology
University of Oklahoma
Norman OK 73019, USA
E-mail: mxue@uoknor.edu

1. INTRODUCTION

A three dimensional nonhydrostatic model, the Advanced Regional Prediction System (ARPS) has been developed at the Center for Analysis and Prediction of Storms (CAPS) for the past several years (Droegemeier, et.al., 1995; Xue, et.al., 1995; CAPS, 1992). The primary goal of CAPS is to demonstrate the feasibility of numerical weather prediction (NWP) on the storm scale. To achieve this, a numerical model is required that meets a number of criteria: the model must be able to accomodate, through various assimilation strategies, new data of higher temporal and spatial density than what have been traditionally used. An example is the WSR-88D Doppler radar data. The model must also serve as an effective tool for studying the dynamics and predictibities of storm scale weather. Finally, the model should take full advantage of parallel computing power so that operational prediction can be done in a timely manner in the near future. Based on these needs, ARPS has been developed. This paper first describes the model dynamics and physics of the current ARPS, then presents two examples of real-time real-data forecast.

2. MODEL DESCRIPTION

2.1. Model Equations

The governing equations of ARPS include the prognostic equations for Cartesian velocity components (u, v, w), perturbation potential temperature ( ), perturbation pressure (p'), mixing ratios for water vapor (qv), cloud water (qc), rainwater (qr), cloud ice (qi), snow (qs) and hail (qh), and the subgrid scale turbulent kinetic energy (E). The equation of state for moist air is also used. These equations are represented in a curvilinear coordinate system ( ) that has a coordinate surface following the terrain at the model lower boundary. In addition to the above equations, ARPS also solves prognostic equations for the surface and deep soil temperature, surface and deep soil moisture, and the canopy moisture in a coupled two-layer soil model.

In ARPS, a special curvilinear coordinate being orthogonal in the horizontal directions is used. The transformation is defined as

(1)

where x, y and z are the independent variables in the Cartesian coordinate. As shown in Sharman et al. (1988), the contravariant velocities can be expressed as the functions of Cartesian velocities u, v and w, which for our special case are:
, , (2)

Here, the transformation Jacobians and are defined as
, , and . (3)

The transformation relations for spatial derivatives from Cartesian coordinate (x,y,z) to the curvilinear coordinate ( ) are
(4)

For the convenience of notation, we define
(5)

where is the air density. In the curvilinear coordinate, the gorverning equations of ARPS for Cartesian momentum components , perturbation pressure , perturbation potential temperature ( ) mixing ratios (q) of water and ice species, and subgrid scale turbulent kinetic energy (E) become:

(6)

(7)

(8)

(9)

(10)

(11)

(12)

where advection operator ADV( ) is defined as

In the above equations, over-barred variables represent base-state atmosphere, and primed variables are deviations from the base state. In ARPS, the base state is defined as a function of the Cartesian coordinate z only, it is, however, dependent on all three indepdendent variables in the transformed coordinate. The base state satisfies the hydrostatic relation. In Eq.(9), , where is the total liquid and ice water.

Terms on the l.h.s. of horizontal momentum equations (6) and (7) are the local time tendency and perturbation pressure gradient force, and on the r.h.s. are, respectively, advection, Coriolis force, and turbulent and computational mixing terms. Here, all components of Coriolis force are included (f is the angular velocity of the earth and is the latitude). The contribution from vertical velocity can be as large as the horizontal contribution for convective scale motion. Since the base-state pressure is assumed to be horizontally homogeneous, its horizontal gradient terms do not appear in the horizontal momentum equations. The absence of these terms reduces the computational error associated with terrain-following coordinates, as discussed by, e.g., Janjic (1977).

The vertical momentum equation (8) has, on the l.h.s., the local time tendency, pressure gradient force, and the buoyancy due to compression, and on the r.h.s., the advection, thermal and water loading buoyancy, Coriolis force and mixing terms. The total buoyancy, defined by , can be derived from the equation of state for moist air
(13)

where g is the acceleration due to gravity, is the ratio of the gas constants for dry air ( ) and water vapor ( ). T is the air temperature. in Eq.(8) is the acoustic wave speed, with being the ratio of the specific heat of air at constant pressure ( ) and constant volume ( ).

The terms involving in Eqs. (6) - (8) are artificial "divergence damping" terms designed to attenuate acoustic waves, where is the density weighted divergence and the damping coefficient. By performing a divergence operation on the momentum equations, one can see that the "divergence damping" terms form a diffusion term in the 3-D divergence equation, thereby damp acoustic wave modes. Skamarock and Klemp (1983) showed that unstable acoustic modes due to mode-splitting time integration scheme can be effectively controlled by divergence damping. Since atmospheric flows are quasi-incompressible, the divergence damping has little adverse effect on meteorologically significant waves.

Pressure equation (9) is obtained by taking the material derivative of the equation of state and replacing the time derivative of density by velocity divergence using the continuity equation. The terms in the pressure equation are, on the l.h.s., the base-state pressure advection and velocity divergence term, and on the r.h.s., the perturbation pressure advection and diabatic heating and moisture/water effect terms. The hydrostatic relation was used to substitute for the vertical gradient of in the advection term. In general, the divergence term is dominant for meteorological applications, and the diabatic and moisture/water effect terms are small and are, as is done in Klemp and Wilhelmson (1978), neglected in the ARPS.

It should be noted that in obtaining the above equations, linearization approximations are made. The state variables that appear in the coefficients of certain terms are replaced by their base-state values. This is true for the density appearing in the pressure gradient force. These approximations are consistent with those made in anelastic systems of Ogura and Phillips (1962).

Equation (10) is a conservation equation for potential temperature. The terms in the equations are, from left to right, the local time tendency, vertical advection of ,advection, subgrid scale and computational mixing and diabatic source or sink terms. Moist and microphysical processes are included in the source/sink term and surface heat flux is handled by the subgrid scale mixing. The horizontal advection of vanishes because is horizontally homogeneous. The absense of this term reduces the noises associated with the calculation of horizontal gradients in transformed coordinates, similar to the treatment of gradient terms.

The conservation equations for the mixing ratios of water vapor qv, cloud water qc, rainwater qr, cloud ice qi, snow qs and hail qh are represented by Eq.(11) for a generic variable q. The right hand side terms are, in order, advection, sedimentation, mixing and source terms. The source term Sq includes all microphysical contributions. The sedimentation term represents the falling of hydrometeors (qr, qs and qh) at their respective terminal speed (Vq). Cloud water and cloud ice are usually assumed to float with the air, therefore their flow-relative terminal velocity is zero.

The turbulent mixing process represented by D in Eqs.(6)-(11) has to be parameterized using a closure scheme. In ARPS, three closure options are available. They are, respectively, the first order Smagorinsky/Lilly scheme (Smagorinsky, 1963; Lilly, 1962), the 1.5 order TKE-based scheme (Deardorff, 1975; Klemp and Wihelmson, 1978; Moeng, 1984) and Germano dynamic closure scheme (Germano et al., 1991; Wong, 1992; Wong and Lilly, 1994; Wong, 1994).

The 1.5 order scheme requires the solution of turbulent kinetic energy (TKE) equation (12). In Eq.(12), the terms on the right hand side are, respectively, the TKE advection, potential-kinetic energy conversion, shear production, dissipation and diffusion of turbulent kinetic energy. The potential-kinetic energy conversion term C can be found in Xue et. al. (1995). The coefficient in the dissipation term is related to the length scale l and model grid scale and is adjusted at the lowest grid level after Moeng (1984) and Deardorff (1980).

For the momentum equations, D is expressed in terms of Reynolds stress tensor , which is in trun related to deformation tensor :
(14)

where index represents the Cartesian coordinate directions and the repeated lower indices denote summation. is the turbulent mixing coefficient for momentum in j direction.

The turbulent mixing for potential temperature and any water substance can be written in terms of their turbulent fluxes , which is further related to their spatial gradients:
, (15)

where is the turbulent mixing coefficient for in j direction. In general, is taken to be the same for heat, moisture or water/ice quantities and where Pr is the turbulent Prandtl number.

In Smagorinsky scheme (Smagorinsky, 1963; Lilly, 1962), is given by
, (16)

where k = 0.21 after Deardorff (1972) and |Def| is the magnitude of deformation. is the a measure of the grid scale. On a model grid with similar grid spacing in all directions, the turbulence is nearly isotropic, for all j. When the grid aspect ratio is large, determined according to (20) and (21) can become too large for the vertical mixing, therefore we use We call this a case of anisotropic turbulence.

In the 1.5 order TKE scheme, the mixing coefficient is related to E and length scale l :
(17)

For isotropic turbulence, l is equal to in unstable or neutral case, but is further limited by length scale that is related to E and static stability N in stable case . The turbulent Prandtl number is given by . For anisotropic turbulence, 2 , and .

In the Germano scheme, the mixing coefficients are also given by (20), but parameter k is spatially and temporally varying and is dynamically determined based on local flow. The detailed implementation can be found in Xue, et.al. (1995).

Apart from the turbulent mixing, terms D in Eqs.(6)-(11) also include optional second order and fourth order computational mixing terms. These terms acts to damp computational noises of small scales.

2.2. Other Model Physics

Parameterization of cloud microphysics processes are important for storm scale prediction. In ARPS, the microphysics include the parameterized Kessler-type two-category liquid water scheme and the three-category ice scheme after Lin et.al. (1983). The implementation in ARPS follows that of Tao and Simpson (1993) and includes the ice-water saturation adjustment procedure of Tao et.al. (1989). A four category ice scheme with the ice particle number concentration will be implemented in the future. In addition to the explicit microphysics representations, Kuo cumulus parameterization scheme is also available in ARPS for meso- scale prediction.

A proper treatment of the surface processes including momentum, heat and moisture fluxes are equally important for the development and maintenance of storm scale weather and other types of local flow. A stability and roughness-length dependent surface flux model has been implemented in the ARPS using a modified Businger formulation (Businger, et. al., 1971). An analytical procedure is used in the flux calculations for a much improved efficiency (Byun, 1990). Businger's formulation is further modified for more realistic results in highly stable or unstable environments (Deardorff, 1972). The model also distinguishes land and water surfaces for the roughness length calculations. The calculated surface heat fluxes are applied at the ground level in stable conditions, and are linearly distributed throughout the planetary boundary layer (PBL) in unstable conditions. This procedure is similar to the treatment in the PBL parameterization scheme of Blackadar (Zhang and Anthes, 1982). In ARPS, the PBL depth is predicted by a rate equation of Nieuwstadt and Tennekes (1981) for stable boundary layer and that of Gryning and Batchvarova (1990) for unstable boundary layer.

The surface flux model is coupled with a soil-vegetation model after Noilhan and Planton (1989) and Pleim and Xiu (1995). It is designed to simulate the essential processes involved in the surface-atmosphere interactions and requires data of soil and vegetation types at the land surface (Wong, et.al., 1994). In the model, soil surface temperature changes as a result of the net surface energy balance among net radiation, surface sensible and latent heat fluxes, and the deep soil heat transfer. The time rate of change in volumetric soil moisture near the soil surface results from the residual of the precipitation reaching the ground, the evaporation from the ground, and the transfer of moisture between surface and deep soil layers. Separate equations describe the heat and moisture budget in deep soil and the amount of water in the canopy. A simple time-space varying radiation model is used to calculate the net surface radiation. The predicted surface soil temperature and moisture are used in the calculations of surface heat and moisture fluxes. The detailed equations of the soil-vegetation model, as well as the calculations of stability-dependent surface fluxes are described in Xue, et.al. (1995).

2.3. Numerical Treatment

The continuous equations described in Section 2.2 are solved using finite difference method on an Arakawa C-grid. The mode-splitting time integration technique of Klemp and Wilhelmson (1978) is employed. The large time-step integration uses leapfrog time differencing scheme. The advection terms use fourth-order-centered differencing while the other terms use second-order differencing. The small time step integration uses Crank-Nicolson scheme which solves the w and p equations implicitly in the vertical direction. Furthermore, ARPS has the option for evaluating the terms responsible for internal gravity waves in small time steps, which removes the limitation on the large time step size by static stability, thereby allowing bigger large time step size whenever possible. In Eqs.(6)-(10), all terms on the left hand side are evaluated in the small time steps.

Rigid wall, periodic, zero gradient, wave-radiating open boundary and externally forced boundary condition options are available in ARPS. Used by the examples given in this paper is the externally forced condition after Davies (1983), in which all prognostic variables except for w and E are relaxed toward an external solution in a boundary zone. ARPS provides a utility to convert coarse grid fields into boundary condition files for a sub-domain, thereby permitting one-way interactive self-nesting. Furthermore, ARPS has a full implementation of the adaptive grid refinement procedure of Skamarock and Klemp (1993), which is based on a generalized software written by Skamarock and Xue. This procedure provides ARPS with the capability of unlimited-level two-way interactive nesting, while allowing the nested grids to be changed in response to the flow evolution during a run (Xue, et.al., 1993).

Optional upper-level Rayleigh damping can be used to control top-boundary wave reflection.

2.4. Computational considerations

The ARPS computer code is developed under a stringent set of rules. Readability, maintainability and portability of the code have been high priorities during the model development. These virtues, together with extensive internal and external documentation, are perhaps what set this code apart from most of the other application codes. The highly modular design makes adding new packages much easier.

3. REAL-TIME PREDICTION OF CONVECTIVE STORMS

As part of its continuing effort to demonstrate the practicality of storm-scale NWP, CAPS conducted a real-time prediction experiment, using ARPS, during the spring 1995, in collaboration with the VORTEX 95 field experiment. This experiment represents the first time that an atmospheric model initialized with non-homogeneous initial and boundary conditions is used for real-time storm-scale weather prediction, and the model was run on a massively parallel Cray T3D.

To summarize, two 6-hour forecasts were made each day. One of them is on a 1200x1200 km squared domain centered over western Oklahoma, and has a 15 km resolution. Initial and boundary conditions were taken from the forecasts of Rapid Update Cycle valid at 18Z, 21Z and 00Z. The version of ARPS used include the 1.5 order TKE turbulence scheme, Kuo cumulus parameterization, the surface energy budget and soil model, and the stability-dependent surface layer and PBL parameterization. Vertical stretched grid was also used. Each 6-hour forecast took about 20 minutes wall-clock time on 2 dedicated processors of a Cray C90. The second forecast used a 3 km resolution over an 336 x 336 km domain. The location of this grid is targeted at the prime area of severe weather each day, as determined by VORTEX forecasters. This domain is one-way nested within the 15 km resolution domain described above, using the same physics except with Kessler microphysics replacing the Kuo cumulus parameterization. Initial data for this run was from OLAPS analysis (Brewster, et.al., 1994). The 6 hour forecast took about 45 minutes of wall-clock time on a 256 node partition of a Cray T-3D.

Two 3 km resolution forecasts are presented here as examples. The first is a case of multiple squall lines that occurred on May 7. During the day, storms originated from eastern New Mexico in the early morning moved northeastwards across Texas and Oklahoma pan-handles, joined up with a line of cells from western Texas and organized into a connected line of storms that moved eastwards across most part of Oklahoma. The southern end of this line broke into an isolated supercell later in the afternoon and spawn an intense tornado in Ardmore, Oklahoma which caused several casualties. A second line of storms formed at around noon local time in western Texas and propagated east and stretched across Texas, Oklahoma and Kansas. It eventually reached Missouri during its more than 10 hour long life cycle. Fig. 1a shows the observed radar reflectivity from a prominent line of storms at 00Z, May 7, 1995 and Fig. 1b shows the corresponding forecast reflectivity fields in contours (intervall 10 dbZ) together with surface wind vectors. The location of the storm cells associated with the main squall line and the surface wind convergence line agrees well with the observation, although the cells on the southern part were too weak.

On June 8, 1995, another cluster of storms occurred at the triple point between a E-W oriented surface front and a N-S oriented dry line at the Texas and Oklahoma pan-handles. By 00Z, June 9, the storms organized into a SW-NE orientation (Fig.2a), and the storm at the southern end spawned an intense tornado, which was intercepted by the VORTEX field experiment team. Again, ARPS was successful in predicting the location of these storms a few hours in advance (Fig.2b). As is common to most of the forecasts, ARPS is somewhat slow in predicting the initial development of storms, presumably due to the spin-up problem associated with the model old start. This is exactly the type of problem being addressed by the CAPS data assimilation team.
(a)

Figure 1. Observed (a) and predicted (b, in contours) low-level radar reflectivity fields at 00Z, May 08, 1995. Panel (b) also includes the predicted surface winds and a zero contour shown in thick lines. The squared box in panel (a) corresponds to the forecast domain in panel (b).

Figure 2. Same as Figure 1, but for 00Z, June 09, 1995.

4. SUMMARY

The Advanced Regional Prediction System was briefly described in this paper. The results from two high-resolution real-time storm forecasts were presented and compared with radar observations. Additional details on ARPS can be found in Xue, et.al. (1995), and more detailed analyses on the example case will be reported elsewhere. Additional information on ARPS and CAPS can also be found via CAPS WWW home page at http://wwwcaps.uoknor.edu. ARPS Version 4.0 was officially released for public use in August, 1995.

ACKNOWLEDGMENT

Many people contributed to the development and testing of ARPS. Limited by space, their names are not listed here, but can be found in Xue, et.al. (1995).

REFERENCES

Brewster, K., F. Carr, N. Lin, J. Straka, J. Krause, 1994: A local analysis system for initializing real-time convective-scale models. Preprints, Tenth Conference on Numerical Weather Prediction, Portland, Oregon, Amer. Meteor. Soc., 596-598.

Businger, J. A., J. C. Wyngaard, Y. Izumi and E. F. Bradley, 1971: Flux profile relationship in the atmospheric surface layer. J. Atmos. Sci., 28, 181-189.

Byun, D. W., 1990: On the Analytical Solutions of Flux-Profile Relationships for the Atmospheric Surface Layer, J. Appl. Meteorol., 29, 652-657.

Davies H., 1983: Limitations of some common lateral boundary schemes used in regional NWP models. Mon. Wea. Rev. 111, 1002-1012.

Deardorff, J. W., 1972a: Numerical investigation of neutral and unstable planetary boundary layers. J. Atmos. Sci., 29, 91-115.

Deardorff, J. W., 1972b: Parameterization of the Planetary Boundary Layer for Use in General Circulation Models, Mon. Wea. Rev., 100, 93-106.

Droegemeier, K. K., M. Xue, K. Johnson, K. Mills, and M. O'Keefe, 1993: Experiences with the scalable-parallel ARPS cloud/mesoscale prediction model on massively parallel and workstation cluster architectures. Parallel Supercomputing in Atmospheric Science, G.R. Hoffman and T. Kauranne, Eds., World Scientific, 99-129.

Droegemeier, K. K., M. Xue, K. Johnson, M. O'Keefe, A. Sawdey, G. Sabot, S. Wholey, and K. Mills, 1994b: Weather Prediction: A scalable stormscale model. In High Performance Computing, Addison-Wesley. Ed. G. W. Sabot, 45-92.

Durran, D. R., and J. B. Klemp, 1983: The effects of moisture on trapped mountain lee waves. J. Atmos. Sci., 39, 2490-2506.

Germano, M., U. Piomelli, P. Moin, and W. H. Cabot, 1991: A dynamic subgrid-scale eddy viscosity model. Phys. Fluids A 3, 1760-1765.

Gryning, S. E. and E. Batchvarova, 1990: Simple Model of the Daytime Boundary Layer Height. The 9th Symposium on Turbulence and Diffusion, 379-382.

Janjic, Z. I., 1977: Pressure gradient force and advection scheme used for forecasting with steep and small scale topography. Beitr. Phys. Atmos., 50, 186-199.

Johnson, K. W., J. Bauer, G. A. Riccardi, K. K. Droegemeier, and M. Xue, 1994: Distributed processing of a regional prediction model. Mon. Wea. Rev. 122, 2558-2572.

Klemp, J. B., and R. B. Wilhelmson, 1978: The simulation of three-dimensional convective storm dynamics. J. Atmos. Sci., 35, 1070-1096.

Lilly, D. K., 1962: On the Numerical Simulation of Buoyant Convection. Tellus, 14, 168-172.

Lilly, D. K., 1992: A proposed modification of the Germano subgrid-scale closure method. Phys. Fluids A 4, 633-635.

Lin, Y. L., R. D. Farley and H. D. Orville, 1983: Bulk parameterization of the snow field in a cloud model. J. Clim. Appl. Meteor., 22, 1065-1092.

McGinley, J. A., S. C. Albers, and P. A. Stamus, 1991, Validation of a composite index as defined by a real-time local analysis system. Wea. Forecasting, 6, 337-356.

McGinley, J. A., 1995, Opportunities for hugh resolution data analysis, prediction and product dissemination within the local weather office. Preprints, 14th Conf. on Wea. Analysis and Forecasting, Amer. Meteor. Soc., 478-483.

Moeng, C.-H., 1984: A large-eddy-simulation model for the study of planetary boundary-layer turbulence. J. Atmos. Sci., 41, 2052-2062.

Nieuwstadt, F. T. W. and Tennekes, H., 1981: A rate equation for the nocturnal boundary-layer height, J. Atmos. Sci., (38) 1418-1429.

Noilhan, J., and S. Planton, 1989: A simple parameterization of land surface processes for meteorological models. Mon. Wea. Rev., 117, 536-549.

Ogura, Y., and N. A. Phillips, 1962: Scale analysis of deep and shallow convection in the atmosphere. J. Atmos. Sci,. 19, 173-179.

Orlanski, I., 1976: A simple boundary condition for unbounded hyperbolic flows. J. Comput. Phys., 21, 251-269.

Pleim, J. E. and A. Xiu, 1995: Development and testing of a surface flux and planetary boundary layer model for application in mesoscale models. J. Appl. Meteorol. 34, 16-32.

Sharman, R. D., T. L. Keller, and M. G. Wurtele, 1988: Incompressible and anelastic flow simulations on numerically generated grids. Mon. Wea. Rev., 116, 1124-1136.

Skamarock, W., and J. B. Klemp, 1992: The stability of time-splitting numerical methods for the hydrostatic and non-hydrostatic elastic systems. Mon. Wea. Rev., 120, 2109-2127.

Skamarock, W., and J. B. Klemp, 1993: Adaptive grid refinement for two-dimensional and three-dimensional nonhydrostatic atmospheric flow. Mon. Wea. Rev., 121, 788-804.

Smagorinsky, J., 1963: General Circulation experiments with the primitive equations. Mon. Wea. Rev. 91, 99-164.

Tao, W. K., and J. Simpson, 1993: Goddard Cumulus Ensemble Model. Part I: Model description. Terrestrial, Atmospheric and Oceanic Sciences, 4, 35-72.

Tao, W. K., J. Simpson, and M. McCumber, 1989: An ice-water saturation adjustment. Mon. Wea. Rev., 117, 231-235.

Tripoli, G. J., and W. R. Cotton, 1982: The Colorado State University three-dimensional cloud/mesoscale model - 1982. Part I: General theoretical framework and sensitivity experiments. J. Rech. Atmos., 16, 185-219.

Wong, V. C., 1992: A proposed statistical-dynamical closure method for the linear or nonlinear subgrid-scale stresses. Phys. Fluids A 4, 1