- Email: [email protected]
Simulation study of the effects of model uncertainty in variational assimilation of radar data on rainfall forecasting
Journal of Hydrology 239 (2000) 85–96 www.elsevier.com/locate/jhydrol
Simulation study of the effects of model uncertainty in variational assimilation of radar data on rainfall forecasting M. Grecu, W.F. Krajewski* Iowa Institute of Hydraulic Research, The University of Iowa, Iowa City, IA 52242, USA Received 5 October 1999; revised 13 September 2000; accepted 21 September 2000
Abstract Recent developments in data assimilation techniques make use of cloud models initialized with radar data an attractive alternative for real-time quantitative precipitation forecasting (QPF). Before such approaches are used operationally, there are a number of aspects that need to be addressed to understand better the benefits and drawbacks of the practical applications of cloud models in QPF. One aspect is the effect of various sources of uncertainty on the forecasting performance. Data assimilation formulations based on variational techniques allow accounting for uncertainty in observations but have no efficient mechanism of accounting for uncertainty in the model on which they are based. To investigate the issue, a simulation-based Monte Carlo methodology, suitable for the analysis of complex nonlinear models, is used. A one-dimensional stochasticdynamic cloud model, derived by considering stochastic terms in a physically based cloud model is used to simulate rainfall and radar reflectivity data. A deterministic version of the model is then initialized by a variational assimilation technique and used for forecasting. The differences between the forecasts and actual realizations of the stochastic cloud model are statistically analyzed to assess the effect of cloud model uncertainty on forecasting. In this paper this methodology is also used to study additional effects of other types of uncertainty, such as those in radar observations and in the description of rain drop size distribution, for a more complete understanding of the impact of uncertainties on rainfall forecasting. Based on the scenarios investigated in the paper, conclusions and recommendations concerning the use of complex cloud models in real-world applications are made. 䉷 2000 Elsevier Science B.V. All rights reserved. Keywords: Model uncertainty; Variational assimilation; Rainfall forecasting
1. Introduction In this paper we investigate the effects of uncertainty in variational assimilation (VA) of radar reflectivity data on quantitative precipitation forecasting (QPF) derived from cloud models. VA has been demonstrated to be an effective technique to retrieve atmospheric microphysical variables based on * Corresponding author. Fax: ⫹1-319-335-5238. E-mail address: [email protected] (W.F. Krajewski).
indirect observations. For example, Sun and Crook (1997) used a VA technique that involved a cloud model and obtained, based on radar reflectivity and Doppler velocity observations, a complete description of the local atmosphere, including its velocity, pressure, temperature, and various species of cloud water. Other convincing examples may be found in Bao and Warner (1993); Verlinde and Cotton (1993) and Sun and Crook (1996). An attractive potential application of VA is in quantitative rainfall forecasting. A cloud model initialized using a VA procedure and radar observations may be
0022-1694/00/$ - see front matter 䉷 2000 Elsevier Science B.V. All rights reserved. PII: S0022-169 4(00)00356-5
86
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
used to make forecasts. The practical application of this methodology for rainfall forecasting encounters numerous difficulties. Aside from high mathematical complexity, a major problem with the VA formulations is that they are able to account for errors in the observations, but cannot efficiently deal with cloud model errors. Usually, the model is assumed to be perfect, which, given the complexity of the physical processes and our limited understanding of them, is never the case. Attempts to account for uncertainty in models have been made (Zupanski, 1997), but they are still in the development phase and a general strategy has not been formulated yet. At present, rigorous alternatives to VA, such as Kalman filterbased techniques, are infeasible because of the prohibitively large number of variables involved. This justifies our study of the effect of random and systematic errors in cloud models on VA of radar reflectivity data, and, consequently, on cloud-modelbased rainfall forecasting. The high degree of non-linearity of cloud models makes it virtually impossible to analytically quantify the effect of errors on assimilation and, subsequently, forecasting performance. In this study, we used the most reliable alternative: Monte Carlo methods. Accordingly, we conducted a large number of assimilation experiments that differ from each other by random terms. We used a stochastic cloud model to generate rainfall events. We then synthesized radar reflectivity data based on these events and carried out VA using a deterministic cloud model. The deterministic cloud model, initialized based on the VA, is then used for forecasting. To evaluate the forecasting performance of the model, we statistically analyzed the differences between the forecasts and the reference events generated by the stochastic cloud model. In Monte Carlo experiments, the bigger the number of experiments conducted, the more statistically significant the results are. We determined that for our case, a reasonable number of experiments would be on the order of 200 (we will explain the reasons in detail later). Performing a similar Monte Carlo experiment of data assimilation using a two- or threedimensional cloud model is not feasible. This is why we opted for a one-dimensional (1D) cloud model, which is significantly less computer-intensive to use. Although the effect of observational errors on VA of radar data has been investigated (Sun and Crook,
1997), the effect of model errors has not yet been documented in the literature. The paper is organized as follows. We begin with a brief description of the 1D cloud model. Then, we describe the assimilation procedure and the Monte Carlo experiment in Sections 3 and 4. We follow with the results in Section 5, and close with a discussion and conclusions in Section 6. 2. 1D cloud model description The 1D cloud model we used is very similar to the one formulated by Ogura and Takahashi (1971). Since then, many other versions have been developed. Among the most recent are those of Ferrier and Houze (1989) and Yau (1980). These versions differ little in terms of momentum and transport equations; only in the microphysics schemes the differences are more significant. In fact, the study of different microphysics schemes was one of the main reasons 1D models were developed, when two (2-) and threedimensional (3D) versions were at hand. But the microphysics is not the point of interest in our study, and, consequently, we used the basic formulation of Ogura and Takahashi (1971). In the model, the momentum, mass, and energy conservation equations are integrated over the cloud cross-section, which is assumed to be circular. This integration results in average values of the cloud variables, but it also introduces additional terms similar to the Reynolds stresses in turbulent flows. These terms need to be modeled to obtain a determined system of equations. The details can be found in Ogura and Takahashi (1971). According to their results, the momentum equation may be written as 2w 2w 2a2 2 ~ u~ ⫹ B ⫺ gqH ⫹w ⫺ w兩w兩 ⫹ w ⫺ w 2t 2z a a 1
where value w is the mean vertical velocity, z is the vertical coordinate, t is time, a 2 is a proportionality constant which is set, based on experiments with air jets, to 0.1, u is the radial velocity, B is the buoyancy (Houze, 1993), g is the gravity acceleration and qH is the hydrometeor mixing ratio. Any variable with a tilde, e.g. A˜, is
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
defined as the circumference mean as follows: 1 Z2p A fdf A~ 2p 0
at r a;
2
where variable r is the distance from the z-axis, which is the symmetry axis of the cloud, and a is the cloud radius. The radial velocity u˜ is determined from the continuity equation 2 1 2 ÿ u~ ⫹ r w 0 r 2z a
3
87
formulation of Ogura and Takahashi (1971) and did not consider the pressure term. This led to simplified numerical calculations, and a less computer-intensive VA formulation. The thermodynamic equation is similar to that of momentum 2T 2T g 2 a2 T0 ⫺ T兩w兩 ⫹w ⫺ ⫹ a 2t 2z Cp 2 ~ u~ ⫹ ST ⫹ T ⫺ T a
(5)
where rÅ is a mean value of air density, which depends on the height only. The other circumference-averaged values are specified as a function of the sign of u˜. For any generic variable A, the following formula is valid ( A~ A0 if u~ ⬍ 0 4 A~ A if u~ ⬎ 0
where ST is a thermal source or sink due to evaporation or condensation of water. The conservation equations for water vapor, cloud and precipitation water take into account four processes. These are: condensation (S), autoconversion of cloud water into rainwater (A), collection of cloud water (C) and precipitation evaporation (E). Thus,
where A is the variable averaged over the cloud crosssection and A0 is its corresponding value for the surrounding environment. If u~ ⬍ 0; the surrounding air is entrained into the cloud. Otherwise, there is detrainment of cloud air towards the environment. The first term on the right-hand side of Eq. (1) represents lateral eddy exchange, and the second term models the dynamic interaction between the cloud and the environment. Unlike the other terms, the first two are specific to 1D cloud modeling. They are not needed in higher-dimension models where the cloud environment interaction is treated implicitly. The 1D model we used in this study may be regarded as “one and a half”-dimensional because of its ability to consider lateral mixing and entrainment/ detrainment. Note that Eq. (1) does not contain a pressure term. This is because there is no momentum equation for the radial velocity u. The radial velocity is determined from the continuity equation, which is the constraint that provides the pressure in 2- or 3D models. Consequently, it is impossible to derive an exact Poisson equation starting from the continuity and only one momentum equation. However, approximate Poisson equations may be derived from the exact formulation of the 2D axially symmetric case by making various assumptions (Holton, 1973; Yau, 1980; Ferrier and Houze, 1989). In our study, we followed the original
2ql 2q 2a 2 2 ⫹w l ql0 ⫺ ql 兩w兩 ⫹ ql ⫺ q~l u~ ⫹ S 2t 2z a a ⫺A⫺C
(6)
2qr 2q 1 2 2 a2 ⫹w r qr0 ⫺ qr 兩w兩 r Vr qr ⫹ 2t 2z r 2z a ⫹
2 q ⫺ q~r u~ ⫹ A ⫹ C ⫺ E a r
7
2qv 2q 2 a2 ⫹w v qv0 ⫺ qv 兩w兩 2t 2z a 2 ⫹ qv ⫺ q~ v u~ ⫺ S ⫹ E a
(8)
where ql, qr, qv, are the mixing ratios for cloud water, rain and water vapor, respectively. We denote raindrop terminal velocity Vr. It may be determined as Vr 21:9 r qr 0:125 with Vr in m/s, rÅ in g/cm 3 and qr in g/g. We model the terms A, C and E according to Kessler’s parameterization (Kesller, 1969): A max 0; ac ql ⫺ qlcrit ; where ac 0:001 s⫺1 and ⫺1 qlcrit 1:5 g=kg; C gql q0:95 r ; where g 0:002 s ; and E max 0; ⫺b qv ⫺ qvs ; where b 0:001 and qvs is the water vapor mixing ratio at the cloud temperature and pressure. Condensation S is determined such
88
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
that the physical constraint qv ⱕ qvs is always satisfied. This is achieved through a moist adjustment process (Sun and Crook, 1997). Based on the S and E, ST may be determined as ST
L S ⫺ E: CP
9
Eqs. (1)–(8) are completely deterministic. For a given set of initial conditions, they yield uniquely determined results. Random, high-intensity fluctuations in cloud variables cannot be predicted by these equations. On the other hand, actual precipitation is characterized by such fluctuations. From the practical point of view, it is useful to know how the random fluctuations sensed by the observing systems influence the data assimilation process performed by means of a deterministic model. To accomplish this, we perturbed the 1D model with random terms to simulate precipitation events and synthesized radar observations. Then, we performed the VA using the unperturbed cloud model, and used that (deterministic) model, initialized by the VA, for forecasting. We conducted a large number of such experiments (runs), and summarized the results statistically. Below we describe this methodology in detail. 3. VA of radar data in cloud models Although rigorous presentations of VA formulations may be found in several papers (see for example Tremolet and Le Dimet, 1997; Chao and Chang, 1992), we describe VA briefly here to identify more easily the uncertainty-related issues. Suppose that the cloud model is written in a generic, discrete form as X i⫹1 F X i
10
i
where X is the state vector variable (it contains the complete information to fully describe the model physical state in space and time, i.e. velocity, temperature and the three species of water fields), and F is the model that describes the system evolution. Superscript i indicates the time. Given several sets of measurements at different times ti Y i H X i ⫹ ei
i 僆 {1; 2; …; N}
11
where Y i are the measured variables, H is the observa-
tion function which relates the state variable to the measurements and ei are Gaussian errors with known characteristics, one may initialize the cloud model at time t1, such that the model output obtained by applying Eq. (10) is repeatedly in agreement with observation (11). Thus, an objective function may be defined as follows: J
N X
Y i ⫺ H X i T Wi⫺1 Y i ⫺ H X i
12
i1
where Wi is a square weight matrix. The problem may be reduced to finding the variable X1 that minimizes Eq. (12) with constraint (10). Objective function (12) may be expressed exclusively as a function of X1, and its gradient with respect to X1 becomes 7J ⫺2
N X
Y i ⫺ H X i T Wi⫺1 7H i 7F i⫺1 …7F 1
i1
13 where, for simplicity of notation 7F X 7F and 7H X i 7H i : Then, a gradient-based optimization technique may be used to determine the initial condition X1 that minimizes Eq. (12). The actual number of observations may be significantly less than N. In other words, Eq. (10) must be applied N ⫺ 1 times from t1 to tN, but observations are not necessarily available at each intermediate time ti, where t1 ⬍ ti ⬍ tN : Then, formula (13) is still valid if the corresponding weight matrices are set to infinity. The expression on the right-hand side of Eq. (13) may be evaluated efficiently only by performing the algebraic operations starting from the left-most factor to the right. To do so, one must develop a special model, called the adjoint, from the forward model (10). The adjoint model may be considered as the transposed of the original model’s gradient but evolving reversely from the final time to the initial time (Tremolet and Le Dimet, 1997). The task of deriving it is quite difficult. In this paper, we used the automatic code of Gering and Kaminski (1998) for deriving the adjoint model from the forward one. Once evaluated, gradient (13) may be used by an optimization procedure to determine the initial condition X1. The number of variables involved in the optimization is quite large, even for a 1D cloud model. Consequently, we used a special optimization procedure, namely LMBFGS of Liu and Nocedal i
i
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
4. Monte Carlo experiment setup The 1D cloud model formulated in Eq. (10) may be transformed into a stochastic model of the type X i⫹1 F X i ⫹ Gi win
14
where G is a coefficient matrix, and wn a sequence of model errors. The stochastic terms introduced in the deterministic model account for the unresolved interactions of the cloud with both large and small scales. Indirect evidence, i.e. the variability observed in rainfall fields (Georgakakos and Krajewski, 1996), suggests that these unsolved interactions may be described by stochastic terms. We consider perturbations caused by unresolved interactions only in the vertical velocity and temperature deviation. Stochastic effects in the other variables are induced by modi~ u=a ~ at random. fying the entrainment terms 2 A ⫺ A This is equivalent to considering that there is a degree of randomness in the cloud/environment interaction. Specifically, we treat the stochastic component as follows. First, we integrate cloud model (10) one
step forward with the only modification being that any entrainment term multiplier becomes 2 ⫹ wn =a (as opposed to 2=a where wn is a normal random variable with mean zero and standard deviation one. Then, the temperature and velocity are perturbed, Ti⫹1 Ti⫹1 ⫹ 0:5wnt and Wi⫹1 Wi⫹1 ⫹ 0:5Wnw ; where wnw, wnt are random variables with mean zero and standard deviation one. We consider the perturbations to be independent. This formulation of the stochastic terms in our model is similar to those of Georgakakos and Krajewski (1996); Tsintikidis and Georgakakos (1999) and Tao and Soong (1987) who provided indirect evidence pertaining to the suitability of such terms. In our numerical experiments, we specified the environmental conditions for the stochastic 1D cloud model by the sounding shown in Fig. 1. The model was initialized by a thermal perturbation that extends to an altitude of 4 km and consists of 3 K temperature excess. We set all other variables to the corresponding values of the environment. We prescribed the boundary conditions as follows. The vertical velocity was zero at both boundaries, and the rain content was zero at the top, and interpolated from above at the bottom. We set the cloud radius to 5 km, and the computational domain from 0 to 17 km, with the vertical grid size of 500 m and the time step of 5 s. With these initial and boundary conditions, we ran the stochastic cloud model 200 times for different random terms, i.e. different initializations of the random number generator. The reason for choosing 100 Temperature Pressure (mbar)
(1989), which is efficient for a large number of variables. Even when model (10) is perfect, if the observations do not contain complete information, estimate X1 may be subject to errors. In that case, there are multiple optimal solutions for objective function (12). Thus, the optimization procedure finds an arbitrary solution that may be far from the true solution (global optimum). If model (10) itself is inaccurate (aside from incomplete information in the observations), X1 will contain additional errors. In this paper, we investigate how the errors induced in X1 by the inaccuracy of the model influence the forecasting performance. A Kalman filter-based formulation is less sensitive to model errors because objective function (12) is modified to account for an inaccurate model, and model (10) is not used explicitly in gradient calculation (13). However, a Kalman filter-approach is infeasible for large problems because it requires operations with large dimension matrices. The VA framework presented here is based mostly on operations with vectors that are affordable on current computer systems, and this justifies its use versus Kalman filters.
89
200 300 400 500 600 700 800 900 -90.0
Dew point
-60.0 -30.0 0.0 o Temperature ( C)
30.0
Fig. 1. Sounding data. The dashed line represents the temperature of an air parcel raised from the surface. The other two lines are for the temperature and dew point.
90
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
16 14
qr (g/kg)
12 Z (km)
qr (g/kg)
1 2 3 4
10
1 2 3 4
8 6 4 2 16 14
Rain rate (mm/h)
Z (km)
12 10
Rain rate (mm/h)
10 20 30 40
10 20 30 40
8 6 4 2
Ground rain rate (mm/h)
40
30
20
10
0 0
20
40 60 Time (minutes)
0
20
40 60 Time (minutes)
Fig. 2. Examples of the stochastic model output. Shown are cloud variables as a function of time. Top panels are contours of rainwater content (g/kg), middle panels are contours of rain rate (mm/h), and the bottom panels are for surface rain rate. The left column represents the mean, and the right column the standard deviation over 200 realizations.
200 as the number of repetitions that yields statistically significant results comes from the central limit theorem (Hogg and Tannis, 1997). According to this theorem, the confidence in estimating the mean of a normal random variable increases proportionally with the square root of the size of sample set used in estimation. A sample size equal to 200 leads to a p confidence interval proportional to 1= 200; which seems to be a reasonable tradeoff between accuracy and computational effort.
In Fig. 2, we summarize the mean over the 200 realizations of cloud variables as functions of time. The corresponding rain variables obtained by running the cloud model with the same initial and boundary conditions but without random perturbations, are shown as a reference in Fig. 3. The differences between the rain related variables yielded by the deterministic model (9) and the means of the same variables obtained from the stochastic model are also shown in Fig. 3. We may conclude based on
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
91
Fig. 3. Deterministic model output. Shown variables are organized as in Fig. 2. The left column is now a single run output, and the right column is the difference between the deterministic and the mean of the stochastic model.
Figs. 2 and 3 that, as anticipated, the average output of the stochastic model is different from that of the deterministic one. It is often assumed that the deterministic model may be considered to represent the modeled reality in some “average sense”. However, applying the stochastic model and taking the average of the results is not the same as taking the mean model and applying it. This is because of the model’s nonlinearity. Considering the mean model instead of the average output of the stochastic one
leads to systematic errors. In our case, the systematic errors induced by the actual deterministic model are relatively small compared to those caused by the stochastic variability. To examine the effect of assimilating radar data in view of such variability, we conducted a Monte Carlo experiment. The steps involved were as follows: (1) generate a sample using the stochastic cloud model;
92
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
Table 1 Uncertainty parameters in assimilation experiments Case
Radar data error
Deterministic model
Drop size distribution
C1 C2 C3 C4 C5
No error No error N(0,9) N(0,9) No error
Exact Inexact Exact Inexact Exact
Known Known Known Known Unknown
(2) synthesize three consecutive, 4.16 min apart (50 times the time step which is 5 s) radar observations; (3) assimilate the radar observations in the deterministic cloud model; (4) forecast the cloud variables using the deterministic cloud model; (5) analyze the forecast performance; (6) repeat steps (1) through (5) 200 times; (7) characterize statistically the results of step (5). We ran the stochastic cloud model for 75 min in step (1). After t 75 min the convection was quite weak, and the rainfall rate was small relative to its peak. Consequently, we did not investigate the process beyond that moment. The assimilation interval started at 16.67 min and ended at 25.00 min. There were three sets of radar observations in the assimilation interval, at t1 16:67; t2 20:83 and
Rain rate error (mm/h)
8
4
0
-4
-8
0
20
40 Time (min tes)
60
Fig. 4. Case C1. The thick line represents the mean error, and the shaded region the error within one standard deviation bounds.
t3 25:00 min: We chose time t 16:67 min as the starting point of the assimilation interval because it provided significant information early enough for the forecasts to have utility. The objective function contained only one term, namely J
3 X
⫺1 i Z i ⫺ Z i;obs T WZ;i Z ⫺ Z i;obs ;
15
i1
where Z i was the vector containing the radar reflectivity factors in all grid points at time ti. Superscript “obs” indicates the “observed” values. The radar observations were synthesized from the rain content assuming a gamma drop size distribution (Ulbrich, 1983), with N0 8 × 104 m3 cm⫺1 and m 1:15 and Rayleigh scattering (Doviak and Zrnic, 1993). In a particular case, which will be specified later, we considered the Marshall and Palmer’s drop-size distribution to investigate the effect of uncertainty model relative to that of other types of uncertainty. It was assumed that the radar sampling volumes were the same as the model’s computational cells and simultaneous with model time steps. In other words, no uncertainty caused by the radar spatial and temporal sampling was considered. The Doppler velocity information could not be projected in the cylindrical coordinate system of the 1D cloud model. Therefore, we did not consider it in the objective function, although in 2- or 3D formulations it is usually considered. We assumed the observational errors were not correlated, and consequently, WZ,i was set to the identity matrix. 5. Results We applied the Monte Carlo experimental setup described in the previous section in a series of cases.
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
Ground rain rate (mm/h)
40 Exact Inexact 30
20
10
0 0
20
40 Time (minutes)
60
Fig. 5. Ground rain rate of the exact and inexact models.
They differ by the degree of noise in the radarsynthesized measurements and the accuracy of the model in assimilation. That is, the radar-synthesized observations are corrupted by noise, and the deterministic model is modified to yield results significantly different from its original output. These cases provide means to investigate the effect of radar errors and model uncertainty on the assimilation and forecasting performance. Another major uncertainty factor that affects the assimilation performance is the lack of information about the rain drop size distribution. We also studied its consequences. Table 1 defines the cases investigated. Normal noise with zero mean and standard deviation 3 dB is expressed as N(0,9). In case C1, the radar observations are considered error free. The stochastic model is used to produce different realizations. Radar observations corresponding to these realizations are produced, and the deterministic model is used for assimilation and forecasting. The mean and standard deviation of the differences between the forecasted and “true” cloud variables (velocity, temperature, water vapor, cloud and rain fields) are then calculated for each realization as a function of time. The results for the rain rate are given in Fig. 4. We can see that the relative errors are small, and there is a fairly good agreement between the forecasted and true rain rate. The wavy appearance of the error might be an indication that the assimilation problem is not well posed, i.e. does not have unique solution, stable to small perturbations. However, the rain rate forecasted globally is quite
93
accurate, although it is underestimated in some periods and overestimated in others. The standard deviation of the error is smaller than the variance of the true rain. This indicates that the inability of a deterministic model to follow the complex cloud environment interactions (mimicked here by the stochastic terms used to generate the reference rain), rather than the assimilation procedure, is responsible for the errors. In case C2, the entrainment terms in the deterministic model are artificially increased 1.4 times. The ground rainfall rates of the initial and altered deterministic cloud models are given in Fig. 5. They are significantly different. The altered model produces about 40% less rainfall. The reason for making this arbitrary modification is to investigate the effect of using an inexact model in assimilation and forecasting on QPF procedure. Even the most complex cloud models may be quite poor in replicating real phenomena. Therefore, it is important to characterize the performance degradation introduced by the model inaccuracy. The standard deviation and mean of the forecasting error, defined as the difference between the forecast and truth, are given in Fig. 6. It is interesting to note that although there is bias, this is significantly smaller (about four times) than that of the deterministic model itself. This is an indication that the assimilation procedure compensates for the inaccuracy of its model. In other words, the optimization procedure finds the initial conditions that make the dynamics of the inaccurate model match with that of the rainfall process. The inaccuracy of the deterministic model can be rigorously accounted for by introducing a set of parameters that influences its dynamics. In the current problem, what caused the negative bias is known, but in real applications only the effects may be observed. In those situations, it is possible to introduce heuristically parameters that make the model more flexible. The parameters may be determined optimally along with model’s initial conditions or state variables (Rajaram and Georgakakos, 1989). This is generally an efficient strategy to mitigate the model inaccuracy. However, in the present study, we did not use it because only one assimilation interval (three sets of observations) was considered for each realization and the requirements of the formalism of Rajaram and Georgakakos (1989) were not met.
94
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
Rain rate error (mm/h)
12 C2
C3
C4
C5
6
0
-6
Rain rate error (mm/h)
-12 12
6
0
-6
-12
0
20
40 60 Time (minutes)
0
20
40 60 Time (minutes)
Fig. 6. Errors of the forecasted rain rate for cases C2, C3, C4 and C5. The thick line represents the mean error, and the shaded region the error within one standard deviation bounds.
Despite this, the assimilation yielded satisfactory results even with an inexact model. Cases C3 and C4 are similar to C1 and C2, but include noise in observations. Gaussian noise with zero mean and 3 dB standard deviation is added to every simulated reflectivity value in the computational domain. The results, given in Fig. 6, are comparable to those of cases C1 and C2, although an increase in errors is noticeable both in means and standard deviations. In case C5, two different gamma drop size distributions are used to express the reflectivity as a function of rainwater content. The reflectivity observations are synthesized based on N0 8 × 104 m3 cm⫺1 and m 1:15 and calculated using the same N0 but m 0:0: The first formula gives a slower growth of reflectivity with rain than does the second one. This implies that for the same reflectivity values, the rain rate estimated
by the second formula is smaller than that of the first one. Therefore, it is expected that the rain content is underestimated in assimilation, and consequently, in forecasting. Fig. 6 confirms this. Note that the error is significantly larger than that of C1, and similar to that of C4. This leads to the conclusion that the uncertainty in raindrop size distribution has an impact as big as that of using a significantly different model in assimilation.
6. Summary and conclusions In this paper, the effect of uncertainty on variational assimilation of radar data in cloud models is investigated. This is accomplished by using Monte Carlo simulations. A 1D cloud model with appropriate modifications to account for the complex interactions
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
between the cloud and environment is used. The 1D model is altered with random terms and run to produce particular realizations of precipitation events. These events are used as references, to synthesize radar observations of reflectivity, and to evaluate the assimilation and forecasting performance. The assimilation is performed in conditions of several kinds of uncertainty, namely, systematic errors in the assimilation model, observation errors and lack of knowledge about rain drop size distribution. The results show that assimilation based the 1D model is quite robust, i.e. it is mathematically well conditioned (it has a solution in conditions of several kinds of uncertainty). We found that the assimilation is able to compensate for systematic errors in the model. Even significant errors in radar observations (normally distributed with zero mean and 3 dB standard deviation) degrade the performance only slightly. However, the assimilation is quite sensitive to the uncertainty in reflectivity introduced by the rain drop size distribution. Based on these results, we conclude that VA of radar data in cloud models may be effective in rainfall forecasting, despite uncertainties in the model. In addition to errors caused by the model uncertainty, other errors may occur in 2- or 3D cloud models, either from the ill-posed character of the problem (multiple solutions for the initial conditions), or from the optimization procedure’s inability to find the best initial conditions (Sun and Crook, 1997). The results in this paper suggest that it may be beneficial for forecasting to reduce the complexity of the models to the level where a unique, optimal estimation of the initial conditions may be obtained. As we demonstrated here, the VA, when mathematically well defined, compensates for the uncertainty of the model used. This issue needs to be further investigated.
Acknowledgements This study was supported by NASA Grant NAG81425 under the auspices of the US Weather Research Program. The radar data were provided under the Cooperative Agreement between the Office of Hydrology of the National Weather Service and the Iowa Institute of Hydraulic Research (NA47WH0495). We thank Drs Konstantine Georgakakos and Mark French for useful discussions.
95
References Bao, J.W., Warner, T.T., 1993. Treatment of on/off switches in the adjoint method: FDDA experiments with a simple model. Tellus 45A, 525–538. Chao, W.C., Chang, L.P., 1992. Development of a four-dimensional variational analysis system using the adjoint method at GLA. 1. Dynamics. Monthly Weather Review 20 (8), 1661–1673. Doviak, R.J., Zrnic, D., 1993. Doppler Radar and Weather Observations. Academic Press, San Diego, CA (458 pp.). Ferrier, B.S., Houze, Jr., R, A., 1989. One-dimensional timedependent modeling of GATE cumulonimbus convection. Journal of the Atmospheric Sciences 46, 330–352. Georgakakos, K.P., Krajewski, W.K., 1996. Statistical-microphysical causes of rainfall variability in the tropics. Journal of Geophysical Research 101 (26), 165–180. Gering, R., Kaminski, T., 1998. Recipes for adjoint code construction. ACM Transactions on Mathematical Software 24 (4), 437– 474. Hogg, R.V., Tannis, E.A., 1997. Probability and Statistical Inference. Prentice Hall, Upper Saddle River, NJ, p. 782. Holton, J.R., 1973. A one-dimensional cumulus model including pressure perturbations. Monthly Weather Review 101, 201– 205. Houze Jr., R.A., 1993. Cloud Dynamics. Academic Press, New York. Kesller, E., 1969. On the distribution of water substance in atmospheric circulations. Meteorology Monographs, Vol. 10. American Meteorological Society, Boston, MA (84 pp.). Liu, D.C., Nocedal, J., 1989. On the limited memory BFGS method for large-scale optimization. Mathematical Programming 45, 503–528. Ogura, Y., Takahashi, T., 1971. Numerical simulation of the life cycle of a thunderstorm cell. Monthly Weather Review 99 (12), 895–910. Rajaram, H., Georgakakos, K.P., 1989. Recursive parameter estimation of hydrologic models. Water Resources Research 25 (2), 281–294. Sun, J.-Z., Crook, N.A., 1996. Comparison of thermodynamic retrieval by the adjoint method with the traditional retrieval method. Monthly Weather Review 124 (2), 308–324. Sun, J.-Z., Crook, N.A., 1997. Dynamical and microphysical retrieval from Doppler radar observations using a cloud model and its adjoint. 1. Model development and simulated data experiments. Journal of the Atmospheric Sciences 54 (12), 1642–1661. Tao, W.-K., Soong, S.-T., 1987. Statistical properties of a cloud ensemble: a numerical study. Journal of the Atmospheric Sciences 44, 3175–3187. Tsintikidis, D., Georgakakos, K.P., 1999. Microphysical and largescale dependencies of temporal rainfall variability over a tropical ocean. Journal of the Atmospheric Sciences 56 (5), 724–748. Tremolet, Y., Le Dimet, F.-X., 1997. Parallel algorithms for variational data assimilation and coupling models. Parallel Computing 22 (5), 657–674. Ulbrich, C.W., 1983. Natural variations in the analytical form of raindrop size distribution. Journal of Climate and Applied Meteorology 22, 1764–1775.
96
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
Verlinde, J., Cotton, W.R., 1993. Fitting microphysical of nonsteady convective clouds to a numerical model: an application of the adjoint technique of data assimilation to a kinematic model. Monthly Weather Review 121, 2776–2793. Yau, M.K., 1980. A two-cylinder model of cumulus cells and its
application in computing cumulus transports. Journal of the Atmospheric Sciences 37 (11), 2470–2485. Zupanski, D., 1997. A general weak constraint applicable to operational 4DVAR data assimilation systems. Monthly Weather Review 125 (9), 2274–2292.
Simulation study of the effects of model uncertainty in variational assimilation of radar data on rainfall forecasting M. Grecu, W.F. Krajewski* Iowa Institute of Hydraulic Research, The University of Iowa, Iowa City, IA 52242, USA Received 5 October 1999; revised 13 September 2000; accepted 21 September 2000
Abstract Recent developments in data assimilation techniques make use of cloud models initialized with radar data an attractive alternative for real-time quantitative precipitation forecasting (QPF). Before such approaches are used operationally, there are a number of aspects that need to be addressed to understand better the benefits and drawbacks of the practical applications of cloud models in QPF. One aspect is the effect of various sources of uncertainty on the forecasting performance. Data assimilation formulations based on variational techniques allow accounting for uncertainty in observations but have no efficient mechanism of accounting for uncertainty in the model on which they are based. To investigate the issue, a simulation-based Monte Carlo methodology, suitable for the analysis of complex nonlinear models, is used. A one-dimensional stochasticdynamic cloud model, derived by considering stochastic terms in a physically based cloud model is used to simulate rainfall and radar reflectivity data. A deterministic version of the model is then initialized by a variational assimilation technique and used for forecasting. The differences between the forecasts and actual realizations of the stochastic cloud model are statistically analyzed to assess the effect of cloud model uncertainty on forecasting. In this paper this methodology is also used to study additional effects of other types of uncertainty, such as those in radar observations and in the description of rain drop size distribution, for a more complete understanding of the impact of uncertainties on rainfall forecasting. Based on the scenarios investigated in the paper, conclusions and recommendations concerning the use of complex cloud models in real-world applications are made. 䉷 2000 Elsevier Science B.V. All rights reserved. Keywords: Model uncertainty; Variational assimilation; Rainfall forecasting
1. Introduction In this paper we investigate the effects of uncertainty in variational assimilation (VA) of radar reflectivity data on quantitative precipitation forecasting (QPF) derived from cloud models. VA has been demonstrated to be an effective technique to retrieve atmospheric microphysical variables based on * Corresponding author. Fax: ⫹1-319-335-5238. E-mail address: [email protected] (W.F. Krajewski).
indirect observations. For example, Sun and Crook (1997) used a VA technique that involved a cloud model and obtained, based on radar reflectivity and Doppler velocity observations, a complete description of the local atmosphere, including its velocity, pressure, temperature, and various species of cloud water. Other convincing examples may be found in Bao and Warner (1993); Verlinde and Cotton (1993) and Sun and Crook (1996). An attractive potential application of VA is in quantitative rainfall forecasting. A cloud model initialized using a VA procedure and radar observations may be
0022-1694/00/$ - see front matter 䉷 2000 Elsevier Science B.V. All rights reserved. PII: S0022-169 4(00)00356-5
86
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
used to make forecasts. The practical application of this methodology for rainfall forecasting encounters numerous difficulties. Aside from high mathematical complexity, a major problem with the VA formulations is that they are able to account for errors in the observations, but cannot efficiently deal with cloud model errors. Usually, the model is assumed to be perfect, which, given the complexity of the physical processes and our limited understanding of them, is never the case. Attempts to account for uncertainty in models have been made (Zupanski, 1997), but they are still in the development phase and a general strategy has not been formulated yet. At present, rigorous alternatives to VA, such as Kalman filterbased techniques, are infeasible because of the prohibitively large number of variables involved. This justifies our study of the effect of random and systematic errors in cloud models on VA of radar reflectivity data, and, consequently, on cloud-modelbased rainfall forecasting. The high degree of non-linearity of cloud models makes it virtually impossible to analytically quantify the effect of errors on assimilation and, subsequently, forecasting performance. In this study, we used the most reliable alternative: Monte Carlo methods. Accordingly, we conducted a large number of assimilation experiments that differ from each other by random terms. We used a stochastic cloud model to generate rainfall events. We then synthesized radar reflectivity data based on these events and carried out VA using a deterministic cloud model. The deterministic cloud model, initialized based on the VA, is then used for forecasting. To evaluate the forecasting performance of the model, we statistically analyzed the differences between the forecasts and the reference events generated by the stochastic cloud model. In Monte Carlo experiments, the bigger the number of experiments conducted, the more statistically significant the results are. We determined that for our case, a reasonable number of experiments would be on the order of 200 (we will explain the reasons in detail later). Performing a similar Monte Carlo experiment of data assimilation using a two- or threedimensional cloud model is not feasible. This is why we opted for a one-dimensional (1D) cloud model, which is significantly less computer-intensive to use. Although the effect of observational errors on VA of radar data has been investigated (Sun and Crook,
1997), the effect of model errors has not yet been documented in the literature. The paper is organized as follows. We begin with a brief description of the 1D cloud model. Then, we describe the assimilation procedure and the Monte Carlo experiment in Sections 3 and 4. We follow with the results in Section 5, and close with a discussion and conclusions in Section 6. 2. 1D cloud model description The 1D cloud model we used is very similar to the one formulated by Ogura and Takahashi (1971). Since then, many other versions have been developed. Among the most recent are those of Ferrier and Houze (1989) and Yau (1980). These versions differ little in terms of momentum and transport equations; only in the microphysics schemes the differences are more significant. In fact, the study of different microphysics schemes was one of the main reasons 1D models were developed, when two (2-) and threedimensional (3D) versions were at hand. But the microphysics is not the point of interest in our study, and, consequently, we used the basic formulation of Ogura and Takahashi (1971). In the model, the momentum, mass, and energy conservation equations are integrated over the cloud cross-section, which is assumed to be circular. This integration results in average values of the cloud variables, but it also introduces additional terms similar to the Reynolds stresses in turbulent flows. These terms need to be modeled to obtain a determined system of equations. The details can be found in Ogura and Takahashi (1971). According to their results, the momentum equation may be written as 2w 2w 2a2 2 ~ u~ ⫹ B ⫺ gqH ⫹w ⫺ w兩w兩 ⫹ w ⫺ w 2t 2z a a 1
where value w is the mean vertical velocity, z is the vertical coordinate, t is time, a 2 is a proportionality constant which is set, based on experiments with air jets, to 0.1, u is the radial velocity, B is the buoyancy (Houze, 1993), g is the gravity acceleration and qH is the hydrometeor mixing ratio. Any variable with a tilde, e.g. A˜, is
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
defined as the circumference mean as follows: 1 Z2p A fdf A~ 2p 0
at r a;
2
where variable r is the distance from the z-axis, which is the symmetry axis of the cloud, and a is the cloud radius. The radial velocity u˜ is determined from the continuity equation 2 1 2 ÿ u~ ⫹ r w 0 r 2z a
3
87
formulation of Ogura and Takahashi (1971) and did not consider the pressure term. This led to simplified numerical calculations, and a less computer-intensive VA formulation. The thermodynamic equation is similar to that of momentum 2T 2T g 2 a2 T0 ⫺ T兩w兩 ⫹w ⫺ ⫹ a 2t 2z Cp 2 ~ u~ ⫹ ST ⫹ T ⫺ T a
(5)
where rÅ is a mean value of air density, which depends on the height only. The other circumference-averaged values are specified as a function of the sign of u˜. For any generic variable A, the following formula is valid ( A~ A0 if u~ ⬍ 0 4 A~ A if u~ ⬎ 0
where ST is a thermal source or sink due to evaporation or condensation of water. The conservation equations for water vapor, cloud and precipitation water take into account four processes. These are: condensation (S), autoconversion of cloud water into rainwater (A), collection of cloud water (C) and precipitation evaporation (E). Thus,
where A is the variable averaged over the cloud crosssection and A0 is its corresponding value for the surrounding environment. If u~ ⬍ 0; the surrounding air is entrained into the cloud. Otherwise, there is detrainment of cloud air towards the environment. The first term on the right-hand side of Eq. (1) represents lateral eddy exchange, and the second term models the dynamic interaction between the cloud and the environment. Unlike the other terms, the first two are specific to 1D cloud modeling. They are not needed in higher-dimension models where the cloud environment interaction is treated implicitly. The 1D model we used in this study may be regarded as “one and a half”-dimensional because of its ability to consider lateral mixing and entrainment/ detrainment. Note that Eq. (1) does not contain a pressure term. This is because there is no momentum equation for the radial velocity u. The radial velocity is determined from the continuity equation, which is the constraint that provides the pressure in 2- or 3D models. Consequently, it is impossible to derive an exact Poisson equation starting from the continuity and only one momentum equation. However, approximate Poisson equations may be derived from the exact formulation of the 2D axially symmetric case by making various assumptions (Holton, 1973; Yau, 1980; Ferrier and Houze, 1989). In our study, we followed the original
2ql 2q 2a 2 2 ⫹w l ql0 ⫺ ql 兩w兩 ⫹ ql ⫺ q~l u~ ⫹ S 2t 2z a a ⫺A⫺C
(6)
2qr 2q 1 2 2 a2 ⫹w r qr0 ⫺ qr 兩w兩 r Vr qr ⫹ 2t 2z r 2z a ⫹
2 q ⫺ q~r u~ ⫹ A ⫹ C ⫺ E a r
7
2qv 2q 2 a2 ⫹w v qv0 ⫺ qv 兩w兩 2t 2z a 2 ⫹ qv ⫺ q~ v u~ ⫺ S ⫹ E a
(8)
where ql, qr, qv, are the mixing ratios for cloud water, rain and water vapor, respectively. We denote raindrop terminal velocity Vr. It may be determined as Vr 21:9 r qr 0:125 with Vr in m/s, rÅ in g/cm 3 and qr in g/g. We model the terms A, C and E according to Kessler’s parameterization (Kesller, 1969): A max 0; ac ql ⫺ qlcrit ; where ac 0:001 s⫺1 and ⫺1 qlcrit 1:5 g=kg; C gql q0:95 r ; where g 0:002 s ; and E max 0; ⫺b qv ⫺ qvs ; where b 0:001 and qvs is the water vapor mixing ratio at the cloud temperature and pressure. Condensation S is determined such
88
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
that the physical constraint qv ⱕ qvs is always satisfied. This is achieved through a moist adjustment process (Sun and Crook, 1997). Based on the S and E, ST may be determined as ST
L S ⫺ E: CP
9
Eqs. (1)–(8) are completely deterministic. For a given set of initial conditions, they yield uniquely determined results. Random, high-intensity fluctuations in cloud variables cannot be predicted by these equations. On the other hand, actual precipitation is characterized by such fluctuations. From the practical point of view, it is useful to know how the random fluctuations sensed by the observing systems influence the data assimilation process performed by means of a deterministic model. To accomplish this, we perturbed the 1D model with random terms to simulate precipitation events and synthesized radar observations. Then, we performed the VA using the unperturbed cloud model, and used that (deterministic) model, initialized by the VA, for forecasting. We conducted a large number of such experiments (runs), and summarized the results statistically. Below we describe this methodology in detail. 3. VA of radar data in cloud models Although rigorous presentations of VA formulations may be found in several papers (see for example Tremolet and Le Dimet, 1997; Chao and Chang, 1992), we describe VA briefly here to identify more easily the uncertainty-related issues. Suppose that the cloud model is written in a generic, discrete form as X i⫹1 F X i
10
i
where X is the state vector variable (it contains the complete information to fully describe the model physical state in space and time, i.e. velocity, temperature and the three species of water fields), and F is the model that describes the system evolution. Superscript i indicates the time. Given several sets of measurements at different times ti Y i H X i ⫹ ei
i 僆 {1; 2; …; N}
11
where Y i are the measured variables, H is the observa-
tion function which relates the state variable to the measurements and ei are Gaussian errors with known characteristics, one may initialize the cloud model at time t1, such that the model output obtained by applying Eq. (10) is repeatedly in agreement with observation (11). Thus, an objective function may be defined as follows: J
N X
Y i ⫺ H X i T Wi⫺1 Y i ⫺ H X i
12
i1
where Wi is a square weight matrix. The problem may be reduced to finding the variable X1 that minimizes Eq. (12) with constraint (10). Objective function (12) may be expressed exclusively as a function of X1, and its gradient with respect to X1 becomes 7J ⫺2
N X
Y i ⫺ H X i T Wi⫺1 7H i 7F i⫺1 …7F 1
i1
13 where, for simplicity of notation 7F X 7F and 7H X i 7H i : Then, a gradient-based optimization technique may be used to determine the initial condition X1 that minimizes Eq. (12). The actual number of observations may be significantly less than N. In other words, Eq. (10) must be applied N ⫺ 1 times from t1 to tN, but observations are not necessarily available at each intermediate time ti, where t1 ⬍ ti ⬍ tN : Then, formula (13) is still valid if the corresponding weight matrices are set to infinity. The expression on the right-hand side of Eq. (13) may be evaluated efficiently only by performing the algebraic operations starting from the left-most factor to the right. To do so, one must develop a special model, called the adjoint, from the forward model (10). The adjoint model may be considered as the transposed of the original model’s gradient but evolving reversely from the final time to the initial time (Tremolet and Le Dimet, 1997). The task of deriving it is quite difficult. In this paper, we used the automatic code of Gering and Kaminski (1998) for deriving the adjoint model from the forward one. Once evaluated, gradient (13) may be used by an optimization procedure to determine the initial condition X1. The number of variables involved in the optimization is quite large, even for a 1D cloud model. Consequently, we used a special optimization procedure, namely LMBFGS of Liu and Nocedal i
i
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
4. Monte Carlo experiment setup The 1D cloud model formulated in Eq. (10) may be transformed into a stochastic model of the type X i⫹1 F X i ⫹ Gi win
14
where G is a coefficient matrix, and wn a sequence of model errors. The stochastic terms introduced in the deterministic model account for the unresolved interactions of the cloud with both large and small scales. Indirect evidence, i.e. the variability observed in rainfall fields (Georgakakos and Krajewski, 1996), suggests that these unsolved interactions may be described by stochastic terms. We consider perturbations caused by unresolved interactions only in the vertical velocity and temperature deviation. Stochastic effects in the other variables are induced by modi~ u=a ~ at random. fying the entrainment terms 2 A ⫺ A This is equivalent to considering that there is a degree of randomness in the cloud/environment interaction. Specifically, we treat the stochastic component as follows. First, we integrate cloud model (10) one
step forward with the only modification being that any entrainment term multiplier becomes 2 ⫹ wn =a (as opposed to 2=a where wn is a normal random variable with mean zero and standard deviation one. Then, the temperature and velocity are perturbed, Ti⫹1 Ti⫹1 ⫹ 0:5wnt and Wi⫹1 Wi⫹1 ⫹ 0:5Wnw ; where wnw, wnt are random variables with mean zero and standard deviation one. We consider the perturbations to be independent. This formulation of the stochastic terms in our model is similar to those of Georgakakos and Krajewski (1996); Tsintikidis and Georgakakos (1999) and Tao and Soong (1987) who provided indirect evidence pertaining to the suitability of such terms. In our numerical experiments, we specified the environmental conditions for the stochastic 1D cloud model by the sounding shown in Fig. 1. The model was initialized by a thermal perturbation that extends to an altitude of 4 km and consists of 3 K temperature excess. We set all other variables to the corresponding values of the environment. We prescribed the boundary conditions as follows. The vertical velocity was zero at both boundaries, and the rain content was zero at the top, and interpolated from above at the bottom. We set the cloud radius to 5 km, and the computational domain from 0 to 17 km, with the vertical grid size of 500 m and the time step of 5 s. With these initial and boundary conditions, we ran the stochastic cloud model 200 times for different random terms, i.e. different initializations of the random number generator. The reason for choosing 100 Temperature Pressure (mbar)
(1989), which is efficient for a large number of variables. Even when model (10) is perfect, if the observations do not contain complete information, estimate X1 may be subject to errors. In that case, there are multiple optimal solutions for objective function (12). Thus, the optimization procedure finds an arbitrary solution that may be far from the true solution (global optimum). If model (10) itself is inaccurate (aside from incomplete information in the observations), X1 will contain additional errors. In this paper, we investigate how the errors induced in X1 by the inaccuracy of the model influence the forecasting performance. A Kalman filter-based formulation is less sensitive to model errors because objective function (12) is modified to account for an inaccurate model, and model (10) is not used explicitly in gradient calculation (13). However, a Kalman filter-approach is infeasible for large problems because it requires operations with large dimension matrices. The VA framework presented here is based mostly on operations with vectors that are affordable on current computer systems, and this justifies its use versus Kalman filters.
89
200 300 400 500 600 700 800 900 -90.0
Dew point
-60.0 -30.0 0.0 o Temperature ( C)
30.0
Fig. 1. Sounding data. The dashed line represents the temperature of an air parcel raised from the surface. The other two lines are for the temperature and dew point.
90
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
16 14
qr (g/kg)
12 Z (km)
qr (g/kg)
1 2 3 4
10
1 2 3 4
8 6 4 2 16 14
Rain rate (mm/h)
Z (km)
12 10
Rain rate (mm/h)
10 20 30 40
10 20 30 40
8 6 4 2
Ground rain rate (mm/h)
40
30
20
10
0 0
20
40 60 Time (minutes)
0
20
40 60 Time (minutes)
Fig. 2. Examples of the stochastic model output. Shown are cloud variables as a function of time. Top panels are contours of rainwater content (g/kg), middle panels are contours of rain rate (mm/h), and the bottom panels are for surface rain rate. The left column represents the mean, and the right column the standard deviation over 200 realizations.
200 as the number of repetitions that yields statistically significant results comes from the central limit theorem (Hogg and Tannis, 1997). According to this theorem, the confidence in estimating the mean of a normal random variable increases proportionally with the square root of the size of sample set used in estimation. A sample size equal to 200 leads to a p confidence interval proportional to 1= 200; which seems to be a reasonable tradeoff between accuracy and computational effort.
In Fig. 2, we summarize the mean over the 200 realizations of cloud variables as functions of time. The corresponding rain variables obtained by running the cloud model with the same initial and boundary conditions but without random perturbations, are shown as a reference in Fig. 3. The differences between the rain related variables yielded by the deterministic model (9) and the means of the same variables obtained from the stochastic model are also shown in Fig. 3. We may conclude based on
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
91
Fig. 3. Deterministic model output. Shown variables are organized as in Fig. 2. The left column is now a single run output, and the right column is the difference between the deterministic and the mean of the stochastic model.
Figs. 2 and 3 that, as anticipated, the average output of the stochastic model is different from that of the deterministic one. It is often assumed that the deterministic model may be considered to represent the modeled reality in some “average sense”. However, applying the stochastic model and taking the average of the results is not the same as taking the mean model and applying it. This is because of the model’s nonlinearity. Considering the mean model instead of the average output of the stochastic one
leads to systematic errors. In our case, the systematic errors induced by the actual deterministic model are relatively small compared to those caused by the stochastic variability. To examine the effect of assimilating radar data in view of such variability, we conducted a Monte Carlo experiment. The steps involved were as follows: (1) generate a sample using the stochastic cloud model;
92
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
Table 1 Uncertainty parameters in assimilation experiments Case
Radar data error
Deterministic model
Drop size distribution
C1 C2 C3 C4 C5
No error No error N(0,9) N(0,9) No error
Exact Inexact Exact Inexact Exact
Known Known Known Known Unknown
(2) synthesize three consecutive, 4.16 min apart (50 times the time step which is 5 s) radar observations; (3) assimilate the radar observations in the deterministic cloud model; (4) forecast the cloud variables using the deterministic cloud model; (5) analyze the forecast performance; (6) repeat steps (1) through (5) 200 times; (7) characterize statistically the results of step (5). We ran the stochastic cloud model for 75 min in step (1). After t 75 min the convection was quite weak, and the rainfall rate was small relative to its peak. Consequently, we did not investigate the process beyond that moment. The assimilation interval started at 16.67 min and ended at 25.00 min. There were three sets of radar observations in the assimilation interval, at t1 16:67; t2 20:83 and
Rain rate error (mm/h)
8
4
0
-4
-8
0
20
40 Time (min tes)
60
Fig. 4. Case C1. The thick line represents the mean error, and the shaded region the error within one standard deviation bounds.
t3 25:00 min: We chose time t 16:67 min as the starting point of the assimilation interval because it provided significant information early enough for the forecasts to have utility. The objective function contained only one term, namely J
3 X
⫺1 i Z i ⫺ Z i;obs T WZ;i Z ⫺ Z i;obs ;
15
i1
where Z i was the vector containing the radar reflectivity factors in all grid points at time ti. Superscript “obs” indicates the “observed” values. The radar observations were synthesized from the rain content assuming a gamma drop size distribution (Ulbrich, 1983), with N0 8 × 104 m3 cm⫺1 and m 1:15 and Rayleigh scattering (Doviak and Zrnic, 1993). In a particular case, which will be specified later, we considered the Marshall and Palmer’s drop-size distribution to investigate the effect of uncertainty model relative to that of other types of uncertainty. It was assumed that the radar sampling volumes were the same as the model’s computational cells and simultaneous with model time steps. In other words, no uncertainty caused by the radar spatial and temporal sampling was considered. The Doppler velocity information could not be projected in the cylindrical coordinate system of the 1D cloud model. Therefore, we did not consider it in the objective function, although in 2- or 3D formulations it is usually considered. We assumed the observational errors were not correlated, and consequently, WZ,i was set to the identity matrix. 5. Results We applied the Monte Carlo experimental setup described in the previous section in a series of cases.
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
Ground rain rate (mm/h)
40 Exact Inexact 30
20
10
0 0
20
40 Time (minutes)
60
Fig. 5. Ground rain rate of the exact and inexact models.
They differ by the degree of noise in the radarsynthesized measurements and the accuracy of the model in assimilation. That is, the radar-synthesized observations are corrupted by noise, and the deterministic model is modified to yield results significantly different from its original output. These cases provide means to investigate the effect of radar errors and model uncertainty on the assimilation and forecasting performance. Another major uncertainty factor that affects the assimilation performance is the lack of information about the rain drop size distribution. We also studied its consequences. Table 1 defines the cases investigated. Normal noise with zero mean and standard deviation 3 dB is expressed as N(0,9). In case C1, the radar observations are considered error free. The stochastic model is used to produce different realizations. Radar observations corresponding to these realizations are produced, and the deterministic model is used for assimilation and forecasting. The mean and standard deviation of the differences between the forecasted and “true” cloud variables (velocity, temperature, water vapor, cloud and rain fields) are then calculated for each realization as a function of time. The results for the rain rate are given in Fig. 4. We can see that the relative errors are small, and there is a fairly good agreement between the forecasted and true rain rate. The wavy appearance of the error might be an indication that the assimilation problem is not well posed, i.e. does not have unique solution, stable to small perturbations. However, the rain rate forecasted globally is quite
93
accurate, although it is underestimated in some periods and overestimated in others. The standard deviation of the error is smaller than the variance of the true rain. This indicates that the inability of a deterministic model to follow the complex cloud environment interactions (mimicked here by the stochastic terms used to generate the reference rain), rather than the assimilation procedure, is responsible for the errors. In case C2, the entrainment terms in the deterministic model are artificially increased 1.4 times. The ground rainfall rates of the initial and altered deterministic cloud models are given in Fig. 5. They are significantly different. The altered model produces about 40% less rainfall. The reason for making this arbitrary modification is to investigate the effect of using an inexact model in assimilation and forecasting on QPF procedure. Even the most complex cloud models may be quite poor in replicating real phenomena. Therefore, it is important to characterize the performance degradation introduced by the model inaccuracy. The standard deviation and mean of the forecasting error, defined as the difference between the forecast and truth, are given in Fig. 6. It is interesting to note that although there is bias, this is significantly smaller (about four times) than that of the deterministic model itself. This is an indication that the assimilation procedure compensates for the inaccuracy of its model. In other words, the optimization procedure finds the initial conditions that make the dynamics of the inaccurate model match with that of the rainfall process. The inaccuracy of the deterministic model can be rigorously accounted for by introducing a set of parameters that influences its dynamics. In the current problem, what caused the negative bias is known, but in real applications only the effects may be observed. In those situations, it is possible to introduce heuristically parameters that make the model more flexible. The parameters may be determined optimally along with model’s initial conditions or state variables (Rajaram and Georgakakos, 1989). This is generally an efficient strategy to mitigate the model inaccuracy. However, in the present study, we did not use it because only one assimilation interval (three sets of observations) was considered for each realization and the requirements of the formalism of Rajaram and Georgakakos (1989) were not met.
94
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
Rain rate error (mm/h)
12 C2
C3
C4
C5
6
0
-6
Rain rate error (mm/h)
-12 12
6
0
-6
-12
0
20
40 60 Time (minutes)
0
20
40 60 Time (minutes)
Fig. 6. Errors of the forecasted rain rate for cases C2, C3, C4 and C5. The thick line represents the mean error, and the shaded region the error within one standard deviation bounds.
Despite this, the assimilation yielded satisfactory results even with an inexact model. Cases C3 and C4 are similar to C1 and C2, but include noise in observations. Gaussian noise with zero mean and 3 dB standard deviation is added to every simulated reflectivity value in the computational domain. The results, given in Fig. 6, are comparable to those of cases C1 and C2, although an increase in errors is noticeable both in means and standard deviations. In case C5, two different gamma drop size distributions are used to express the reflectivity as a function of rainwater content. The reflectivity observations are synthesized based on N0 8 × 104 m3 cm⫺1 and m 1:15 and calculated using the same N0 but m 0:0: The first formula gives a slower growth of reflectivity with rain than does the second one. This implies that for the same reflectivity values, the rain rate estimated
by the second formula is smaller than that of the first one. Therefore, it is expected that the rain content is underestimated in assimilation, and consequently, in forecasting. Fig. 6 confirms this. Note that the error is significantly larger than that of C1, and similar to that of C4. This leads to the conclusion that the uncertainty in raindrop size distribution has an impact as big as that of using a significantly different model in assimilation.
6. Summary and conclusions In this paper, the effect of uncertainty on variational assimilation of radar data in cloud models is investigated. This is accomplished by using Monte Carlo simulations. A 1D cloud model with appropriate modifications to account for the complex interactions
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
between the cloud and environment is used. The 1D model is altered with random terms and run to produce particular realizations of precipitation events. These events are used as references, to synthesize radar observations of reflectivity, and to evaluate the assimilation and forecasting performance. The assimilation is performed in conditions of several kinds of uncertainty, namely, systematic errors in the assimilation model, observation errors and lack of knowledge about rain drop size distribution. The results show that assimilation based the 1D model is quite robust, i.e. it is mathematically well conditioned (it has a solution in conditions of several kinds of uncertainty). We found that the assimilation is able to compensate for systematic errors in the model. Even significant errors in radar observations (normally distributed with zero mean and 3 dB standard deviation) degrade the performance only slightly. However, the assimilation is quite sensitive to the uncertainty in reflectivity introduced by the rain drop size distribution. Based on these results, we conclude that VA of radar data in cloud models may be effective in rainfall forecasting, despite uncertainties in the model. In addition to errors caused by the model uncertainty, other errors may occur in 2- or 3D cloud models, either from the ill-posed character of the problem (multiple solutions for the initial conditions), or from the optimization procedure’s inability to find the best initial conditions (Sun and Crook, 1997). The results in this paper suggest that it may be beneficial for forecasting to reduce the complexity of the models to the level where a unique, optimal estimation of the initial conditions may be obtained. As we demonstrated here, the VA, when mathematically well defined, compensates for the uncertainty of the model used. This issue needs to be further investigated.
Acknowledgements This study was supported by NASA Grant NAG81425 under the auspices of the US Weather Research Program. The radar data were provided under the Cooperative Agreement between the Office of Hydrology of the National Weather Service and the Iowa Institute of Hydraulic Research (NA47WH0495). We thank Drs Konstantine Georgakakos and Mark French for useful discussions.
95
References Bao, J.W., Warner, T.T., 1993. Treatment of on/off switches in the adjoint method: FDDA experiments with a simple model. Tellus 45A, 525–538. Chao, W.C., Chang, L.P., 1992. Development of a four-dimensional variational analysis system using the adjoint method at GLA. 1. Dynamics. Monthly Weather Review 20 (8), 1661–1673. Doviak, R.J., Zrnic, D., 1993. Doppler Radar and Weather Observations. Academic Press, San Diego, CA (458 pp.). Ferrier, B.S., Houze, Jr., R, A., 1989. One-dimensional timedependent modeling of GATE cumulonimbus convection. Journal of the Atmospheric Sciences 46, 330–352. Georgakakos, K.P., Krajewski, W.K., 1996. Statistical-microphysical causes of rainfall variability in the tropics. Journal of Geophysical Research 101 (26), 165–180. Gering, R., Kaminski, T., 1998. Recipes for adjoint code construction. ACM Transactions on Mathematical Software 24 (4), 437– 474. Hogg, R.V., Tannis, E.A., 1997. Probability and Statistical Inference. Prentice Hall, Upper Saddle River, NJ, p. 782. Holton, J.R., 1973. A one-dimensional cumulus model including pressure perturbations. Monthly Weather Review 101, 201– 205. Houze Jr., R.A., 1993. Cloud Dynamics. Academic Press, New York. Kesller, E., 1969. On the distribution of water substance in atmospheric circulations. Meteorology Monographs, Vol. 10. American Meteorological Society, Boston, MA (84 pp.). Liu, D.C., Nocedal, J., 1989. On the limited memory BFGS method for large-scale optimization. Mathematical Programming 45, 503–528. Ogura, Y., Takahashi, T., 1971. Numerical simulation of the life cycle of a thunderstorm cell. Monthly Weather Review 99 (12), 895–910. Rajaram, H., Georgakakos, K.P., 1989. Recursive parameter estimation of hydrologic models. Water Resources Research 25 (2), 281–294. Sun, J.-Z., Crook, N.A., 1996. Comparison of thermodynamic retrieval by the adjoint method with the traditional retrieval method. Monthly Weather Review 124 (2), 308–324. Sun, J.-Z., Crook, N.A., 1997. Dynamical and microphysical retrieval from Doppler radar observations using a cloud model and its adjoint. 1. Model development and simulated data experiments. Journal of the Atmospheric Sciences 54 (12), 1642–1661. Tao, W.-K., Soong, S.-T., 1987. Statistical properties of a cloud ensemble: a numerical study. Journal of the Atmospheric Sciences 44, 3175–3187. Tsintikidis, D., Georgakakos, K.P., 1999. Microphysical and largescale dependencies of temporal rainfall variability over a tropical ocean. Journal of the Atmospheric Sciences 56 (5), 724–748. Tremolet, Y., Le Dimet, F.-X., 1997. Parallel algorithms for variational data assimilation and coupling models. Parallel Computing 22 (5), 657–674. Ulbrich, C.W., 1983. Natural variations in the analytical form of raindrop size distribution. Journal of Climate and Applied Meteorology 22, 1764–1775.
96
M. Grecu, W.F. Krajewski / Journal of Hydrology 239 (2000) 85–96
Verlinde, J., Cotton, W.R., 1993. Fitting microphysical of nonsteady convective clouds to a numerical model: an application of the adjoint technique of data assimilation to a kinematic model. Monthly Weather Review 121, 2776–2793. Yau, M.K., 1980. A two-cylinder model of cumulus cells and its
application in computing cumulus transports. Journal of the Atmospheric Sciences 37 (11), 2470–2485. Zupanski, D., 1997. A general weak constraint applicable to operational 4DVAR data assimilation systems. Monthly Weather Review 125 (9), 2274–2292.