Changes to the NCEP Meso Eta Analysis and Forecast System : Assimilation of satellite radiances and increase in resolution
Eric Rogers, Thomas Black, William Collins, *Geoffrey Manikin, Fedor Mesinger, David Parrish and Geoffrey DiMego
Mesoscale Modeling Branch and *General Sciences Corporation
Environmental Modeling Center, National Centers for Environmental Prediction
1. Introduction
A series of changes to the NCEP Mesoscale Eta Analysis and Forecast System are described. These changes are designed to take advantage of the parallel architecture and increased computing power from the NCEP's new IBM-SP supercomputer. When the operational 32 km Eta model (hereafter referred to as Eta-32) was implemented on the IBM-SP in January 2000, it was converted to run 4 times per day at 6 hour intervals (00Z, 06Z, 12Z, 18Z) with a uniform forecast length of 48-h for all forecast cycles. Each forecast is initialized from the fully cycled Eta Data Assimilation System (EDAS) of 12 hour duration. Further changes were implemented in March 2000 (Manikin et al., 2000), including the extension of the 00Z and 12Z forecasts to 60 hours and modifications to the cumulus convection scheme.
The planned changes consist of the following:
1) Increase in resolution from 32 km / 45 layers to 22 km / 50 layers; increase in the size of the grid domain.
2) Use of a new parallel version of the Eta 3-dimensional variational (3DVAR) analysis, with direct assimilation of GOES and TOVS-1B radiance data, and modifications to better fit rawinsonde moisture data.
3) Inclusion of a non-linear quality control algorithm within the Eta 3DVAR analysis.
4) Addition of vertical advection of cloud water / ice and minor modifications to the Betts-Miller-Janjic cumulus convection scheme.
5) Reduction in horozontal diffusion.
2. Increase in resolution
When the Eta-32 was first implemented into NCEP's operational suite in February 1998, it replaced the 48 km / 38 level version of the EDAS and Eta model. Given the constraints at that time on NCEP's computing resources, it was necessary to reduce the size of the integration domain so that operational products from the Eta forecast which are used by field forecasters were not delayed (see Figure 1 of Rogers et al., 1997). With the increase in computing power provided by the IBM-SP, not only can the resolution be increased from 32 km/45 layers to 22 km/50 layers, but the previous domain size used at 48 km resolution can be restored. This is seen in the plot of the Eta-32 and Eta-22 domains in Figure 1. The greatest impact of the new domain is along the western boundary, where the Eta-22 domain is approximately 900 km further west of Hawaii than the Eta-32 domain.
Figure 2 shows the vertical distribution of the 45 and 50 eta layers. With 50 layers, there is a 10-30% increase in vertical resolution through most of the troposphere, with the largest increase seen between 500 and 750 mb.
3. Use of Satellite Radiances in the Eta 3DVAR Analysis
The 3DVAR analysis for the Eta-22 implementation makes direct use of radiances from the GOES satellites and the NOAA polar orbiting satellites. The algorithm developed for the use of radiances in the NCEP Global 3DVAR analysis has been adapted to work in the EDAS. This code (referred to as OPTRAN) was developed by Tom Kleespies of NOAA's National Environmental Satellite, Data, and Information Service (NESDIS) and is applied without change to the EDAS. OPTRAN requires a 3-dimensional field of ozone, which is not an Eta model prognostic variable. To accommodate OPTRAN, the latest available global model ozone forecast is interpolated to the eta model grid.
The satellite radiances are used to modify temperature and moisture in the analysis. In the EDAS, the Eta model background (a 3-h Eta forecast) is used by OPTRAN to compute a simulated brightness temperature for each observed radiance channel. The brightness temperature difference (observed - simulated) is translated by the 3DVAR analysis through a sensitivity matrix (adjoint) provided by OPTRAN, which indicates how to change temperature and moisture so that the simulated and observed brightness temperatures are more nearly in agreement. For this process to work well, a brightness temperature bias must also be estimated. This is because the radiative transfer model used by OPTRAN is highly parameterized (e.g., it neglects the effect of aerosols). It is important that this bias estimation is always up to date with the current analysis. This will be carefully watched, until a more automated bias correction monitoring scheme is in place. Future upgrades to the global 3DVAR radiance processing algorithm will be implemented as quickly as possible into the Eta 3DVAR code.
The microwave channels, which are most sensitive to moisture, are currently used only over water, as in the global 3DVAR analysis. However, unlike the global 3DVAR, the Eta 3DVAR has and will continue to use SSM/I total precipitable water retrievals over the ocean, since radiances from this satellite are not processed. Although the Eta 3DVAR will use radiances from the GOES moisture channels over the ocean, it will continue to used the GOES 4-layer precipitable water retrievals over land. Poor knowledge of the land surface emissitivity precludes the use of the moisture channels over land.
An evaluation of the impact of radiances has been performed using a lower resolution version of the Eta / EDAS system. Since May 2000 a twice-daily run of a fully cycled EDAS at 48 km / 45 level resolution with satellite radiances has been done, along with a control Eta-48 system identical to the operational Eta-32. Figure 3 shows the 24, 36, and 48-h forecast precipitation skill scores (bias and equitable threat) for control and parallel Eta-48 forecasts for the eastern and western U.S. The eastern / western U.S. regions are defined on the map given in Figure 4. Although the overall impact of the assimilation of radiances on precipitation scores over the contiguous U.S. is small and positive, the regional breakdown shows a slightly greater positive impact in the western U.S. With improved treatment of satellite information by using radiances instead of retrievals, one would expect greater improvement in the western U.S., where satellite observations are the only upstream data source.
The relatively small impact of satellite radiances in the 48-km test, as seen in the precipitation scores and by RMS forecast errors vs. rawinsondes (not shown), is probably due to the combined effects of 1) the relatively low weight given to vertical profiles from satellites in the Eta 3DVAR and 2) the weaker synoptic patterns during the summer months. Additional tests of the EDAS with satellite radiances during winter show that the use of radiances can have a significant impact on the Eta forecast. Figure 5 shows the 24-h accumulated precipitation over northern California valid 1200 UTC 4 February 2000. A strong short-wave trough which came on-shore around 0000 UTC 4 February resulted in copious precipitation in the region, with over a foot of snow (> 1.00 inches water equivalent) observed in the mountains west and northwest of Lake Tahoe. Figure 6 shows the 24-h Eta-32 forecast of precipitation and 500 mb height / vorticity valid 1200 UTC 4 February. The Eta-32 severely underpredicted the intensity of the 500 mb trough, resulting in less forecast precipitation (0.01-0.10 in / 24 h) then what was observed.
Since we expect some improvement in western U.S quantitative precipitation forecasts (QPF) with the use of satellite radiances, the operational EDAS was rerun from 1200 UTC 30 January - 1200 UTC 3 February using satellite radiances in the Eta 3DVAR, followed by a 48-h Eta-32 forecast from 1200 UTC 3 February. Figure 7 shows the 24-h forecast of precipitation and 500 mb heights / vorticity from this rerun. Although the rerun with radiances did not capture the high amounts seen in the mountains, it did predict more precipitation to fall inland, especially north of the Monterey Peninsula. The 500 mb short-wave trough responsible for this event was stronger in the Eta-32 rerun with radiances. A comparison of the operational and experimental Eta-32 first guess and analyses of total column precipitatable water (TPW) valid 1200 UTC 3 February is shown in Figure 8. In both runs analysis of the available moisture data produces a negative moisture increment (i.e, the analysis has lower TPW than the guess) along the offshore frontal boundary. However, the assimilation of radiances for 4 days has resulted in higher TPW in the EDAS first guess with radiances than in the operational EDAS. Although far from perfect, this example shows that use of satellite radiances in the EDAS has the potential to improve QPF forecasts in the western U.S.
A recently observed problem in the operational Eta 3DVAR has been unrealistic maxima / minima in the analysis profile of moisture. An example of this behavior is seen in the comparison of the 0000 UTC 25 July 2000 Nashville, TN sounding and the Eta-32 analysis profile (taken from the eta grid point closest to the station) shown in Figure 9. The sounding shows a mixed layer below 770 mb with constant specific humidity, while the 3DVAR analysis has a wave-like vertical profile of moisture below cloud base. The problem is caused by the uneven distribution of significant levels in a rawinsonde report. If a rawinsonde report has a close grouping of moisture reports, the 3DVAR analysis has a tendency to extrapolate or "overshoot" as it iterates to obtain a solution. Since there can be strong vertical gradients of moisture, the analysis specific humidity profile is more prone to this type of behavior. To overcome this, the 3DVAR analysis has been modified to combine the mandatory and significant level rawinsonde moisture data into a merged profile with a fixed pressure interval of 10 mb, using linear interpolation with respect to the natural logarithm of pressure. The impact of this change can be seen in Figure 10, which shows the Eta-22 sounding at the same time. The use of the merged moisture profile has greatly reduced the irregularities in the analysis moisture sounding. Figure 11, which shows the same two soundings as in Figure 10 plus the 22-km EDAS first guess, reveals that some of the structure in the analysis moisture sounding was brought in by the first guess and was not the result of extrapolation by the analysis.
4. New Quality Control Algorithm
The present NCEP Eta Data Assimilation System uses elements of quality control (qc) that are specific to particular data types. It uses individual qc codes for radiosonde temperatures and heights, aircraft observations, profiler winds, and VAD winds. The EDAS does not use a general optimal interpolation QC (oiqc) as does NCEP's Global assimilation system, but rather presently only has a check on gross data errors.
Variational analysis operates by searching for a minimum to an objective function which measures the square difference of the analysis from all observations and the background (6-hour forecast) state, where the fit is weighted according to the observation and background error statistics. All data mutually influence each other in both the analysis and qc through the observation and background objective functions. The solution for the objective function is iterative, allowing the possibility for qc of the data which adapts to the estimate of the analysis at each iterative step (Ingleby and Lorenc, 1993).
The qc is based upon Bayesian probability theory (Lorenc and Hammon, 1988). The data is assumed to contain a non-Gaussian (gross) component in addition to the usual Gaussian error distribution. This gross error component is assumed to have a flat distribution within some number of standard deviations from the mean. This leads to a modification to the observation part of the objective function which effectively gives data which are far from the analysis a reduced weight (Schyberg and Breivik, 1997, Andersson and Järvinen, 1998). Figure 12 shows a sample assumed data error distribution (dashed curve) and the resulting modified analysis weight (solid curve) due to gross errors, as a function of the normalized data error.
At each iteration in which a minimum to the objective function is more closely approached, each datum is given a weight that reflects its present fit to the analysis. As can be seen from Figure 12, its weight will be nearly 1.0 unless its fit is bad. Even data which receive small weight at one iteration, may later be given large weight if the analysis is brought closer to its value by the influence of surrounding data. Therefore, data are not rejected by this QC strategy, but can regain influence by the support of surrounding data.
One advantage of this QC method is that it does not require the design of a complicated decision making algorithm for its operation, and may be applied equally to all data types. The application of an oiqc is likewise not needed.
The impact of this scheme can be seen from examination of the difference between analyses from runs with and without the non-linear quality control. Such a comparison is seen in Figure 13 for the 1200 UTC 24 January 2000 analyses, which shows the sea level pressure and 500 mb height differences. It is seen that the differences are barotropic in the vertical. Observations (not shown) support the lower sea level pressure analysis using the non-linear quality control. In this particular case, these differences did not lead to any forecast improvement; in fact, there is little difference in the downstream forecast, even at 12 hours.
5. Eta Model Changes
5.1 Changes to Convective Parameterization
Prior to the March 2000 changes, a notable bias in the Eta model was its tendency to predict too much precipitation over the coastal regions in the southeastern U.S. This was caused by the use of drier reference humidity profiles over land in the model's Betts-Miller-Janjic convective scheme. The original intention of this distinction was to attempt to generate more convective precipitation over land. However, this discontinuity at coastlines would cause problems in situations where low-level air parcels with a long trajectory over water achieve a state of convective quasi-equilibrium, consistent with the moister sea profiles. Once these parcels reached land and "sensed" the drier land profiles, the excess moisture is removed in the form of heavy convective precipitation on the coast. This would not only lead to erroneously high precipitation at the coast, but it would frequently prevent the model from producing sufficient precipitation further inland due to the removal of moisture along the coastline.
As described in Manikin et al. (2000), this problem was partially alleviated by using the moister sea profiles at all grid points. This change caused a significant reduction in forecast precipitation bias for amounts at or less than 1 inch / 24 hours over the southeastern U.S. in parallel tests prior to the implementation of this change in March 2000. However, when the land/sea mix of humidity profiles was put into the model, additional features were also added which were intended to reduce heavy convective precipitation over land. These so-called "legacy" features were not removed in the March 2000 implementation:
1) Over land, the humidity profiles revert back to the drier land profiles as long as the intensity of the convective precipitation remains less than 0.25 in / 24 h.
2) The "starting" efficiency parameter in the convection scheme was set at the beginning of each forecast to one over land and to 0.6 over water. It would then evolve during the model integration.
To truly "unify" the convection profiles, in this implementation the convection threshold is reset to zero and the efficiency parameter which evolved during the EDAS is passed into the forecast rather than being reset at the start of the forecast. An independent test of these changes was done in a 48 km version of the Eta model during May 2000. The impact on forecast precipitation (not shown) was small, with modest improvement in equitable threat score for 24-h accumulated precipitation and a 5% increase in precipitation bias score.
5.2 Vertical Advection of Cloud Water
Another Eta model precipitation bias frequently seen is widespread light (~ 0.01 mm / 12 h) precipitation over the eastern North Pacific (Figure 14). This precipitation falls from a persistent and shallow low-level status deck produced by the Eta model cloud scheme over the Pacific. In the operational Eta-32, cloud water is carried in all model layers but is not subjected to vertical advection. When the vertical advection of cloud water was added to the Eta-22 forecast, the amount of light precipitation is reduced (Figure 14). Thus, vertical advection of cloud water will be included in this implementation.
5.3 Horizontal Diffusion
While "horizontal diffusion" is a standard component of atmospheric as well as ocean models, there is no uniformity of viewpoints as to what its objective in a model should be. At least three distinct objectives of "horizontal diffusion" are postulated : controlling small scale noise, maintaining numerical stability, and parameterizing the effects of unresolved on the resolved scales. It is the latter that we see as the purpose of the horizontal diffusion in the Eta. Note that with the resolutions in use horizontal turbulence stresses are unimportant compared to the vertical ones so that horizontal diffusion used by most models is not a representation of turbulence but of the unresolved horizontal advective processes (Mellor 1985). As pointed out by Mellor diffusivities in use are "many orders of magnitude larger" than those which would be appropriate for horizontal turbulence closure, supporting the statement just made.
The variety of formulations in use seem to bear no clear relationship to the objectives stated above. A Smagorinsky-like, deformation dependent formulation may be the most popular; with the so-called 4th order, or biharmonic, schemes used more frequently than the 2nd order (e.g., a summary of 11 models in Doyle et al. 2000; see also Griffies and Hallberg 2000). What we wish to point out is that it is the second-order formulation that is consistent with the view of the exchanges due to unresolved horizontal advective processes; thus, this is the formulation we use in the Eta model.
In the Eta model, a Smagorinsky-like deformation-dependent coefficient is in place (Janjic 1990). Inertial subrange considerations have led Smagorinsky to postulate the coefficient to be proportional to the resolved horizontal deformation times the squared grid distance (e.g., Smagorinsky 1993). Experiments performed at about the time of the operational implementation of the Eta in 1993 with the 80-km and the 40-km resolutions (e.g., Black 1994, Mesinger et al 1993) have led us to believe that this scaling by the squared grid distance was resulting in insufficient diffusion at resolutions higher than the 80 km. As a result, it was abandoned for our 48- and 32-km models, keeping the coefficient multiplying the deformation constant instead.
It has become obvious that the diffusion with no grid-distance scaling was too intense at our higher resolutions; thus, a smaller coefficient was chosen for our experimental 10-km runs of several years ago. For the current 22-km experiments and implementation we are choosing scaling proportional to the grid distance. As can be easily verified from the stability condition of the customary forward time differencing of the diffusion equation (e.g., Mesinger and Arakawa 1976, p. 20) scaling proportional to the grid distance when used with a constant Courant number results in the stability condition independent of the grid distance. Thus, this clearly is a choice of some appeal. It results in a 22-km diffusion coefficient of about 0.275 of that which we would have had we retained our original no-grid-distance-scaling coefficient of the 48- and the 32-km resolutions.
6. Forecast example : 2-3 August 2000 Bering Sea cyclone
One example which shows the benefits of the Eta-22 system over the Eta-32 is discussed below. During 1-3 August 2000 a cyclone formed south of the Aleutian Islands and explosively deepened as it moved into the Bering Sea. A comparison of the forecasts of this cyclone by the NCEP AVN model and the Eta-32 (not shown) revealed that the AVN did a much better job of both initializing the storm and forecasting the explosive deepening. The Eta-AVN central pressure difference was large as 15-20 mb for the 36-48 h forecasts from the 0000 UTC 1 August cycle.
To paraphrase Carven Scott, the SOO at the NWSFO in Anchorage, the Eta-32 analysis did a much poorer job in capturing details in the antecedent jet stream pattern (personal communication). The forecast discussions from the Anchorage NWSFO made it clear that the AVN was the "model of choice" based on run-to-run consistency and a consistently deeper (and better) sea level pressure forecast.
Figure 15 shows the 12-h forecast from the Eta-32 and the Eta-22 valid at 0000 UTC 2 August 2000, while Figure 16 shows the 24-h and 36-h forecasts valid at 0000 UTC 3 August. The observed surface winds and sea level pressure are shown in Figure 17. In each forecast example the Eta-22 predicted the cyclone to be 6-8 mb deeper than the Eta-32 and tracked it closer to the observed position. Since a parallel Eta-48 run which assimilates satellite radiances (not shown) also had similar improvement in its forecast of this storm, it is likely that the improved prediction of this event in the Eta-22 was largely due to assimilation of satellite radiances in the 22 km EDAS.
7. Verification : Precipitation skill scores and fits to rawinsondes
A real-time parallel test of the Eta-22 is ongoing as this bulletin is being written and the quantitative evaluation of its performance is still under way. This will also include a cool-season retrospective test of the Eta-22 for the period 17 January - 3 February 2000. During this interval several significant events occurred, including the Eastern U.S. blizzard of 25 January, a significant precipitation event in California on 25-26 January, and the Eta model forecast bust of 4 February discussed in Section 3. When it is completed the results of the evaluation will be available online at http://wwwt.emc.ncep.noaa.gov/mmb/mmbpll/eta22tpb.stats/ .
8. Conclusions
A series of changes to the operational Meso Eta analysis and forecast system has been described. These changes have been designed to take better advantage of the recent and future upgrades to NCEP's computing power and to address problems seen in the EDAS analysis and Eta forecasts.
The target date for this implementation is pending, awaiting formal approval from the NWS Committee on Analysis and Forecast Techniques Implementation (CAFTI), completion of the cool season retrospective and the quantitative evaluation. It is targeted for the latter half of September 2000.
9. External Qualitative Evaluation of the Eta-22
I. NCEP's Hydrometeorological Prediction Center, late June - mid July 2000 (compiled by Peter Manousos)
1. Introduction
In order to gather data comparative data on the operational performance and utility of the Eta 22km model vs the operational Eta, EMC invited HPC to evaluate the Eta 22km during its parallel production during the period of mid June through mid July.
2. Methodology
HPC assessed the operational performance of the Eta 22km vs the operational Eta by viewing the model using 2 methods:
a. Hard copy 4 panel printout of the Eta 22km mimicking that of the operational Eta. The 4 panel hard copies provided the following information in 12 hour increments.
- 500 HPA Geopotential Height and Absolute Vorticity
- Mean Sea Level Pressure and 1000-500 HPA Thickness
- 700 HPA Geopotential Height/Vertical Motion and Mean 1000-500 HPA Relative Humidity
- 12 hour accumulated precipitation in inches
b. Nominal use of NMAP restore files
Note: Due to the timeliness of the Eta 22km, true parallel run to run comparisons to the operational Eta were not able to be made.
3. Results
The following is a subjective evaluation of the Eta 22km by the operational staff at HPC.
A. Positives
1. The model has performed as well as the operational Eta with respect to capturing the synoptic scale pattern via the mass fields (the mass fields are more noisy as compared to the operational model but are still capturing the pattern).
2. The model has not had any convective feedback problems
3. The model has forecast small scale precipitation events in the short range (6 - 18 hour range) for both heavy and light events
4. The 00Z July 14th run did EXTREMELY well handling the convection in Ohio the afternoon of July 14th.
B. Negatives
1. The model tends to not depict a large enough area of precipitation (seems to be dry around the core of a precipitation area associated with strong forcing mechanisms)
2. The model seems to be a tad too dry with smaller scale convection associated with weak forcing mechanisms.
3. The model has a tendency to over generate light amounts of spotty precipitation near terrain in a moist environment (i.e.- the Appalachians).
4. The model has on tended to be too far north with nocturnal convection.
C. Other
1. The model's high resolution depicts numerous small scale features on the mass fields which the HPC staff had difficulty determining which was a dominant feature. However, this can be accommodated for by mapping the mass fields to a lower resolution grid. HPC would want some fields to remain on the high resolution grid to retain such features (QPF, temperature, wind, etc.)
II. COMET (Steven Jascourt)
I should first note that my evaluation includes only a small sample size and not a rigorous validation, mostly just a comparison of the most noticeable characteristics of the forecast fields.
Most days, most model fields show only minor subtle differences between the 32 km operational run and the 22 km parallel run. The most noticeable and always present difference is speed maxima in jet cores at 250 mb are 10-20 knots faster in the 22 km by 60 hours - less difference early in the forecast period. The faster speeds agree with the AVN valid at the same time. The shapes and locations of the maxima are consistent with the Eta fields (not the AVN fields), so the pattern
looks like the 32 km runs, just stronger jet cores. 850mb temperatures over the southwest U.S. and into the adjacent high plains often show more detail, noticeable even on the coarsely contoured plots on the web site. Seems plausibly realistic though I don't know how well it verifies Some features, especially showing up in relative humidity, are more filamented and stringy in the 22 km instead of more oval-shaped blobs in the 32 km. This too is visible on the web plots, not just on the detailed native-grid NAWIPS plots. This is more characteristic of features we actually see on satellite images and is a natural consequence of deformation and nonlinear advection in the atmosphere. How the features actually match up against observations I don't know - the more detailed structures look more realistic but may or may not be accurately predicted. The subtle differences can and sometimes do alter where the model triggers convection, and this seems to be doing a better job in the 22 km, though the sample size is small so far.
All in all, it seems there are some positives and no identifiable negatives.
References
Andersson, E. and J. Järvinen, 1998: Variational quality control, ECMWF Technical Memorandum, No. 250, 31 pp.
Black, T. L., 1994: The new NMC mesoscale Eta Model: Description and forecast examples. Wea. Forecasting, 9, 265-278.
Doyle, J. D., D. R. Durran, and Coauthors, 2000: An intercomparison of model-predicted wave breaking for the 11 January 1972 Boulder windstorm. Mon. Wea. Rev., 128, 901-914.
Griffies, S. M., and R. W. Hallberg, 2000: Biharmonic friction with a Smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models. Mon. Wea. Rev., 128, 2935-2946.
Janjic, Z. I., 1990: The step-mountain coordinate: physical package. Mon. Wea. Rev., 118, 1429-1443.
Ingleby, N.B., and A.C. Lorenc, 1993 : Bayesian quality control using multivariate normal distributions. Q. J. R. Meteor. Soc., 119, 1195-1225.
Lorenc, A. C., and P. Hammon, 1988 : Objective quality control of observations using Bayesian methods : Theory and a practical implementation. Q. J. R. Meteor. Soc., 114, 515-543.
Manikin, G., and others, 2000 : Changes to the NCEP Meso Eta runs : Extended range, added input, added output, convective changes. NWS Technical Procedures Bulletin 465, NOAA/NWS. [ Available at http://www.nws.noaa.gov/om/tpb/465.htm and from the National Weather Service, Office of Meteorology, 1325 East-West Highway. Silver Spring, MD 20910.
Mellor, G. L., 1985: Ensemble average, turbulence closure. Advances in Geophysics, Issues in Atmos. and Ocean Modeling. Part B: Weather Dynamics, S. Manabe, Ed., Academic Press, 345-357.
Mesinger, F., and A. Arakawa, 1976: Numerical Methods used in Atmospheric Models. WMO, GARP Publ. Ser. 17, Vol. I, 64 pp.
Rogers, E., and others, 1997 : Changes to the NCEP Operational "Early" Eta Analysis / Forecast System. NWS Technical Procedure Bulletin 447, NOAA/NWS. [ Available at http://www.nws.noaa.gov/om/tpb/447.htm and from the National Weather Service, 1325 East-West Highway, Silver Springs, MD 20910 ]
Schyberg, H., and L.-A. Breivik, 1997 : Objective analysis combining observation errors in physical space and observation space. Research Report No. 46, Det Norske Meterologiske Institutt, Oslo, 35 p.
Smagorinsky, J., 1993: Some historical remarks on the use of nonlinear viscosities. Large Eddy Simulations of Complex Engineering and Geophysical Flows, B. Galperin and S. Orszag, Eds., Cambridge Univ. Press, 1-34.