Final Report for 2004 CWB project

Mesoscale and Microscale Meteorology Division

National Center for Atmospherics Research

Boulder, Colorado, USA

Submitted to Central Weather Bureau, Taiwan, ROC

30 November 2004

1. Introduction

Radio occultation (RO) of the terrestrial atmosphere is now possible through the use of signals transmitted by satellites of the Global Positioning System (GPS) and received by one or more satellites in low Earth orbit. The Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) mission will be launched in late 2005, and about 2,500 ~ 3,000 RO soundings will be provided on a daily basis. This is an important source of data for numerical weather prediction (NWP) in the near future. Meanwhile, NCAR/MMM (National Center for Atmospheric Research, Mesoscale and Microscale Meteorology Division) has developed an advanced 3D-Var (3-Dimensional Variational) data assimilation system over the past few years. Currently, this system has several main features: 1) assimilation of 19 types of observations, including GPS RO local refractivity, radar radial velocity and reflectivity, etc.; 2) multiple choices of control variables (CV) and the corresponding background error statistics (BES): cv_options=2 (NCAR approach) and cv_options=3 (NCEP approach), and the global BES files are provided as default; 3) choices of the minimization algorithms: Quasi-Newton and conjugate-gradient; 4) outer-loop implemented for the incremental approach accounting for the adaptive quality control (QC) and nonlinear effects of the observation operators; 5) ability of being applied to both MM5 model and the WRF model; and 6) code developed in the WRF framework and parallelized with high computing efficiency, and can be run on multiple platforms: IBM-SP, DEC, PC-Linux, Cray X1, Fujitsu VPP5000(?). [Kuo: Why do we put a question mark here?]

Although the GPS refractivity data are not yet available from the COSMIC mission, GPS RO data are available from two existing missions (CHAMP and SAC-C) through CDAAC (COSMIC Data Analysis and Archival Center) in Boulder. The Refractivity data from CHAMP are available from 19 May 2001 to present, and the data from SAC-C are available from 2 August 2001 to 14 November 2002. In preparation for the arrival of the COSMIC data, CWB (Central Weather Bureau) established collaboration with NCAR MMM Division to start the initial experimentation on the assimilation of GPS occultation with the versatile WRF 3D-Var data assimilation system, using CWB observation database and NWP model. It is anticipated that this system will be implemented into CWB operational environment. During the first year (2004) of the CWB-NCAR collaboration, six tasks were performed. This includes establishing the WRF 3D-Var system on the CWB computers, followed by a series of preliminary experiments with CWB operational model domain, CWB observation data base, and the NCEP AVN first guess for an interesting case.

In this report, the accomplishments for the CWB-NCAR collaborative project (through joint efforts between NCAR and CWB staff) are summarized task by task below. The preliminary results from a series of the 3D-Var experiments with WRF 3D-Var/WRF model for Typhoon Dujuan case (30 to 31 August 2003) are presented. Based on these results, the future directions of the work are also suggested.

2. Review of the project tasks and accomplishments

A, Task 1: NCAR/MMM provided the web site to access the WRF 3D-Var code and test data

The WRF 3D-Var code and test data can be downloaded from

http://www.mmm.ucar.edu/wrf/users/download/get_source2.html

From this address, users can also find many links to get the information on the NCAR 3D-Var system, such as documentation, tutorial presentation, and activities, etc.

Ms. Hui-Chuan Lin, CWB staff, has successfully downloaded, implemented, and run 3D-Var system on CWB computers: PC-Linux, Fujitsu VPP5000(?). [Kuo: Why question mark?]

B, Task 2: NCAR/MMM provided the technical consulting support for CWB staff

NCAR staff has had many discussions with Ms. Hui-Chuan Lin during her visit to NCAR from June to September 2004, and helped her run the 3D-Var system: both the observation preprocessor (3DVAR_OBSPROC) and the 3D-Var code (wrf3dvar) itself. In addition to the public release version on the web site, Ms. Hui-Chuan Lin also received the updated versions of the 3D-Var system from NCAR staff.

C, Task 3: Assimilation of GPS radio occultation data

From the COSMIC web site

http://www.cosmic.ucar.edu

the GPS refractivity data (level2/wetPrf) can be downloaded from CDAAC. For the Typhoon Dujuan case, the CHAMP data on the Julian day 242 (30 August 2003, 170 soundings) and 243 (31 August 2003, 168 soundings) has been downloaded and processed. A decoder program: wetprf_decoder, which converts the downloaded file to the ASCII LITTLE_R format file, has been delivered to CWB staff, and it works properly.

D, task 4: 3D-Var observation preprocessor (3DVAR_OBSPROC) upgraded

1) Generalization of the domain definition

Originally, the 3DVAR_OBSPROC was developed for MM5 application, i.e. the domains are defined by a coarse domain and its nests. The conversion of (longitude, latitude) to domain (x, y) is completed by subroutine LLXY. However, the CWB 45-km model domain is a WRF-type domain with the central longitude different from the standard longitude, no information on the coarse domain and nesting were provided. To obtain the correct conversion between (lon, lat) to (x, y), a module map_utils adapted from wrfsi was introduced. With the help of this module, the correct (x, y) coordinates, the coarse domain and nested information were obtained (Fig. 1).

2) Observation errors for GPS refractivity

WMO code FM=116 is assigned to GPS refractivity N. In order to minimize the code changes, the data, N, are placed in the field dew_point in both LITTLE_R and 3D-Var OBS format files. As the initial implementation, the observation error specification is based on Huang et al (2004)

where N_bottom = 10 and N_top = 3 N-units are N error at the bottom and top of the atmosphere. p₀ = 1000 hPa is the pressure at the bottom and p_t = 10 hPa is the pressure at the top of atmosphere. p is the pressure at the observed level. The observation errors are placed in the field dew_point%error. Figure 2 shows the vertical profile of the N errors.

This single profile of the N observation errors is not a good global representation of the RO errors (see below). Kuo et al. (2004) showed that the GPS RO refractivity observation errors can vary with latitudes and seasons. Next year, we may refine the N error specification based on the observation error estimates described in Kuo et al (2004). A data quality control procedure may also be introduced to remove or truncate bad data

3) GPS refractivity data distribution for Typhoon Dujuan case

The distributions of the CHAMP level2 wetPrf data within 6-h time window centered at 1200 UTC 30 and 31 August 2003 are plotted by 3DVAR_OBSPROC (Fig. 3).

E, Task 5: Use of the WRF 3D-Var analysis for CWB’s regional forecast system

1) NCAR staff reviewed the converter code (NFStoWRF, WRFtoNFS, etc.) developed by CWB staff and the resulted data files. In October 2004, we found that the data files were not correct and informed CWB. This helped CWB to fix the database problem. Now, the CWB regional NFS model has been successfully run with WRF 3D-Var analyses.

2) At the request of CWB, NCAR developed the code to output the WRF 3D-Var analysis increments file, including the variables of u, v, t, p, q, and ph, for converting to CWB model initial condition.

3) To reduce the conversion errors, NCAR proposed a procedure through adjustment of water vapor to ensure that the surface pressure P_sfc (total air mass) will not be changed before and after the conversion (this is not yet fully implemented).

§ NFS surface pressure and water vapor

Water vapor pressure Pv in a column is defined as

where r is the total air density, q is the specific humidity. Note in CWB NFS model, the vertical coordinate s is defined by the total pressure p and p_sfc,

§ WRF surface pressure and water vapor

In WRF the vertical coordinate h is defined by the dry hydrostatic pressure. The water vapor pressure in a column is defined as

where r_d is the dry air density, r is the mixing ratio; p_d and pd_sfc are the dry pressure and dry surface pressure.

§ Procedure of the conversion

a) To obtain m for WRF from NFS data

b) To obtain the dry air pressure p_d at the model h levels

c) To obtain the NFS dry pressure p_d(nfs) at the NFS s levels.

The water vapor pressure Pv_k at the level s_k can be calculated as:

d) To obtained the mixing ratio r at the WRF h levels based on p_dk, p_dk(nfs) and r_k(nfs)

Use the linear or the log-p interpolation in the vertical.

e) To obtain the WRF water vapor pressure Pv_wrf in a column based on the mixing ratio r and m, dh, and the WRF surface pressure:

At this point, we may say that the NFS to WRF conversion was completed. But due to interpolation truncation errors, the total water vapor pressure Pv from NFS may be different from Pv_wrf, thus the surface pressure P_sfc from NFS may be slightly different from P_sfc(WRF).

§ Water vapor adjustment

Because Pv_wrf only depends on r, m, and h, and since m and h cannot be modified, the only way to make Pv_wrf=Pv is to adjust the mixing ratio r at the WRF model levels.

Because Pv is a linear function of the mixing ratio, with the new mixing ratio r_adj, the new Pv_wrf will be equal to Pv, and P_sfc(WRF) = P_sfc.

The advantages of the procedures of “NFS->WRF->xb” and “xa->WRF->NFS” are that the analysis from WRF 3D-Var analysis can be used to initialize both the CWB NFS model and the WRF model. However, it will introduce some truncation errors from the conversion between the NFS and WRF. In the future, we should consider the procedures of “NFS->xb” and “xa->NFS” directly.

F, Task 6: Analysis of the WRF 3D-Var experimental results

To ensure the WRF 3D-Var system working properly, the following tasks were conducted.

1) The innovations of the GPS refractivity soundings were plotted. We found that in most cases, the innovations are less than the observation errors assigned here. As we discussed before, the observational errors may be too large, as defined in Huang et al. 2004.

2) The Background Error Statistics (BES) tuning was conducted. Since we do not have BES file from the CWB model, the global NCEP-derived BES has to be used for the 3D-Var experiments. The subjective tuning is performed through the single obs tests for GPS refractivity.

3) A series of the 3D-Var/WRF forecast experiments were carried out. Preliminary results are reported here. For Typhoon Dujuan’s track forecast, the impact of the GPS Ref. data is moderate, but 3D-Var analysis usually gave superior results over the control run.

· Observation operator (Local refractivity observation operator)

The refractivity N is defined as:

where p is the pressure in hPa, T in the temperature in K, and q is the specific humidity in kg/kg (Zou et al 1995). This is a one-step formula for refractivity calculation because p, T, and q are the analysis variables in wrf3dvar.

· The innovations (O-B) of GPS refractivity with NCEP analysis as the first guess and the observation errors

The (O-B) and observation errors for 7 GPS soundings within the 6-h time window centered at 1200 UTC 30 August 2003 are shown in Fig. 4. As pointed out earlier, the observation errors based on Huang et al. (2004) may be too large. The large values of (O-B) located below 8-km and the large negative values near the surface indicate potential data problems. This suggests that certain QC (Quality Control) should be applied prior to the assimilation.

· Background error statistics tuning

Several steps can be taken to tune the background error statistics.

a) Tuning the BES subjectively by looking at the analysis increments plots (when collected enough (O-B) data, tuning the parameters is possible based on the Hollingsworth and Lonnberg (1986) method);

b) Increasing the variances of BES will reduce the (O-A) or increase the magnitudes of the increments;

c) Reducing the scale-length will reduce the influence range of the observations.

The tuning factors for cv_options=3 and cv_options=2 are listed in Table 1.

Table 1 Tuning factors for cv_options=3 and cv_options=2

CV_OPTION = 3				CV_OPTION = 2
Control	Variance Default / Tuned	H-length Default / Tuned	V-length Default / Tuned	Control	Variance Default / Tuned	H-length Default / Tuned
y	0.25 / 0.25	0.75 / 1.0	1.5 / 1.5	y	1.0 / 1.35	1.0 / 0.2
cu	0.25 / 0.25	0.75 / 1.0	1.5 / 1.5	c	1.0 / 1.35	1.0 / 0.2
Tu	0.25 / 0.35	0.75 / 1.0	1.5 / 1.5	pu	1.0 / 1.35	1.0 / 0.2
qp	0.25 / 0.10	0.75 / 1.75	1.5 / 1.5	q	1.0 / 1.00	1.0 / 1.0
Psfcu	0.25 / 0.35	0.75 / 1.0	1.5 / 1.5

For cv_options=3, the background variance factors for the mass variables, T_u and Psfc_u are increased (background weightings decreased) but for moisture, q_p, is decreased (background weighting increased). This means that increments for T_u and Psfc_u will be increased and for q_p decreased. The horizontal scale-lengths are increased for all of the control variables, especially for moisture (1.75), because the default value (0.75) is too small. As a result, the influence range for the mass variables (T and P) is only about 1,500-km (Fig. 5d). With the tuned value (1.0), the influence range became more reasonable, about 2,300-km (Fig. 5b). This helps to spread the GPS data impact more widely in a 45-km resolution CWB domain (Fig. 6a).

For cv_options=2, the default background variance factors are too small, and the scale-length are too large for wind and pressure (Fig. 5c). With the tuned values, more reasonable increments are obtained (Fig. 5a).

· Single OBS tests for GPS refractivity

The cross-sections of the temperature and pressure increments for single GPS refractivity OBS tests for Typhoon Dujuan case in CWB domain are shown in Fig. 5. After tuning, the influence range of an OBS is more reasonable by subjective inspection as mentioned above. A single GPS refractivity OBS (data=1.0 N-unit with the error=1.0 N-unit) is located at the middle level in the center of the domain (see explanation in Fig. 5).

These figures for the single OBS tests are easily plotted with the analysis increment output file from WRF 3D-Var by using a utility program: da_3dvar/utl/convert_to_mm5v3.f90 and MM5/GRAPH.

· Background error for GPS refractivity

Because there are many transforms in 3D-Var system and the GPS refractivity is not one of the direct model variables, it is not straightforward to know what the background errors of the GPS refractivity are in the analysis system. However, based on the following derivation, the background error of the GPS refractivity can be calculated from the (O-B), (O-A), and the observation error s_o. This has already been implemented in the latest WRF 3D-Var code.

For single OBS

B is the background value, O is the observation value, and the

s_b² is the background variance and s_o² is the observation variance.

Then, the analysis A should be

When the (O-B), (O-A), and s_o² of GPS refractivity are known, the background error for GPS refractivity can be computed.

From the single GPS refractivity tests, we obtained the background error, s_b (Table 2). It is clear that the bigger s_b is, the better fit to the observation, with smaller (O-A). With the tuned factors for cv_options=3, the (O-A) = 0.324 and s_b = 1.444 are more realistic in comparing with the case of (O-B) = 1.0 and s_o = 1.0.

Table 2, The background error s_b for (O-B)=1.0 and s_o=1.0

BES	CV_OPTION=3		CV_OPTION=2
	(O-A)	sb	(O-A)	sb
Default	0.171	2.204	0.494	1.012
Tuned	0.324	1.444	0.487	1.027

From the single GPS refractivity tests, we obtained the background error, s_b (Table 2). It is clear that the bigger s_b is, the better fit to the observation, with smaller (O-A). With the tuned factors for cv_options=3, the (O-A) = 0.324 and s_b = 1.444 are more realistic in comparing with the case of (O-B) = 1.0 and s_o = 1.0.

Even when more sophisticate approaches (NMC-method with the ensemble forecasts, the advanced tuning techniques, etc.) are used, subjective inspection for the analysis increment responses to the observation may still be a necessary step prior to the 3D-Var experiments.

· Experiment design

CNTRL : No assimilation, initiated with NCEP AVN analysis

3DREF : Assimilation of GPS refractivity only

3DCWB : Assimilation of CWB observations (SOUND, AIREP, SYNOP. SHIPS, PILOT, SATOB, SATEM, METAR).

3DCWBREF : Assimilation of CWB observations + GPS refractivity.

In the WRF 3D-Var experiments, the first guess fields are obtained from NCEP AVN, and the background error statistics are interpolated from the NCEP global statistics (192 latitudes with resolution of 0.945^o in horizontal and 42 h levels in vertical). The control variables, y, c_u, T_u, q_p and Psfc_u, are selected (CV_OPTION = 3). The initial times are 1200 UTC 30 and 1200 UTC 31 August 2003.

The analysis increments from 3DREF and 3DCWBREF are shown in Fig. 6. From Figure 6, the increments from WRF 3D-Var are reasonable as a result of the assimilation of GPS refractivity and all other types of data.

· Experiment results

Table 3 the averaged track forecast errors for Typhoon Dujuan

Experiment	2003083012Z			2003083112Z
	24-h(km)	48-h(km)	72-h(km)	24-h(km)	48-h(km)	72-h(km)
CNTRL	69.2	67.1	221.1	110.2	115.7	_
3DREF	69.8	85.8	214.2	115.1	120..7	_
3DCWB	58.5	89.3	259.1	103.6	104.9	_
3DCWBREF	52.5	79.8	260.0	105.6	108.1	_

Here we just present some preliminary results of the Typhoon track forecast. Figure 7 shows the track forecasts from the different experiments for Typhoon Dujuan initiated at 1200 UTC 30 and 31 August 2003. From the figures, there are only small differences between the experiments. When the averaged track forecast errors are calculated based on the 3-hourly model output for the periods of 3~24h, 27~48h, and 51~72h, we can see that the 3D-Var experiments, in general, gave better results than the control experiment. The results are listed in Table3.

Bearing in mind that no bogus vortex was used, and the resolution is 45-km for all these experiments, and there are many weakness pointed out earlier (GPS refractivity observation error specification, global background error statistics, etc.) in the current sets of 3D-Var experiments, we cannot expect significant positive impacts from the few GPS refractivity profiles. However, it is encouraging to find that there is moderate impact of assimilating a relatively small number of GPS data. Also, the use of 3D-Var gives positive impact. There is much room for improvements for the assimilation of the GPS data with WRF 3D-Var system. As the number of GPS data is increased (with the availability of the COSMIC data), data quality is improved and the 3D-Var system is fine tuned, we believe that the assimilation of the GPS data will have a positive impact on numerical weather prediction and typhoon forecast in the vicinity of Taiwan.

3, Future work

The basic WRF 3D-Var system has been implemented into the CWB environment, and the preliminary results of assimilating the GPS refractivity are encouraging. In order to assess and improve the performance of the WRF 3D-Var system at CWB, the following tasks should be performed in the near future:

1) Perform semi-operational run with 3D-Var system in CWB

To develop the direct data converters between NFS model and WRF 3D-Var.

To improve the OBS preprocessor incorporating the sophisticate observation error specification and quality control procedure.

To improve the efficiency of 3D-Var system in CWB operational environment (Fujitsu VPP 5000 ?)

To develop the post processing (display and verification) software

2) Continuation of the Typhoon Dujuan case study

To do more experiments with both higher resolution (15-km, 5-km) CWB model and WRF model by using CWB archived data.

To introduce the Typhoon bogus technique in 3D-Var system

To improve the method of determining the forecast Typhoon track, and othermore diagnostic programs.

3) Perform better estimate of the background error and observation error statistics by using more data

To establish the archive system for the semi-operational run (observations, background, innovations, etc.).

To estimate the s_b and s_o for CWB model and observation data.

To derive the CWB-specified BES for different seasons, different resolutions, etc.

4) Assimilation more types of observations, such as satellite radiances, etc.

To develop the BUFR format observation preprocessor for satellite radiance data.

To develop the radiance assimilation code in 3D-Var system

To assess the impact of radiance assimilation by comparing with the assimilation of the retrieved satellite data.

The selection of these tasks for 2005 will be determined jointly by CWB and NCAR based on available resources and project priorities.

References

Huang, C.-Y, Y.-H. Kuo. S.-H. Chen, and F. Vandenberghe, 2004: Improvements in Typhoon Forecasting with Assimilated GPS Occultation Refractivity. Weather and Forecasting, Submitted.

Kuo, Y.-H., T.-K. Wee, S.Sokolovskiy, C. Rocken, W. Schreiner, D. Hunt, and R. A. Anthes, 2004: Inversion and Error Estimation of GPS Radio Occultation data. Special issue of the Journal of the Meteorological Society of Japan, Vol. 82, No. 1B, 507-531.

Zou, X., Y.-H. Kuo, and Y.-R. Guo, 1995: Assimilation of Atmospheric Radio Refractivity Using a Nonhydrostatic Adjoint Model. Mon. Wea. Rev. 123, 2229-2249.

Figure 1, CWB 45-km resolution domain

Figure 2, the vertical profile of the GPS refractivity observation error.

2003083009Z to 2003083015Z 2003083109Z to 2003083115Z

Figure 3, GPS refractivity data from CHAMP level2 wetPrf for Typhoon Dujuan case

Figure 4, Innovations (O-B) of GPS refractivity with NCEP AVN analysis as the background, and the observation errors.

Figure 5, The T/P increments cross-section (West-East about 4,500 km) for Single GPS Ref OBS tests for Typhoon Dujuan case in CWB domain (a) cv_options=2 with the tuning factors: 1.35, 1.35, 1.35, 1.00 for VAR and 0.2, 0.2, 0.2, and 1.0 for LEN; (b) cv_options=3 with the tuning factors: as1=(0.25, 1.0, 1.5), as2=(0.25, 1.0, 1.5), as3=(0.35, 1.0, 1.5), as4=(0.10, 1.75, 1.5), and as5=(0.35, 1.0); (c) cv_options=2 default; (d) cv_options=3 default. The OBS value and error of the Ref. is 1.0 located at x=112, y= 65, z=15. The CWB domain is 222x128x30 with 45-km grid distance.

Figure 6, wind and pressure increments at the lowest model level for Exp. 3DREF and 3DCWBREF.

Figure 7, Typhoon Dujuan track forecast for different experiments initiated at 1200 UTC 30 and 31 August 2003.