Method of Solving Moist Thermodynamic Equations in NTU-Purdue Non-hydrostatic Model and Tests on 2D Moist Mountain Waves

A diagnosing method to solve pressure and temperature with the pre­ dicted values of equivalent potential temperature, total water content and density is developed for the NTU-Purdue Nonhydrostatic Model. The effi­ ciency and accuracy of the method is examined for a wide range of atmo­ spheric conditions. We found that very few iteration steps are required because of the fast conver gent rate of the method. The diagnosin g method is tested on moist mountain waves, and the results are consistent with our present understanding of the phenomenon. The model seems to be extremely. stable, as it endures very long integration without showing signs of deterio­ ration that is common in a numerical solution. It is also very accurate in depicting many detailed structures of the mountain waves.


INTRODUCTION
Purdue Mesoscale Model (PMM) developed by Dr. Wen-Yih Sun's group at Purdue Uni versity has been proved accurate and useful in many different research topics. These research works include studies of shallow convection (Hsu and Sun 1991), Mei-Yu front (Hsu and Sun 1994), mountain lee vortices (Sun andChem 1993, 1994) and regional climate, etc. One of the major factors attributing to its success is the choice of thermodynamic prognostic variables in its governing equations. The bases of the numerical atmospheric model are, of course, the conservation laws of mass, energy, and momentum. Many models use temperature (Dudhia 1993;Juang 1992) or potential temperature (Durran and Klemp 1983) as one of the prognostic variables. Yet there is no conservation law that would apply directly to temperature and po tential temperature. Additional terms must be calculated to include the heat changes due to moist processes. This is also the case for pressure. The prognostic equation for pressure has been derived from other conservation laws under some approximation of moist thermodynam ics, which usually involves a definition of equivalent potential temperature (e.g., Wilhelrnson 1 Department of Atmospheric Sciences, National Taiwan University, Taipei, Taiwan, ROG 1977). PMM, however, predicts equivalent potential temperature and total (liquid plus vapor) water content. They have the nice property of being nearly conservative. Modification is required only for precipitation process and deep convection. Temperature, specific humidity and other intensive properties are determined diagnostically by moist thermodynamics for this approach. Deardorff (1976) first developed this method in his shallow-cloud model. Ooyama (1990) extended the concept further; he presented a thermodynamic foundation for modeling the moist atmosphere, which includes the ice phase in the formulation. He claimed that his method returns to the 'primitive' form of moist thermodynamics, in which prediction is made strictly in terms of conservative properties, such as mass and entropy.
A joint effort by the first author of this paper and Dr. Wen-Yih Sun of Purdue University has been made to develop a fully compressible, nonhydrostatic numerical model. This Na tional Taiwan University (NTU)-Purdue Nonhydrostatic Model is mainly derived from the original hydrostatic PMM. The two models are similar except that two additional prognostic equations for density and vertical momentum are added in the non-hydrostatic model. ·Pres sure field in a hydrostatic model can be easily obtained with the calculated temperature field and boundary values of pressure through the hydrostatic equation. However, for a compress ible, nonhydrostatic model, pressure field is either directly integrated or indirectly retrieved from the integrated density field. With a new set of prognostic thermodynamic variables, the method for diagnosing temperature field must be redesigned and pressure must also be deter mined in a different manner for the new NTU-Purdue Nonhydrostatic Model.
This paper presents a method of retrieving temperature and pressure from the predicted values of equivalent potential temperature, total water content and density.
We used Newton's iteration method to solve the set of moist thermodynamic equations, instead of using empirical formula as in Deardorffs (1976) model and PMM. The empirical formula is usually valid for only limited range of temperature and pressure, The accuracy is often lost in very high or very low temperature condition. The Newton's method, however, is more accurate in these situations. The accuracy of our newly developed method is investi gated and the method is tested on moist mountain waves.

GOVERNING EQUATIONS
The NTU-Purdue Nonhydrostatic Model uses a terrain-following <5 coordinate. a is defined as: where p0, pressure of a reference atmosphere, is strictly a function of height. The six prognos tic equations are: (2) av av .av The vertical velocity in (j coordinate is diagnosed as: . aa (aa J ifz ax , • and p', defined asp-pJz), can be diagnosed through: p' = function of ( ee, q w, p).

307
(3) (5) The purpose of this paper is to show the method of solving (9) in our model. Coriolis force is not calculated due to the small scale of the mountain. Since the present study aims at linear waves, which do not break into small eddies, diffusion term are also ignored. The symbols are conventional.

MODEL
Since the emphasis of the present paper is on the moist thermodynamic equations, in addition, the details of the integrating schemes for dry conditions are already presented in Hsu and Sun (1998), we shall only repeat some of most basic numerical framework here for conve nience to the readers. The compressible equations of motion (2)-(5) include fast propagating sound waves. It imposes severe limit to the time step that can be used in explicit numerical integration schemes. Our model is explicit. In order to save computing time, we employ a time-splitting (Gadd 1978) technique. The integration is divided into two stages, in which the length of the time step varies according to the corresponding physical time scale of the terms calculated.
The terms that are responsible for rapid sound wave propagation are isolated on the left hand side of the equations (2)-(5). They are integrated separately in the sound-wave stage with a smaller time step than that used for other processes. We use forward integration for momentum fields and backward integration for density.
The other terms are calculated in the gravity wave stage using a much larger time step.
Another forward-backward procedure is devised in gravity-wave stage as well. Thus, there is no unstable computational mode resided in the whole model. The algorithm for the advection terms follows Sun (1993). This advection scheme is more accurate than the Crowley 41h order scheme (Crowley, 1968). As for the boundary conditions, the lateral conditions are opened for both sides, the top 1/3 of the domain is a sponge layer to reduce reflection. The ground is the only physical boundary. An inviscid condition is imposed.

METHOD FOR RETRIEVING PRESSURE AND TEMPERATURE
Pressure and temperature are determined with the predicted values of (} , q and p. For (15) where q . , the saturated specific humidity, is expressed as a function of pressure and tempera ture with the aid of Tetens' formulae (Tetens 1930). The two quantities, q w and q 1 in (13) are total water content, liquid water content, respectively. There are three constants in (11). The constant, a, is set to be 0.98 to allow condensation in a grid volume where the relative humid ity has not yet reached 100%. The constant, e, is Rd /Rv and e,0 is the saturate vapor pressure at triple point. This set of equations can be manipulated to yield two nonlinear equations: and (17) where a, band care known with the predicted values of density, equivalent potential tempera ture and total water content: By evaluating the partial derivatives ofjand g with respect top and Tseparately, (16) and (17) can be solved using Newton's method. · As for the unsaturated condition, q 1 vanishes and q , in (12)-( 14) is replaced by q w . This set of equation can be rearranged into a single nonlinear equation: where h is known: Potential temperature in (21) can then be solved using Newton's method. Once we have obtained potential temperature, pressure and temperature can be easily calculated though (10), (12) and (14).

ACCURACY TESTS
The retrieval method presented in the last section is tested for a wide (both pressure and temperature) range of atmospheric thermal states. We have chosen 20 constant-lapse-rate atmospheric profiles of temperature as functions of pressure. The lapse rate is 6.5 K km·' and the temperature at 1000 mb level varies from 284 to 303 K for the 20 profiles. Since the algorithms for saturated and unsaturated air should yield the same result for marginally satu rated condition, the same thermal states are tested for both methods assuming saturated air with no liquid water. Figure 1 shows the distribution of temperature and saturated specific humidity of all 20 profiles for the test.  We first calculate density, equivalent potential temperature and saturated specific humid ity from the chosen thermal states. Then, both temperature and pressure are offset deliberately by 1 % of their true values as the first guess of Newton's method with known values of density, equivalent potential temperature and saturated specific humidity. That is equivalent to initial errors of 1 to 10.5 mb in pressure and 1.8 to 3.1 Kin temperature. The results are shown in Figures 2 and 3. The convergent rate is very fast for the unsaturated case ( Figure 2). The error in temperature reduces to less th an 0.002% in just one calculation. With two iterations, the error is close to the machine truncation error, even though the calculation is performed in double precision. The error in pressure (not shown) is very much the same. As for saturated air (Figure 3), the convergent rate is not as fast. Still, even in the high-temperature and high pressure condition (lower right-hand-side comer of the diagram), each iteration reduces error by nearly 10% of the error of the previous iteration. The method is very accurate at least for the thermal states of our concern. The number of iterations required in this method depends on how fast the thermal state changes with time and how short the integration time step is used.
For the moist mountain waves simulation in the next section of the paper, only one calculation is sufficient to achieve the desired accuracy.

SIMULATION OF MOIST MOUNTAIN WAVES
In order to show the latent-heat effect and test the diagnosing method of the NTU-Purdue  ground wind velocity (10 m s-1) and uniform static stability (Brunt Vaisala frequency of 10-2 s-1) are imposed for simplicity. The air is forced to pass over a very narrow bell-shaped ridge with half-width length of 1 km and maximum height of 10 m. Queney (1948) solved the corresponding mountain waves solution analytically, and we can therefore validate our nu merical model in this particular simulation. Grid interval is 200 m in both the x and z direc tions. Open lateral boundary and inviscid bottom boundary conditions are assumed. Figure 4 shows the vertical velocity distribution after 10 hr of integration. Short lee waves appear and they propagate upward with little dissipation in energy except in the upper portion of the domain, where strong damping terms are artificially added to avoid reflection from the top boundary. The property of the upward transport of wave energy can be demon strated in Figure 5. In the lowest 12 km, wind stress varies only very slightly with height, as should be the case for small-amplitude, vertically propagating mountain waves (Eliassen and Palm 1961). Wind stress reduces to near zero in the upper damping layer as expected. This figure also shows that simulated waves approach near steady state in about 4 hr, and the nu merical model is well behaved. Moisture and liquid water are added in the lowest 10 km of the domain for moist moun tain waves simulation. Initial liquid water content is 0.2 g kg-1, and the other environmental conditions are the same as in the dry run. The moist mountain waves are much weaker ( Figure   6). This phenomenon is due to the reduction of vertical stability of the disturbed upward/ downward airflow by the condensation/evaporation process. The reduction can be shown quantitatively by the moist Brunt Vaisala frequency, Nm: Figure 7 shows the distribution of Nm2 with height. It is very small close to the surface, where the latent-heat effect is most pronounced.

km
In addition to amplitude, the wavelength is also affected. For the linear, dry mountain waves, the wavelength equals to 21n1/N (Queney 1948), which corresponds to 6.28 km for the present settings of u and N. Our model results show consistency (almost perfectly) with the theoretical value. The vertical wavelength of the mountain waves shortens with decreasing degree of stability. The location, marked 'A' in Figure 6, where ve,rtical velocity changes to positive value above the top of the mountain, is at the level of 5.4 km. As for the dry simula tion (Figure 4), the corresponding point is only 4.2 km above the top of the mountain. The wavelength is larger for the moist run near the surface due to the reduction of the vertical stability. These characteristics of the moist mountain waves have been discussed in Durran and Klemp (1983) and Jusem and Barcilon (1985).  Labels of x-axis are in 10·5s·2•

315
The amount of water that changes phase due to the gravity waves process is depicted in Figure 8. Its distribution is consistent with perturbation virtual potential temperature field shown in Figure 9. The amount of liquid water is larger than its initial value of 0.2 g kg-1 for the negatively buoyant air parcel with neg ative virtual potential temperature perturbation. The· reverse is also true for the downward motion, since warmer temperature favors more evapora tion from the initial cloud water. Another interesting feature is that although the magnitude of the temperature perturbation is about the same at all vertical levels, but the amount of water condensed or evaporated is much greater at a lower level than at an upper level. This feature is shown by the packed contour lines of q1 near the mountain in Figure 8. The simulation demonstrates that the diagnosing method described in Section 3 is very effective in solving intensive properties of the air with the presence of the moist process. The contour lines in Figure 8 are very smooth, even though very small plotting intervals are used.
This implies that there are no undesirable noises induced by our newly developed numerical method. In addition, the numerical solution is very stable (Figure 10), it remains steady and accurate for very long integration. The method will enable the NTU-Purdue Nonhydrostatic Model to become a useful tool in studying a wide range of atmospheric problems.
201�mmnmnm1 1l l1 1Tl llm mm mmrm mmrnmr rm mm mTm mmTm mmTnm mr rm mrm mml ll1 11 11l lnm mn nnm mn nm mm mmr rm mmrm mmrm mmrm mml ll1l ll1l"[l l km H.oooa In order to verify our numerical method, a case of moist, linear, hydrostatic mountain waves simulation presented in Durran and Klemp (1983) is repeated. The saturated (but no clouds) air flows over a bell-shaped mountain with a half-width length of 10 km. The back ground conditions are uniform with N=0.0132 s·• and U=20 m s·1• The surface potential tern-

317
perature is 285 K and pressure is 1020 Pa. In this simulation the domain contains 90 points in the horizontal as in Durran and Klemp (1983), and 91 vertical grid points (25 grids more than theirs). The grid intervals are the same (2 km in x-direction and 333 min z-direction).
The results are presented in Figures 11 and 12. Figure 11 shows Durran and Klemp's (1983) results at the 30,000 s. We show the results at the 1Qth hr for comparison. Both models reach a near steady state after a very long integration. It is evident that we can fairly well reproduce their results. Both the magnitude and pattern of the mountain waves matches very closely. Our method is also valid in the hydrostatic situation.

CONCLUDING REMARKS
We have developed a diagnosing method to solve pressure and temperature with the pre dicted values of equivalent potential temperature, total water content and density for the NTU Purdue Nonhydrostatic Model. After arranging the system of equations into only two nonlin ear equations for the more complicated saturated condition, we can solve them without too much trouble. The efficiency and accuracy of the method is examined for a wide range of atmospheric conditions. We found that very few iteration steps are required because of the fast convergent rate of the method.
The diagnosing method is tested on moist mountain waves, and the results are consistent with our present understanding of the phenomenon. The model seems to be extremely stable, as it endures very long integration without showing signs of deterioration that is common in a numerical solution. It is also very accurate in depicting many detailed structures of the moun tain waves. The method can be applied in many research problems related to strong latent heating effects. The present version of the model, however, does not treat ice, which can be quite important in some problems. Future works are needed in that area.  Durran and Klemp (1983). The back ground air is saturated. Cloudy regions are shaded. The perturbations have been multiplied by 1000 for display; as such they constitute a lin ear numerical solution for a 1 km high mountain. The contour interval is 2m s-1•