II. 14.
COUPLED RAINFALL-FLOW FORECASTING MODELS
Konstantine P. Georgakakos
Almost ten years ago Georgakakos (1986a-b) showed that coupling a hydrologic flow forecast model with a rainfall prediction model results in better forecasts of flood peak magnitude and timing, and in better overall skill as measured by quadratic performance indices. The convective rainfall prediction model he used was based on the cloud- and rain-water mass balance and on microphysical parameterizations (Georgakakos and Bras, 1984a-b). The hydrologic model used consisted of the Sacramento soil-water balance model (e.g., Peck, 1974) and of a kinematic flood routing model (Georgakakos and Bras, 1982). The coupled model was a stochastic-dynamical one in that it generated state updates for all components from available observations of rainfall and flow discharge in real time. A state estimator suitable for use with nonlinear dynamical formulations was used (i.e., a form of the Extended Kalman Filter). The improved performance of the coupled model was attributed to the generation of skillful mean areal rainfall estimates for the study catchment by the rainfall prediction component. It was also attributed to the ability of the stochastic-dynamical coupled model to produce improved state updates using both rainfall and flow observations in real time. Later, Georgakakos (1987), using data from several catchments in the U.S., showed the high utility of such coupled models for the real time prediction of flash floods when the catchment response time is short. Georgakakos and Foufoula-Georgiou (1991) elaborated further on this issue by showing that the utility of coupled models diminishes as the ratio of the catchment response time to the forecast lead time increases beyond the value of 1. The same authors also showed that apart from the direct coupling of the rainfall and flow formulations due to mass conservation, there is indirect coupling in real time due to the updating procedure, in particular, due to non-zero state estimator gains. They further show that such non-zero gains exist because of the physical requirement for water continuity at the interface between the atmosphere and the land surface. Finally, recently Bae et al. (1995) confirmed the good performance of a stochastic-dynamical coupled rainfall-flow model in an operational environment. They showed that, for short-term flow prediction, the particulars of the soil-water balance formulation do not affect significantly the performance of the model, as long as the parameters of the formulation are reasonably well estimated from historical data.
Although implicit in previous work as outlined above, the effects that errors in meteorological parameters have on the predictions of flow and flow variance has not been demonstrated. The purpose of the paper is to use coupled rainfall-flow formulations in numerical experiments to show the influence of certain thermodynamical and microphysical meteorological parameters on the channel flow estimates during storm periods. In particular, the influence of meteorological-parameter errors on flow predictions is studied. For the generation of numerical results and in lieu of actual time series, the coupled model is driven by stochastic input processes simulating convective available potential energy (CAPE) and precipitable water (in Monte Carlo numerical experiments), which preserve reported statistical properties of the corresponding natural processes. Apart from its theoretical interest, this study is pertinent to the development of design criteria for establishing hydrometeorological observation networks in support of real-time hydrometeorological prediction systems. All the analyses presented are for an idealized tropical storm-catchment system and the model components are physically-based and observable to allow inference for other catchments with similar parameters.
Section 2 describes the formulation of the rainfall, soil water, and channel flow model components. Section 3 contains the sensitivity analyses and discusses important results. Section 4 presents concluding remarks.
The components to the coupled rainfall-flow model to be studied are: (a) a rainfall generation component, (b) a soil water mass balance component, and (c) a kinematic routing component for channel flow propagation and attenuation. No snow computations are made, as numerical results are for an idealized tropical situation with rainfall-driven floods. The formulation of the model components presented here differs substantially from those used in previous applications. This is mainly due to the different use of the coupled model here to allow for realistic simulations of a particular system and to produce results that may be later compared to those obtained for other physical situations. Parameterizations depend on quantities which have physical meaning and in most cases are observable.
Rainfall
A tropical cluster is used as the idealized tropical storm. As described in Ogura (1986), a tropical non-squall cluster is a common tropical storm, with the squall cluster being a more severe, longer-lived but rarer phenomenon. For example, during the Global Atmospheric Research Program's (GARP) Atlantic Tropical Experiment (GATE), which took place over the eastern Atlantic and the African continent during the period June-September 1974, the vast majority of cloud clusters observed was of the non-squall type. Figure 1 shows a schematic representation of a cluster, indicating the convective core region and the trailing anvil region. The arrows indicate the relative path of moist air of high moist static energy as it enters the storm and feeds the anvil updrafts and the downdraft region of intense convective rainfall. The motion of the storm system is from left to right in the Figure. Typical pressure heights and approximate horizontal scales are also indicated. It is noted that the structures of non-squall and squall clusters possess similar features. A principal difference is the sharper distinction between the stratiform and the convective core region. Such a sharper distinction is primarily due to the enhancement of the convective downdraft by dry environmental air from levels between 900 to 600 mbar. Another important difference for hydrologic applications is that the non-squall clusters possess generally much slower propagation velocities than squall clusters, and in such a way are more likely to generate flood-producing rainfall in tropical catchments. For this paper we use a generic formulation for both types of clusters.
The mathematical formulation for the idealized storm of Figure 1 is based on the mass balance of the condensed water equivalent for the convective and the anvil regions, taken separately. The anvil region is modeled separately because of its importance in producing substantial volume of precipitation (studies reviewed by Ogura, 1986, report total volume of anvil precipitation that is 30-50 % of the total storm precipitation). Typical microphysical and convective-thermodynamic parameterizations are used to express the water mass fluxes as functions of the cloud- and rain-water states.
Fig 1 - Schematic of an Idealized Non-Squal Cluster Tropical Storm (Features after Ogura, 1986)
Denote by X the mass of condensed water equivalent per unit area in the convective region and by Y the analogous quantity in the anvil region. Also, denote by Ic the influx of condensed liquid water equivalent into the storm system, by Oc the precipitation flux in the convective downdraft region assuming negligible sub-cloud layer, by Oa the precipitation flux at the anvil base, and by Mc and Ma the water mass-loss flux due to mixing with environmental air in the convective and anvil regions, respectively. With these definitions, the expressions of mass balance for cloud- and rain-water in the convective and anvil regions are:
| dXdt = k[(1-b) Ic - Oc - Mc] | (1) |
dYdt = bk'Ic - Oa -Ma (2)
with b denoting the fraction of the condensation influx which flows into the anvil region from the convective region, with k denoting the ratio of the area of the region of condensation influx to the area of the region of convective downdraft (in Figure 1 it is 20/20=1), and with k' denoting the ratio of the area of the region of condensation influx to the area of the region of the anvil (in Figure 1 it is 20/40=0.5). The mean areal precipitation P at the surface over a catchment is then:
P = a1Oc + a2Pa (3)
with a1 and a2 being the fractions of the catchment area covered by the convective core and the anvil regions, respectively, and with Pa being the precipitation rate at the surface due to anvil precipitation processes. The later rate differs from Oa because of the evaporation in the substantial sub-cloud layer under the anvil (estimated height of sub-cloud unsaturated air is about 4 km for the idealized storm of Figure 1).
The parameterizations of the various fluxes in Eqs. (1)-(3) follows earlier work by Georgakakos and Bras (1984a-b) (heretofore called GB84), with microphysical formulations having different parameter values for the convective core and the anvil regions. In the following formulas, subscript i denotes either the convective core (subscript c) or the anvil (subscript a) region. Thus, the updraft velocity for region i is given by
wi = ei ÷`2Ec (4)
with Ec being the convective available potential energy of surface moist air, which is responsible for the convective updraft. For the convective region, denote by Dwc the difference between the surface mixing ratio minus the cloud-top mixing ratio, and rc the average air density in the clouds. Then,
Ic = wc Dwc rc (5)
The size distribution of the hydrometeors in both the convective core and anvil region is assumed exponential:
Ni(D) = No e-D/Di (6)
where, D is the hydrometeor effective diameter, Di is the average diameter, and No is a measure of the number density of very small drops. The terminal velocity of precipitation particles is taken as vT = cD, with c being a constant estimated from empirical data (see GB84). The precipitation rate Oc in the convective region may then be estimated to be:
Oc = 4cDcXHc (7)
with Hc denoting the height of the convective-core cloud cluster. The result is obtained assuming invariant drop-size distribution parameters with cloud height. Following the derivation of GB84, the precipitation rate at the base of the anvil region is estimated to be:
Oa = 4cDaYHa [ 1 + 3/4 Nv + 1/4 Nv2 + 1/24 Nv3eNv ] (8)
where Ha denotes the cloud height in the anvil region, and Nv is a dimensionless parameter defined as:
Nv = wacDa (9)
and is an indicator of the mesoscale updraft strength in the anvil region. Note that the flux Oa is computed below the melting layer so that the terminal velocity of the drops can be characterized by the same parameter c used for the convective region. The precipitation flux in the anvil region, Pa, is given by:
Pa = 4cDaYHa [ x(NdNv ) (1 - 1/4 Nv)(1 + Nd + 1/2 Nd2) + 1/8 Nd3eNd +
+ (1 - x(NdNv ) ) 1 + 3/4 Nv + 1/4 Nv2 + 1/24 Nv3 - 1/24 Nd3eNv ] (10)
with x(u) denoting the step function being equal to 1 for u≥1 and equal to 0 everywhere else. The dimensionless number Nd is an indicator of the strength of evaporation of drops in the subcloud layer of the anvil region and it was shown by GB84 to be equal to:
Nd = DoDa (11)
where Do is the diameter at the anvil-cloud base of the drop which completely evaporates upon reaching the ground. Do may be expressed as a function of the atmospheric state variables for subcloud air as done in GB84. However, this requires knowledge of the temperature and humidity profiles of the subcloud layer. In lieu of such time series, the parameter Do is set to a constant value corresponding to a reference state of the ambient atmosphere, which largely makes the sub-cloud air in the anvil region. On the basis of results in Pruppacher and Klett (1980) and for the tropics, the diameter Do is fixed to 0.2 mm, assuming a relative humidity equal to 80 % for the reference state.
To close the formulation of the cloud- and rain-water mass balance, the mass-mixing terms Mi (i=c,a) are parameterized in analogy to the momentum mixing terms in Adler and Mack (1986):
Mi = mi Xi max{|wi|,cDi} (12)
where mi is a mass-mixing parameter, Xi denotes X for i=c and Y for i=a, and max{ } denotes the largest (in magnitude) of the two terms in brackets.
In summary, the input of the rainfall component consists of the time series of Ec and Dw, and the output is the mean areal precipitation rate P. It is noted that for SI units, the state variables of the precipitation formulation may be expressed in kg/m2 or equivalently, in mm of water depth, and the precipitation flux P in mm/hr by multiplying P of (3) by the constant 3600.
Several physically-based soil water models have been proposed in the literature. They range from those that emphasize soil physics and plant transpiration processes for climatic studies (e.g., Mintz and Walker, 1993) to those that use parameterizations suitable for forecasting applications over natural catchments (e.g., Georgakakos, 1986a). In the present study we use a one-layer soil water balance formulation, which is most suitable for mountainous catchments with thin soils and steep slopes, prone to flash flooding.
The water balance expression gives the model state (Z - water stored in soil column) as a function of the water influx (P - precipitation) and outflux (E - evapotranspiration, R - surface runoff, Bf - baseflow, and L - seepage loss to deep aquifers):
dZdt = P - E - R - Bf - L (13)
with the following parameterizations,
E = Ep ZZo (14)
R = P ( ZZo )m (15)
Bf = g 11+l Z (16)
L = g l1+l Z (17)
where Ep is the potential evapotranspiration rate, Zo is the soil water capacity, m is a coefficient for surface runoff generation, g is the baseflow recession rate, and l determines the fraction of baseflow lost to deep aquifers. The rate Ep is a function of the atmospheric forcing (i.e., net radiation, near-surface wind speed, near-surface relative humidity, etc.) and of the vegetation state (e.g., growing season). In hydrologic simulations, Ep is either estimated from pan evaporation data and then adjusted monthly with empirical coefficients (e.g., Peck, 1976 for the Sacramento model), or it is based on combination energy-aerodynamic formulas which contain free parameters and which are estimated using historical data (e.g., Brutsaert, 1991). During storm periods with heavy rain (e.g., convective region of storm in Figure 1, Ep is small due to saturated near-surface air with respect to water vapor, and thus, E may be neglected. For inter storm periods or for periods with moderate and light rain, Ep my be significant and cannot be neglected a priori. For the case study of this paper, Ep was assumed constant and equal to 4 mm/day.
The input of the soil water model consists of the time series of P. The output total channel inflow is computed from:
uT = R + Bf (18)
A physically-based kinematic flow routing model is used to simulate flow in upland catchments with moderate to steep slopes. The main channel is subdivided in channel reaches. Each reach is assumed to possess a wide rectangular cross-section of top width B, hydraulic depth Dh, local bottom slope So, and with a local roughness coefficient n. The length of the reach is denoted by L. Note that the wide rectangular shape was chosen for simplicity in the formulation, but other shapes may be used when appropriate (e.g., triangular, trapezoidal). For length L small compared to the flood wavelength (mild assumption for natural channels), the water volume S in the channel is:
S = BDhL (19)
The cross-sectional flow velocity is well approximated by Manning's formula in uniform flow (see Henderson, 1966, pg. 90-101), and the outflow U from the reach is:
U = So1/2 Dh5/3 Bn (20)
From (19) and (20) it follows:
U = d Sd (21)
with
d = 1n So1/2B2/3L5/3 (22)
d = 5/3 (23)
Expression (21) characterizes kinematic routing models for which the energy slope, the water surface slope and the channel bottom slope are equal, and the relationship between water depth and discharge is a single valued function.
The water balance relationship for reach i (i=1,2,...N) in the main channel is:
dSidt = Ui-1 + pi uT - Ui (for all i and for Uo = 0) (24)
with pi being the fraction of total channel inflow that corresponds to the local area drained by channel reach i, and Ui-1 denoting the upstream inflow to channel reach i from channel reach i-1. The channel flow model input is the total channel inflow uT computed by the soil water component. The catchment outflow, UN, is the outflow of the last downstream reach.
In lieu of actual data for particular natural storm-catchment systems and to obtain statistically significant results, baseline runs of the coupled model were made driven by stochastic process models of Ec and Dw, and with nominal values of all model parameters. The input stochastic-process models were constructed to preserve observed statistical features of relevant observations. Thus, Ec was modeled as a periodic process with given period, and amplitude, and with superimposed correlated noise of a given mean, variance and lag-1 correlation coefficient. Periodicity in Ec has been observed in the few past studies in existence, but with substantial noise superimposed (e.g., Fig. 4 in Zawadzki et al. 1981 with data from Canada; Figure 8.12 in Ogura, 1986 using tropical storm data from MONEX; Fig. 2 in Randall and Wang, 1992 using GATE data from the tropics). The input Dw was modeled in a similar fashion. Data sets for Dw are very sparse in the tropics, and the form of the stochastic process was based on a few reported observations from both tropical and extra tropical regions. Characteristic of these observations is the variability of the water vapor mixing ratio over long (months) and short (minutes) time scales. Relevant results are reported by Hogg et al. (1983, for Denver, Colorado), by Ogura (1986, for the tropical experiment MONEX), by Hense et al. (1988, for several tropical regions) and by Wang et al. (1995, for Wallops Island, Virginia, U.S.A.). Preliminary results from the TOGA-COARE experiment (e.g., Webster and Lukas, 1992) show a similar character of variability. Figure 2 presents plots of observed Ec and Dw for ten days in February 1993. The interval between observations was in most cases 3 hours. The soundings were taken
Fig. 2 - CAPE and Mixing ratio difference time-series data from TOGA-COARE.
from the ship Natsushima while it remained between latitudes [0o0' - 0o47'N] and longitudes [155o56'E - 156o03'E]. The Ec and Dw data exhibit a remarkable degree of co-variability. There are cycles in both data sets that span the period of a few days with fluctuations of a few hours duration.
The mathematical formulation of the forcing was as follows:
Ec(ti) = Eco + (Eco + eE(ti)) sin(2ptiTE ) and Ec≥0 (25)
with
eE(ti) = rE eE(ti-1) + sE÷`1-rE2 hE(ti) (26)
and
Dw(ti) = mw + (mw + ew(ti)) sin(2ptiTw ) and Dw≥0 (27)
with
ew(ti) = rw ew(ti-1) + sw÷`1-rw2 hw(ti) (28)
where, 2Eco is the mean amplitude of Ec's periodic fluctuations, TE is the period of the Ec periodicity (e.g., 1 day for diurnal cycle computations); rE and sE are the lag-1 correlation coefficient and the standard deviation of the Ec-amplitude random fluctuations; 2mw, sw, and rw are the mean amplitude of Dw, the standard deviation and lag-1 autocorrelation coefficient of Dw-amplitude random fluctuations; and hE and hw are zero-mean, variance-one, Gaussian random deviates. The discrete time ti is used with a 5-minute interval between generated values for both hE and hw. It is noted that for the experiments reported here, it was assumed that there was always a trigger mechanism to initiate convection as long as CAPE was available. This somewhat restricting assumption may be lifted by using a random arrival process of trigger mechanisms superimposed on the processes of (25) and (26). Nominal values of the parameters, which were used in the numerical experiments, are listed in the following. The models proposed offer much flexibility in the specification of parameters and this is but one example. It is also possible (and advisable in several cases) to impose random fluctuations on selected parameters (e.g., the average convective and anvil hydrometeor diameters) to study the effect of parameter estimation errors to coupled-system states and fluxes.
Forcing: TE=Tw= 4 days; Eco= 600 J/Kg; sE= 1500 J/Kg; rE=rw= 0.9 (1-hr lag); mw= 0.010 Kg/Kg; sw= 0.015 Kg/Kg.
Rainfall: b= 0.30; k= 1; k'= 0.5; a1= 0.3; a2= 0.7; rc= 0.6 Kg/m3; ra= 0.4 Kg/m3; ec= 0.5; ea= 0.05; Dc= 200.10-6 m; Da= 100.10-6 m; Do= 200.10-6 m; Hc= 12.103 m; Ha= 7.103 m; c= 3500 s-1; mc= 10-7 m-1; ma= 5.10-7 m-1.
Soil Water: Zo= 50 mm; Ep= 4 mm/d; m= 8; g= 0.005 d-1; l= 0.
Channel Flow: N=1; p1= 1; So= 0.005; L= 50.103 m; B= 50 m; n= 0.05
The numerical experiments had initial conditions: X(to)= 0; Y(to)= 0; Z(to)= 40 mm; S(to)= 106 m3. The area of the idealized catchment was estimated as a function of channel length from the stable empirical relationship: A(mi2) = [L(mi)/1.40]1.76. Three months of data were generated for all cases, and analysis was focused on the last month to avoid initial transients. Figure 3 presents the rainfall (mm/day) and the channel outflow (mm/day) for the nominal parameters. The nonlinear nature of the rainfall-runoff relationship is apparent.
Fig. 3 - Rainfall and channel flow for the idealized storm-catchment system using nominal parameters.
Results from two types of sensitivity experiments are presented. In both cases, sensitivity was measured by the root mean quadratic difference from the nominal channel flow trajectory. During the first class of experiments, parameters Dc and ec were perturbed from their nominal values. During the second class of experiments the standard deviations sE and sw were perturbed. The first class of experiments illustrates the effect of microphysical-parameter errors on flow predictions. The second class of experiments shows the effect that erroneous measurements of the temporal variability of Ec and Dw have on flow predictions. This last class of experiments is pertinent to determining the allowable measurement error in the thermodynamical input for achieving a set root mean square error in flow prediction. Results of these sensitivity experiments are presented in Figures 4a-4d. For each case, the figures show the root mean square difference from nominal flow, QRMSE, in mm/day. It is apparent that errors in the microphysical parameters Dc and ec
Fig. 4 - Root mean square difference from nominal flow for indicated changes in parameters: (a) Dc, (b) ec, (c) sE, and (d) sw.
produce greater QRMSE than those in the thermodynamical parameters Ec and Dw. (To establish a reference for the magnitude of the QRMSE values, the reader is referred to Figure 3.) While the microphysical parameters and sw result in very-nearly symmetric error behavior, asymmetry is observed in the case of sE, with underestimation being worse than overestimation.
Past studies using coupled meteorological-hydrological models in real time prediction were reviewed. Such models have demonstrated good performance when used with actual observations and when they are implemented in an operational environment, especially for flash flood prediction. It is suggested that coupled models can also be used to model natural storm-catchment systems in order to establish sensitivity of flow predictions to individual meteorological and hydrological parameters. Their use in the design of hydrometeorological measurement networks is also advocated. As an example, a flexible, physically-based formulation of a tropical storm-catchment system was used to elucidate some of these ideas. It was shown through numerical experimentation, that the formulation is capable of producing realistic simulations. It was further shown that estimation of microphysical parameters (such as the average drop diameter in the convective-core region of the storm) is very important for reliable flow prediction; more so than the thermodynamical input variables of CAPE and water vapor mixing ratio.
Future studies should explore this use of coupled models more for several natural systems. Using numerical experiments, the ability of stochastic-dynamical models with assumed erroneous parameters to reproduce the nominal simulations should be studied. In addition, by introducing temporal dependence in the parameters a1 and a2 of (3), the effect of errors in storm motion velocity on flow predictions may be studied. In both these last two cases, radar-rainfall measurements may be simulated with given error statistics to serve as the rainfall measurement for the numerical experiment. Finally, the formulation of spatially-distributed coupled rainfall-flow models remains an important challenge. The spatial scale of soil water variability is a critical parameter in such a case (e.g., Guetter and Georgakakos, 1995).
Adler, R.F., and R.A. Mack, 1986: "Thunderstorm Cloud Top Dynamics as Inferred from satellite Observations and a Cloud Top Parcel Model," Journal of the Atmospheric Sciences, 43(18), 1945-1960.
Bae, D.H., Georgakakos, K.P., and Nanda, S.K., 1995: "Operational Forecasting with Real-Time Databases," in ASCE J. Hydraulics Division, 121(1), 49-60.
Brutsaert, W., 1991: "The Formulation of Evaporation from Land Surfaces," Recent Advances in the Modeling of Hydrologic Systems, eds. D.S. Bowles and P.E. O'Connell, D. Reidel Publishing Company, Boston, Mass., 67-84.
Georgakakos, K.P. and R.L. Bras, 1982: "Real-Time, Statistically Linearized, Adaptive Flood Routing," Water Resources Research, Vol. 18, No. 3, 513-524.
Georgakakos, K.P. and R.L. Bras, 1984a: "A Hydrologically Useful Station Precipitation Model, 1 Formulation," Water Resources Research, Vol. 20, No. 11, 1585-1596.
Georgakakos, K.P. and R.L. Bras, 1984b: "A Hydrologically Useful Station Precipitation Model, 2 Case Studies," Water Resources Research, Vol. 20, No. 11, 1597-1610.
Georgakakos, K.P. and Foufoula-Georgiou, E., 1991: "Real Time Coupling of Hydrological With Meteorological Models for Flood Forecasting", in Recent Advances in the Modeling of Hydrologic Systems, ed. Bowles, D.S. and O'Connell, P.E., D. Reidel Publishing Company, Boston, Mass., 169-184.
Georgakakos, K.P., 1986a: "A Generalized Stochastic Hydrometeorological Model for Flood and Flash-Flood Forecasting, 1 Formulation," Water Resources Research, Vol. 22, No. 13, 2083-2095.
Georgakakos, K.P., 1986b "A Generalized Stochastic Hydrometeorological Model for Flood and Flash-Flood Forecasting, 2 Case Studies," Water Resources Research, Vol. 22, No. 13, 2096-2106.
Georgakakos, K.P., 1987: "Real-Time Flash Flood Prediction," Journal of Geophysical Research, Vol. 92, No. D8, 9615-9629.
Guetter, A.K., and K.P. Georgakakos, "Large Scale Properties of Soil Water Variability," J. Geophysical Research, 1995 (accepted, available from the authors).
Henderson, F.M., 1966: Open Channel Flow, MacMillan Publishing Co., New York.
Hense, A., Krahe, P., and H. Flohn, 1988: "Recent Fluctuations of Tropospheric Temperature and Water Vapor Content in the Tropics," Meteorol. Atmos. Phys., 38, 215-227.
Hogg, D.C., Guiraud, F.O., Snider, J.B., Decker, M.T., and E.R. Westwater, 1983: "A Steerable Dual-Channel Microwave Radiometer for Measurement of Water Vapor and Liquid in the Troposphere," Journal of Applied Meteorology, 22(5), 789-806.
Mintz, Y., and G.K. Walker, 1993: "Global Fields of Soil Moisture and Land Surface Evapotranspiration Derived from Observed Precipitation and Surface Air Temperature," J. Applied Meteorology, 32, 1305-1334.
Ogura, Y., 1986: "Tropical Convection," Thunderstorm Morphology and Dynamics, Ed. E. Kessler, U. of Oklahoma Press, Norman, Oklahoma, 153-186.
Peck, E.L., 1976: "Catchment Modeling and Initial Parameter Estimation for the National Weather Service River Forecast System," NWS HYDRO-31, NWS/NOAA, Silver Spring, Maryland.
Pruppacher, H.R., and J.D. Klett, 1980: Microphysics of Clouds and Precipitation, D. Reidel Publishing Company, Boston, Massachusetts (second printing) 714 pp.
Randall, D.A., and J. Wang, 1992: "The Moist Available Energy of a Conditionally Unstable Atmosphere," Journal of the Atmospheric Sciences, 49(3), 240-255.
Wang, J.R., Melfi, S.H., Racette, P., Whitemen, D.N., Chang, L.A., Ferrare, R.A., Evans, K.D., and E.J. Schmidlin, 1995: "Simultaneous Measurements of Atmospheric water Vapor with MIR, Raman Lidar, and Rawinsondes," J. Applied Meteorology, 34, 1595-1607.
Webster, P.J., and R. Lukas, 1992: "TOGA-COARE: The Coupled Ocean - Atmosphere Response Experiment," Bulletin of the AMS, 73(9), 1377-1416.
Zawadzki, I., Torlaschi, E., and R. Sauvageau, 1981: "The Relationship Between Mesoscale thermodynamic variables and Convective Precipitation," Journal of the Atmospheric Sciences, 38, 1535-1540.