An Observing System Simulation Experiment (OSSE) to Study the Impact of Ocean Surface Observation from the Micro Unmanned Robot Observation Network (MURON) on Tropical Cyclone Forecast

Abstract

The Micro Unmanned Robot Observation Network (MURON) is a planned in-situ observation network over the surface of West Pacific Ocean, and it is designed to sample high spatial and temporal resolution observations of wind and mass fields over the ocean surface. The impacts of MURON observations for Tropical Cyclone (TC) intensity forecast are investigated using Observation System Simulation Experiments. The regional Ensemble Kalman Filter (EnKF) system of Gridpoint Statistical Interpolation is used with the Advanced Research version of the Weather Research and Forecasting model to conduct OSSEs for typhoon Haiyan (2013) while Haiyan goes through rapid intensification. Assimilating MURON observations improves the TC structure and intensity analysis and forecast. The intensity forecast is improved largely due to the correction of initial vorticity and vertical transport of mass flux. The improvement of intensity forecast is attributed largely to the assimilated MURON wind observations when Haiyan is at the tropical disturbance stage, and then by the MURON mass observations when Haiyan enters the tropical storm stage. In addition, our results show that the quality of moisture analysis is sensitive to the choice of the moisture control variable (CV) in the EnKF system. Using the default pseudo relative humidity (PRH) as the moisture CV degrades the accuracy of the moisture analysis. This is likely due to the neglect of updated temperature field during the nonlinear conversion from the PRH CV to the mixing ratio variable and due to the larger deviation of the PRH from Gaussian distribution. The use of mixing ratio moisture CV mitigates these problems.

1. Introduction

In the past several decades, tropical cyclone (TC) track forecasts have improved dramatically [1,2,3] due to the continuous development of numerical modeling and data assimilation (DA) that improved the representation of the meso- and large-scale environment [4]. While TC track forecasts have improved significantly, TC intensity forecasts have shown less improvement over the same periods [5]. Past studies suggest TC intensity is sensitive to the internal dynamics and moist convection in the TC inner core [6,7,8,9,10]. Due to the chaotic nature of these processes, the errors inherent in the initial conditions can lead to rapid growth of intensity forecast error [11,12,13].

Recent studies have shown that the physical processes in the TC boundary layer play an important role in the development and intensification of a TC. The air–sea exchange of momentum, moisture, and heat flux is an essential source to determine TC intensity [14,15,16]. For example, the surface wind can exchange energy through frictional drag and through the evaporation of ocean water [17,18,19,20], and the surface moisture can modulate formation of a deep convective vortex structure in the boundary layer [20,21].

Therefore, in order to better understand the TC intensification process and improve the TC intensity forecasts, high-resolution surface observations of the wind and the thermodynamic fields within the boundary layer near the TC inner core are required [22]. Conventional observations from buoys are available to measure the surface and near surface conditions within the ocean boundary layer routinely. However, these observations typically have a 500–1000 km spacing, which is spatially too coarse to resolve the inner core structure of a TC. Due to flight range constraints, the use of remote sensing and dropsondes onboard reconnaissance aircraft to sample the near surface ocean conditions are limited to storms near landfall.

Infrared surface sensitive channels from instruments onboard satellites can derive valuable information near the ocean surface. However, due to the difficulty of retrieving information under cloudy or rainy conditions [23,24], deriving surface information beneath widespread clouds and precipitation in TCs remains a challenge. Scatterometers on aircraft or polar orbiting satellites are used to derive the wind speed and direction at the ocean surface, which is helpful in estimating the intensity of a TC in most weather conditions [25,26,27]. However, polar orbiting satellites re-visit the same spot only roughly every 12 h, leading to limited temporal resolution. The Cyclone Global Navigation Satellite System (CYGNSS) mission uses scatterometers onboard eight satellites to observe the ocean surface wind, which provides a wide coverage with a high spatiotemporal resolution when compared to individual scatterometer platforms [28,29]. However, scatterometers like CYGNSS are not able to measure the mass field information such as surface air moisture, temperature, and pressure over the ocean [30].

To provide frequent sampling of dynamic and thermodynamic fields near the ocean surface, the Weather News Incorporation (WNI) is developing the Micro Unmanned Robot Observation Network (MURON) system, which can directly measure the surface wind and thermodynamic variables at high resolution over the open water in the West Pacific Ocean. The instrument is expected to function under strong precipitation, such as that associated with TCs. The data collected may improve the accuracy of TC forecasts over the West Pacific Ocean. This study aims to evaluate the effect of the MURON ocean surface observation on TC intensity forecasts using Observing System Simulation Experiments (OSSEs) before mass production of the MURON network designs.

The paper is organized as follows: Section 2 describes the MURON observation system in detail. In Section 3, we describe the configuration of the nature run, and how the MURON observation is simulated using the nature run. Section 4 provides the experimental design, including the model configuration, the DA method, and the experimental configuration. Section 5 describes the OSSEs results for Typhoon Haiyan in 2013, and Section 6 summarizes the results and discusses future work.

2. The MURON Ocean Surface Observation

WNI’s MURON system is an observation network composed of several thousands of micro robots with meteorological sensors to measure surface wind, air temperature, pressure, and humidity (Figure 1). MURON is expected to be deployed in the broad area from 5~25 N° and 130~160 E° with about 30–50 km horizontal resolution and 60 min temporal resolution. The planed deployment region of MURON was determined from accumulated storm track data from 1945 to 2013 from the Japan Meteorological Agency (JMA). This deployment region coincides with the majority of TC genesis and transits. The propulsion system of each robot has been designed to keep the robot at a stationary location in wind speeds up to 40 m s⁻¹ and tidal current up to 3.08 m s⁻¹. These thresholds correspond to Category 1 TCs on the Saffir–Simpson hurricane wind scale. Ocean waves can cause the robots to vibrate and contaminate wind measurements. MURON has been designed to prevent this by applying orientation sensor data with sub-Hertz sampled wind measurements to determine the direction error and remove it. MURON is designed to operate for about three years using its batteries, which are recharged from a solar panel located on top of the robot. Compared to buoy observations, which are existing meteorological sensors over the ocean surface, the design of MURON is optimized in aerodynamics because the only structures exposed to wind are the sensor mast and radio antenna that are thin. Therefore, MURON is less affected by the wind load friction and has less errors in the sensor than the buoy. On the other hand, due to the sophisticated design and advanced sensors of the MURON instrument, individual MURON is expected to be more expensive to maintain than the buoy.

3. OSSE Framework—Nature Run and Simulated Observations

The objective of an OSSE is to estimate the potential impact of new observation platforms on weather and climate numerical predictions before the observing network is fully deployed [31,32]. A well-designed OSSE includes four critical components: the nature run, which is a numerical simulation that is considered the “true” state of the atmosphere in an OSSE, the simulated observations, the assimilation of the simulated observations, and an assessment of the impact of the synthetic observations by comparing forecasts with and without the assimilation. In this section, our case study of Typhoon Haiyan (2013) is discussed, followed by the description of the first two components of the OSSE. The third and fourth components of the OSSE are described in Section 4 and Section 5, respectively.

3.1. Case Description: Tropical Cyclone Haiyan (2013)

To evaluate the impacts of the potential new ocean surface observations from MURON on TC intensity forecasts, an OSSE is conducted for Typhoon Haiyan (2013). Haiyan was selected because its rapid intensification (RI) period could be fully sampled by the simulated MURON observations. Typhoon Haiyan (2013) was the most intense TC and the strongest landfalling TC during the 2013 Northwestern Pacific typhoon season. In addition, Haiyan maintained high intensity over an exceptionally long period, which caused serious human and economic devastation in the Philippines. The Joint Typhoon Warning Center (JTWC) first identified it as a tropical disturbance within the monsoon trough at the south of Pohnpei in the Federated States of Micronesia on 2 November 2013. Because of the favorable large-scale environmental conditions with high SST, high moisture content, and weak shear [33], the tropical disturbance further intensified into a low-level cyclonic circulation within 24 h while moving westward. After the system developed into a tropical storm at 0000 UTC on 4 November, it underwent RI. The system obtained typhoon intensity with a surface maximum wind of 36.01 m s⁻¹ at 0000 UTC on 5 November and quickly intensified to 66.88 m s⁻¹ at 0000 UTC on 6 November while tracking westward [34]. The JTWC assessed Haiyan as a Saffir–Simpson category 5 TC on 6 November and reported that the one-minute sustained wind reached 87.45 m s⁻¹ just before Haiyan made landfall in Eastern Samar of the Philippines at 1800 UTC on 7 November. Haiyan maintained category 5 strength during its passing of the Philippines, then weakened slightly on its way across South China Sea as a 61.73 m s⁻¹ cyclone. After the northward deflection, it continued moving inland into southern China where it rapidly dissipated. Recent studies on Haiyan have shown the influence of equatorial wave disturbances on its genesis [33], sensitivity of intensification to the model physical process [35], climatological characteristics [36], contributions of subsurface ocean warming to the RI [37], and the air–sea enthalpy flux conditions contributing to the RI [38].

3.2. Configuration of the Nature Run

We ran the Advanced Research version of the Weather Research and Forecasting (ARW-WRF) model version 3.6.1 [39] to simulate TC Haiyan for the OSSE. Three domains of dimensions, 361 × 221 (D1), 763 × 421 (D2), and 1420 × 922 (D3), are used for the simulation with horizontal grid spacings of 27 km, 9 km, and 3 km, respectively (Figure 2). The innermost domain (D3) is set to cover the MURON observation area. Both domains have 40 vertical levels with the model top at 50 hPa. The model was configured with WDM6 microphysics [40], Kain-Fritch cumulus parameterization for D1 and D2 [41,42], Yonsei University (YSU) planetary boundary layer (PBL) scheme [43,44], Noah land surface model [45], Rapid Radiative Transfer Model for Global Climate Models (RRTMG; [46,47]), and longwave and shortwave radiation schemes. It is noted that a convection is explicitly resolved in the innermost domain (D3) without cumulus parameterization. To improve the realism of the nature run, a one-dimensional mixed layer model implemented in WRF was used for the nature run simulation [48,49]. The mixed layer model aims to represent realistic sea surface temperatures by capturing the oceanic upwelling and surface cooling after the passage of TC [49]. The one-dimensional mixed layer model, implemented in WRF3.6.1, treats a grid as a single ocean column [49]. The stratification at the top of the thermocline (Γ) and initial mixed-layer depth (h) should be specified as inputs for this model. We examined the sensitivity of the nature run’s TC intensity to these inputs and set h = 30 m and Γ = 0.1 K m⁻¹, which showed the smallest error in TC intensity simulation. These values are comparable to the previous study [49]. Additionally, to prevent drifting of the simulated nature run track from the analysis, Four-Dimensional Data Assimilation (FDDA) nudging [50], which is implemented in WRF3.6.1, was performed on the outer domain (D1) using the National Centers for Environmental Prediction (NCEP) Global Forecast System (GFS) Final (FNL) analysis (NCEP GFS FNL analysis; http://dss.ucar.edu/datasets/ds083.2, accessed on 30 March 2022) for the first 12-h of the simulation. The FDDA nudging is a gridded analysis nudging that relaxes the nature run fields of horizontal wind, temperature, and humidity to the global GFS FNL analysis during the 12-h nudging period. The purpose of this is to improve the TC initialization by reducing the mismatch between the nature run simulation and the global analysis near the boundaries [49,50]. Previous study found that the FDDA nudging improved the hurricane path simulation especially near the domain boundaries [49]. The nudging time of 12-h was determined to best match TC tracks of the nature run and observations. The NCEP GFS FNL analysis provides initial and boundary conditions for the outer domain (D1), and the model is integrated forward for 144-h starting at 0000 UTC on 2 November 2013 (Figure 3). The nature run is generated by the one-way nested simulation on the 3 km innermost domain (D3).

Because this study focuses mostly on the RI of Haiyan, the nature run was configured for maximizing the accuracy of the onset timing and intensification change during RI. Figure 4 shows the intensity of the nature run and the best-track data from the JMA typhoon center for Haiyan. The nature run catches the onset of the RI process at 0600 UTC 4 November and simulates the intensification of Haiyan measured in terms of the minimum Sea Level Pressure (MSLP) until 0600 UTC 6 November. The nature run simulation shows weaker MSLP than the best track data after 0600 UTC 6 November, likely due to the coarse resolution of the simulation (9 km). Although at a relatively coarse resolution, the nature run shows similar RI timing, intensity change, and storm track to the best track data.

3.3. Simulation of Observation

Simulated MURON ocean surface observations and conventional observations including land surface data, marine surface data, radiosonde data from Global Telecommunication System (GTS), and satellite wind data from National Environmental Satellite data and Information Service (NESDIS) are assimilated in this OSSE. The simulated conventional observations are extracted from the nature run and are interpolated to the corresponding locations and times of the real conventional observations. To obtain simulated MURON observations, the zonal and meridional components of wind, temperature, pressure, and relative humidity of the nature run are interpolated to 2 m above the ocean surface, which is assumed to be the elevation of the MURON sensor in this study. Relative humidity is converted to specific humidity using the temperature and pressure observations following that used in the Gridpoint Statistical Interpolation (GSI) DA system. Hereafter, MURON zonal and meridional wind observations are referred to as “wind” observations, and MURON temperature, pressure, and specific humidity observations are referred to as “mass” observations. Random noise is then added to the interpolated fields. The random noise is drawn from a Gaussian distribution with zero mean and a standard deviation equal to the observational errors. For the conventional observations, we adapt the predefined observation errors used in the operational GSI DA system [51,52]. We estimate the MURON representativeness observation error by inflating the MURON instrument error, similar to [53]. The total MURON observation error is calculated by adding the instrument error and representativeness error. The MURON observational error variance was evaluated in terms of the perspective of the DA experiment in Section 4.4.

Table 1. Observation errors of GSI conventional ocean surface observation (GSI_R), MURON (MURON_R), and diagnosed observation error from Desrozier’s method (R_D).

	GSI_R	MURON_R	Diagnosed Observation Error, RD
Temperature	0.92823 °C	0.2 °C	0.2 °C
Barometric pressure	0.536 hPa	1 hPa	0.46 hPa
Specific humidity	34.1 mg/kg	20 mg/kg	17.5 mg/kg
Zonal wind	2.628 m/s	2 m/s	0.48 m/s
Meridional wind	2.628 m/s	2 m/s	0.63 m/s

summarizes the observation error variance of MURON (hereafter, MURON_R). Compared to the variance of conventional ocean surface observation used in the GSI system (hereinafter, GSI_R), MURON has larger observation error for pressure, and smaller observation error for temperature, humidity, and wind. The synthetic MURON dataset is encoded into the Global Data Assimilation System (GDAS) prepbufr files for inclusion in the DA cycles.

4. OSSE Framework—Data Assimilation Experiment

4.1. Model Configuration

The same WRF-ARW model and domain configuration as the nature run are used for the DA control runs. To account for the model errors associated with the physical parameterizations, the nature run and the DA run use a different set of physics schemes to the largest extent possible. Previous studies have shown that PBL scheme, microphysics scheme, and cumulus scheme have a significant influence on the simulated TC intensity [54,55,56]. In this study, because the MURON observations measure ocean surface conditions, the DA and control runs are configured with a different PBL scheme than the nature run to accommodate the model errors caused by the lower-tropospheric thermodynamic and kinematic processes. While the nature run used the complex double moment microphysics (WDM6; [40]), the DA and control runs utilized a simple single moment microphysics scheme [37] to represent model errors associated with the formation of clouds and precipitation. The control run and DA runs have the relatively coarse 9 km horizontal resolution and adopt the Kain-Fritsch cumulus scheme [41,42] while the nature run explicitly resolves convection without the cumulus scheme (D3). The control run and DA runs use the Mellor–Yamada–Janjic (MYJ) PBL scheme [57,58,59,60]. The FDDA nudging and the ocean mixed layer configuration are not used by the control and DA runs. The summary of model configurations for the nature run, the control run, and the DA run are presented in

Table 2. Experimental configuration for nature run, control run, and DA runs.

	Nature Run	Control Run	DA_CONV	DA_CONV_MURON
Domain	27/9/3 km	27/9 km	27/9 km	27/9 km
ICs/BCs	NCEP FNL	NCEP FNL	NCEP FNL/ Perturbed using statistic forecast error covariance	NCEP FNL/ Perturbed using statistic forecast error covariance
Physics parameterization	WDM6	Lin	Lin	Lin
	Kain-Fritsch (only for D1 and D2), None (D3)	Kain-Fritsch (D1, D2)	Kain-Fritsch (D1, D2)	Kain-Fritsch (D1, D2)
	YSU	MYJ	MYJ	MYJ
Initialization	FDDA	no	no	no
Ocean cooling	Mixed layer model	no	no	no
Data assimilation	no	no	EnKF Conventional observations	EnKF Conventional, MURON observations
Forecast time	00 UTC 02~ 00 UTC 8 November (144 h)	00 UTC 03~ 00 UTC 8 November (120 h)	00 UTC 03~ 00 UTC 8 November (120 h)	00 UTC 03~ 00 UTC 8 November (120 h)

.

4.2. Data Assimilation Method

The GSI Ensemble Kalman Filter (EnKF) system was used for this OSSE study, which uses an ensemble square root filter (EnSRF) algorithm [61]. This algorithm has been used in numerous DA studies for TCs [62,63,64,65,66]. The GSI-based EnKF has been adopted as part of the hybrid EnVar system that is used in the operational GFS model [67,68] and the operational HWRF model [65,66]. The EnSRF algorithm is composed of forecast and analysis steps. During the forecast step, 40 ensemble members are integrated for 6-h to estimate a forecast error covariance that will be used in the analysis step.

At the analysis step, the ensemble mean and ensemble perturbations are updated following Equations (1) and (2) [61], respectively.

$${\overline{x}}^a={\overline{x}}^b+K\left(y-\mathit{H}\left({\overline{x}}^b\right)\right)$$

(1)

$$x{\prime }^a=x{\prime }^b-\tilde{K}H\left(x{\prime }^b\right)$$

(2)

where the analysis state vector $x^a$ and the first-guess (background state) vector $x^b$ are composed of the ensemble with a size of K. The ensemble mean of the state vector is denoted by the overbar and the deviation from the mean is denoted by the prime. Here, y denotes the observation, H is the observation operator, and H is the tangent linear of H. K is the Kalman gain defined as $K=P^b{H}^T{\left(H{P}^b{H}^T+R\right)}^{-1}$ where $P^b$ is the background error covariance that is estimated by the background ensemble, and R is the observation error covariance matrix. Here the matrix R is assumed diagonal. $\tilde{K}$ is the modified Kalman gain defined as $\tilde{K}=\mathbf{\alpha }K$ where $\mathbf{\alpha }={\left(\mathbf{1}+\sqrt{R{\left(H{P}^b{H}^T+R\right)}^{-1}}\right)}^{-1}$ [61].

Sampling error associated with a small ensemble size and a misrepresentation of model error can cause an underestimation of the error variances in the EnKF system, which can result in substantial degradation of the quality of the EnKF analysis and subsequent forecast. To alleviate the problem, covariance localization and inflation techniques are introduced to the GSI EnKF system. The covariance localization technique is applied to the Kalman gain in Equation (1) using the [69] function to remove spurious correlation. The localization technique tapers the analysis increment to zero if the distance between a grid point and an observation is larger than 480 km in the horizontal distance and 1 scale height (i.e., the natural log of the pressure) in the vertical distance for both the outer (D1) and inner domains (D2). The horizontal and vertical cutoff radii of localization are determined to minimize the analysis error averaged over the first 8 cycles for the outer domain (D1). These localization radii are comparable with the previous TC study by [65]. A 3-dimensional multiplicative inflation factor is applied to alleviate the sampling error due to the limited ensemble size. The inflation relaxes the posterior ensemble spread to be 85% of the priori ensemble spread [70,71].

4.3. Data Assimilation Experimental Design

For the control and DA runs, initial conditions (ICs) and lateral boundary conditions (LBCs) for ensemble forecast are generated by adding 40 initial ensemble perturbations to the NCEP GFS FNL analysis at the very beginning of the DA cycles valid at 1800 UTC 2 November 2013. The initial and lateral boundary ensemble perturbations for the outer domain (D1) are sampled from a multivariate normal distribution whose covariance is the pre-generated static forecast error covariance of the WRF 3-dimensional variational DA system (3DVar) [72,73]. They were interpolated to the inner domain (D2) to create ensemble perturbations of ICs for the inner domain (D2). The ensemble of LBCs for the later DA cycles is generated in the same manner. The initial ensemble conditions are integrated for 6 h, and the ensemble mean of the 6-h ensemble forecast is used as the initial condition of the control run valid at 0000 UTC 3 November (Figure 3). The DA run assimilates the simulated conventional observations and MURON observations, while the control run does not. Four sets of DA experiments are conducted to evaluate the effect of different combinations of MURON observed variables on the TC intensity forecasts. The DA_CONV experiment assimilates the simulated conventional observations. The DA_CONV_MURON experiment assimilates both the conventional observations and all observing MURON variables, including surface wind, temperature, specific humidity, and surface pressure. The conventional observations and MURON mass variables, such as surface pressure, temperature, and specific humidity, are assimilated in the DA_CONV_MURON_MASS experiment, and only the MURON horizontal wind components, together with the conventional observations, are assimilated in the DA_CONV_MURON_WIND experiment.

4.4. Choice of Observation Error (GSI_R, MURON_R)

Because proper observation error variance is one of the most important factors for successful assimilation of the new observation, the determination of the MURON observation error variance is based on two diagnostics in observation space. One approach is the Desrozier’s method [74],

$$\mathit{E}\left[d_{\mathit{a}}^{\mathit{o}}{\left(d_{\mathit{b}}^{\mathit{o}}\right)}^{\mathit{T}} \right]=R_{\mathit{D}},$$

(3)

where $d_{\mathit{a}}^{\mathit{o}}=y-\mathit{H}\left({\overline{x}}^a\right)$ and $d_{\mathit{b}}^{\mathit{o}}=y-\mathit{H}\left({\overline{x}}^b\right)$ are the observation-minus-analysis residual vector and the observation-minus-background (innovation) residual vector, respectively. E is the statistical expectation operator, and $R_{\mathit{D}}$ is the diagnosed observation error variance matrix from Desrozier’s method. Table 1 shows MURON_R, GSI_R, and $R_{\mathit{D}}$. Compared to $R_{\mathit{D}}$, both MURON_R and GSI_R overestimate the observation errors, except for the pressure observations. However, MURON_R is more consistent with the diagnosed $R_{\mathit{D}}$ for almost all variables, except for the wind observations compared to GSI_R.

The other approach is based on the following innovation statistics [75].

$$\mathit{t}\mathit{r}\left(\mathit{E}\left[d_{\mathit{b}}^{\mathit{o}}{\left(d_{\mathit{b}}^{\mathit{o}}\right)}^{\mathit{T}} \right]\right)\approx \mathit{t}\mathit{r}\left(H{P}^b{H}^{\mathit{T}}\right)+\mathit{t}\mathit{r}\left(\mathbf{R}\right),$$

(4)

which is commonly used to verify the estimation of error variances from an ensemble. The right side can be defined as the total variance and the left side is as the innovation variance. Different observation error variance estimations are evaluated by examining the degree of agreement between the right and left sides of Equation (4). Figure 5 shows the sawtooth plot using the analysis and 6-h forecast at each cycle. MURON_R results in smaller root-mean-square (RMS) of innovations and a more consistent relationship between the RMS of innovations and the total variance than GSI_R for all variables. Based on these two diagnostics, we have configured all DA experiments to use MURON_R as the observation error covariance R in Equations (1) and (2).

5. Results

5.1. Impact of Assimilating MURON Observation on the Accuracy of the Analysis and Forecast

In this section, we evaluate the overall impact of the MURON observation on the analysis and subsequent forecast for all DA cycles. Figure 6 shows the RMS errors (RMSEs) of intensity forecasts averaged over all DA cycles (Figure 6a,b) and for each cycle (Figure 6c,d). In Figure 6c,d, the RMSEs of intensity forecasts are averaged over all the 00-h to 120-h lead times with a 6-h interval at each cycle. Compared to the control run, the RMSE of MSLP and wind speed forecasts of DA_CONV_MURON averaged over all cycles are reduced by 30.2% and 27.5% while the RMSE of the DA_CONV forecasts are reduced by 20.2% and 7.8%, respectively (Figure 6a,b). Assimilating both the conventional and MURON observations improves MSLP and wind speed forecasts for most DA cycles, with a more substantial RMSE reduction at later cycles as the control run, and the DA_CONV has diverged from the nature run at a later time (Figure 6c,d).

The vertical profiles of the RMSE of the control run and DA run analyses are shown in Figure 7. The RMSE is averaged within the 500 km radius from the TC center over all DA cycles. Assimilating the conventional observations reduces the wind and temperature analysis errors at almost all levels (Figure 7a,b). Assimilating MURON observations reduces the wind analysis error significantly below about 400 hPa compared to the conventional observations (Figure 7a). The temperature analysis error is reduced primarily at the lower levels below about 700 hPa compared to the DA_CONV analysis (Figure 7b). Figure 7c shows that assimilating conventional and MURON observations has a neutral impact on the mixing ratio analysis or even degrades the mixing ratio analysis from the surface to the middle level (600 hPa) compared to the control run. This result, in regard to the moisture variable, is unexpected because the surface moisture field and its spatial variability play a crucial role in the intensification of TCs [76], and the assimilation of such observed information should improve the analysis. Sensitivity tests that assimilated each MURON variable show that all observed variables except specific humidity increased the error of the surface mixing ratio analysis (not shown). This result therefore implies that the ensemble estimated cross-error correlation between the moisture variable and non-moisture variables observed by MURON is inaccurate. This issue will be discussed in Section 5.4 where the impact of selection of different moisture control variables is discussed.

5.2. Diagnostics of the Impact of MURON Observation on the RI Forecast

To further diagnose the impact of the assimilation of MURON observations on the RI forecast of Haiyan, the DA cycle initialized at 1800 UTC 3 November 2013 was analyzed. Typhoon Haiyan developed into a tropical storm and continued to intensify and began its RI at 0000 UTC 4 November 2013. This DA cycle is chosen because the 6-h forecast initialized by the DA_CONV_MURON analysis at 1800 UTC 3 November develops into the tropical storm while the control run does not (not shown in the figure). The cyclone valid at 1800 UTC 3 November is not well established in the control run (Figure 8b) and is mis-located in the DA_CONV analysis (Figure 8c). However, the DA_CONV_MURON analysis establishes a primary circulation with a position and an intensity more consistent with the nature run (Figure 8d). In addition, the DA_CONV_MURON analysis shows stronger winds to the north of the TC center and captures the asymmetric wind structure of the nature run (Figure 8d) which is commonly observed in TCs [77,78,79]. Assimilating the conventional and MURON observations increases the magnitude of the wind field below about 700 hPa to the north of the TC center (Figure 8h) while the DA_CONV analysis shows similar vertical structure of the wind with the control run (Figure 8f,g). Assimilating MURON leads to less tilt below 600 hPa, which is more consistent with the nature run (Figure 8e). However, MURON observations have a limited impact on the analysis at mid- to upper-levels above 400 hPa, where nature run shows strong upper-level flow (Figure 8e). In addition, assimilating MURON observations improves the pattern of the upward increase of the equivalent potential temperature (θ_e) near the TC center, which is consistent with structures of stronger TCs (Figure 8h). Although Figure 8 shows that the impacts of MURON are limited at the lower level below about 500 hPa, they affect the vertical tilting of wind and the vertical transport of moisture near the TC center which are favorable for TC intensification within the boundary layer [80].

Surface latent heat flux is a crucial factor determining the vertical extent of TC and plays an important role in TC formation and intensification [81]. Figure 9 shows that the DA_CONV_MURON analysis shows improvement in the spatial distribution and the magnitude of the surface latent heat flux compared to the control run and the DA_CONV analysis (Figure 9b,c). The MURON observations intensify the surface circulation, and the strong surface convergence causes a warm and moist surface air to rise above the ocean surface as shown by the upward increases of θ_e (Figure 8h). The DA_CONV_MURON run condenses and releases more latent heat near the inner core of the storm, which is more consistent with the nature run when compared to the control run.

Figure 10 shows the azimuthally averaged tangential wind, radial wind, and boundary layer height of the analyses valid at 1800 UTC 3 November. The nature run has a wind speed maximum within 40 km and 80 km from the TC center and between 200 m and 1000 m in altitude (Figure 10a), which is similar to that obtained from the composite analysis from past studies [82,83]. The tangential winds of the control run and the DA_CONV analysis are much weaker and the radius of maximum wind is much larger than the nature run (Figure 10b,c). In comparison, both the magnitude and the radius of maximum wind of DA_CONV_MURON are more consistent with the nature run (Figure 10d). In addition, assimilating MURON results in a stronger radial inflow at the lower level below about 1 km, which is more consistent with the nature run when compared to the control run.

The impact of MURON observations on boundary layer structure is evaluated by comparing boundary layer heights against the nature run. The boundary layer height for each column is defined as the level of the maximum wind speed following [83]. The boundary layer height has been identified to have a significant influence on TC intensification through controlling energy transport from the surface layer to the boundary layer above [16,84]. In the nature run, the boundary layer height increases from 10 km to 70 km, levels off between 80 km and 240 km, and increases beyond 250 km (Figure 10a). Assimilating MURON observations enhances the horizontal inflow and improves the vertical wind distribution (Figure 10d), resulting in a boundary layer height closer to the nature run when compared to the control run and the DA_CONV analysis (Figure 10b,d). These results imply that assimilating the MURON observations improves the analysis not only near the surface but also through the lower levels, affecting the TC intensification by regulating the energy transport process from the surface layer to the boundary layer above.

To investigate how the improved analysis leads to more accurate intensification forecast of Haiyan compared to the control run, the mass flux and vorticity are calculated within the TC inner region. The divergence term of the vorticity tendency equation depends on the absolute vorticity and divergence of the flow; therefore, the area averaged absolute vorticity (Figure 11a) and the area averaged vertical transport of mass flux (Figure 11b) are calculated for the analysis and 6-h forecast within 500 km from the TC center. The vertical mass flux is used to estimate the divergence of flow.

Assimilating conventional observations has little impact on the absolute vorticity analysis and the 6-h forecast compared to the control run (Figure 11a). On the other hand, the DA_CONV_MURON analysis has stronger absolute vorticity that is closer to the nature run than the control run and the DA_CONV analysis for all model levels at the initial time valid at 1800 UTC 3 November 2013 (Figure 11a). During the 6-h forecast, like the nature run, the vorticity of DA_CONV_MURON intensifies at almost every level, with its magnitude closer to the nature run than the control run and the DA_CONV forecast (Figure 11a). For the mass flux, the DA_CONV analysis is stronger than the control run below about 500 hPa and is weaker than the DA_CONV_MURON for all levels (Figure 11b). The DA_CONV_MURON analysis at the initial time is stronger than the control run for all levels (Figure 11b), matching closer to the nature run. Below 900 hPa, the DA_CONV_MURON analysis has nearly identical mass flux as the nature run. At the 6-h forecast time, the mass flux of the DA_CONV_MURON run increases, especially at the mid-levels between 700 hPa and 400 hPa, with its magnitude more consistent with the nature run and stronger than the control run (Figure 11b). Overall, the DA_CONV_MURON run has a stronger circulation at almost all levels with greater convergence at low to mid-levels in the analysis than the control run and the DA_CONV. This can contribute to a more rapid evolution of the vorticity than the control run in the subsequent forecast.

While near-ocean surface observations are assimilated, the DA_CONV_MURON analysis has larger vorticity (Figure 11) and higher moisture (Figure 8h) than the control run and the DA_CONV analysis from low to mid-levels below 500 hPa, which are integral for the development of deep convection [85,86]. As shown by the outgoing longwave radiation (OLR) in Figure 12, the DA_CONV_MURON run produces deep convection that is more consistent with the nature run than the control run and the DA_CONV. Specifically, the 6-h forecast of DA_CONV_MURON valid at 0000 UTC 04 shows deep convection cells organized near the TC center (Figure 12d) that are consistent in magnitude and distribution with the nature run than the control run and the DA_CONV (Figure 12). This result suggests that assimilating the MURON surface observation is able to induce vertical mixing to transport energy from the surface to the level above, which aids the development of deep convection. This result suggests that assimilating the MURON surface observations makes vortex stronger, which is a consistent convective structure with the nature run compared to the other.

5.3. Influence of MASS and WIND Observations of MURON

In this subsection, the impact of mass and wind observations from MURON are individually examined. Figure 6 shows that assimilating the conventional observations and the MURON wind observations (DA_CONV_MURON_WIND) improves the intensity forecast more than assimilating the conventional observations and the MURON mass observation (DA_CONV_MURON_MASS) at the early DA cycles before 0000 UTC 4 November 2013, when the simulated Haiyan is at the tropical disturbance (TD) stage (Figure 6c,d). This is likely because the assimilation of wind observations is more effective in correcting an initial TC vortex structure compared to that of mass observations [87]. However, during the later DA cycles after 0000 UTC 4 November 2013 when the simulated Haiyan has developed into a tropical storm (TS) stage followed by the RI, the assimilation of the MURON wind observations has a less positive impact than the assimilation of the MURON mass observations (Figure 6c,d). To investigate the potential source of the intensity forecast errors of DA_CONV_MURON_WIND during the TS stage, we examine the degree of physical balance between the mass and wind fields of the analysis and short-term forecasts. The maximum 10-m surface wind (VMAX) of the analysis, 6-h, 12-h, 18-h, and 24-h forecasts, are compared to the maximum 10-m surface wind (VMAX_AH) which is derived from the pressure–wind relationship of [88]:

$$\mathit{V}\mathit{M}\mathit{A}{\mathit{X}}_{\mathit{A}\mathit{H}}=\mathbf{3.4}{\left(\mathbf{1010}-MSLP\right)}^{\mathbf{0.644}}$$

(5)

where MSLP is the minimum sea level pressure of the analysis, 6-h, 12-h, 18-h, and 24-h forecasts, and 1010 is a constant environmental pressure (hPa) valid for the Western North Pacific [89]. The equation was derived from an empirical relationship between the maximum surface wind and MSLP based on long-term observations. If the analysis and short-term forecasts maintain the balance between mass and wind fields, they should generally follow the relationship of Equation (5) in all the DA cycles. Figure 13 shows the scatter plot of VMAX (red dot) and VMAX_AH (blue dot) as a function of MSLP from the TD (filled circle) and TS (open circle) stages of the DA cycles.

DA_CONV_MURON_WIND produces a larger imbalance between pressure and wind fields than DA_COVN_MURON_MASS, especially when the simulated Haiyan is at the TS stage (Figure 13). As the DA cycle continues, the imbalance between pressure and wind in DA_CONV_MURON_WIND accumulates, contributing to errors in the intensity forecast due to the generation of unrealistic inertia-gravity and its adjustment process during the model integration [90,91] (Figure 6c,d). The positive impact of mass observations during Haiyan’s TS stage demonstrates the potential advantage of the MURON observation platform which samples both wind and mass fields compared to the wind alone remote sensing platforms such as CYGNSS and SATWIND for TC forecasts. These results are consistent with a previous study that showed the benefit of assimilating surface pressure observations together with the wind observations for mesoscale convective system simulations [92].

5.4. Impact of Moisture Control Variable

As discussed in Section 5.1, assimilating the conventional and MURON observations has a neutral or negative effect on the moisture analysis from the surface to mid-levels (Figure 7c). The degradation of the moisture analysis also occurs even when the non-moisture variables of the MURON observations are assimilated. Additional diagnostics to identify the cause are discussed in this section.

Pseudo RH (PRH) is the moisture control variable in the GSI EnKF DA system [51,52]. That is, a background ensemble PRH is updated by assimilating the MURON observations, as shown in Equations (1) and (2). The updated ensemble PRH from the DA is then converted to the WRF prognostic variable, water vapor mixing ratio, before being used to initialize the next WRF background forecast. Mathematically, the PRH ensemble member is defined as

$$\mathit{P}\mathit{R}{\mathit{H}}_{\mathit{i}}=\frac{{\mathit{w}}_{\mathit{i}}}{{\overline{\mathit{w}}}^{\mathit{s}\mathit{b}}}=\frac{{\mathit{w}}_{\mathit{i}}}{{\overline{\mathit{w}}}^{\mathit{s}}\left({\overline{\mathit{T}}}^{\mathit{b}},{\overline{\mathit{p}}}^{\mathit{b}}\right)}$$

(6)

where w, p, and ws are the water vapor mixing ratio, pressure, and saturation mixing ratio, respectively. ${\overline{\mathit{w}}}^{\mathit{s}\mathit{b}}$ is the saturation mixing ratio with respect to the background ensemble mean temperature (${\overline{\mathit{T}}}^{\mathit{b}}$) and pressure (${\overline{\mathit{p}}}^{\mathit{b}}$) where the overbar denotes the ensemble mean. Previous studies have shown that the PRH control variable can better fit background humidity to the observation in a global variational radiance DA system [93,94]; however, the impact of using PRH as the control variable in EnKF-based DA for TC forecasts remains to be explored.

To highlight the issues associated with using PRH as control variable and its transformation, we assume there is a central sea level pressure (SLP) observation that is available for assimilation. Figure 14 shows the DA_CONV_MURON 6-h background ensemble cross-correlation between the SLP at the TC center (marked by black dot) and the near surface water vapor mixing ratio, RH, and temperature throughout the domain valid at 1800 UTC 3 November 2013. Because the analysis ensemble mean mixing ratio has a large error northwest of the TC center (Figure 14d), this focused region (indicated by the black box) is selected for detailed diagnostics. The background SLP has a negative correlation with the background mixing ratio and temperature ensemble northwest of the TC center (indicated by the black box) (Figure 14a,b). This suggests that when the SLP observation innovation is negative, indicating a stronger observed TC, both the temperature and the mixing ratio are higher. This is consistent with known TC structure. On the other hand, the correlation of RH and SLP has the opposite sign. Because the elevated temperature increases the saturation mixing ratio (w^s), RH will relatively decrease, despite the increase in the mixing ratio northwest of the TC center (Figure 14b).

To further reveal how using the PRH moisture control variable can lead to large moisture analysis error at the northwest of the TC center, a single SLP observation is assimilated at the TC center using the EnKF DA with the PRH moisture control variable valid at 1800 UTC 3 November 2013. Hereafter, we refer to this configuration as EnKF_PRH. Figure 15 shows the scatter plots and ensemble estimated probability density function (PDF) of the SLP (y-axis) and the surface moisture variable (x-axis) in the region enclosed in the black box. The red and blue colors represent the analysis ensemble and background ensemble, respectively. The SLP ensemble was sampled from 49 points around the TC center, so a total of 1960 samples (NS) represent the distribution. The same number (NS) of the ensemble of the moisture is sampled in the black box of Figure 14. By assimilating a single SLP observation of 1002 hPa, which has a lower value than the background ensemble mean of 1004.6 hPa, the background ensemble mean of both the SLP and PRH decrease because of the positive cross-correlation of the errors, as indicated by the blue scattered dots (Figure 15a). Figure 15b (red dots) shows the analysis ensemble mixing ratio that is converted from the updated ensemble PRH following Equation (6). The analysis ensemble mean mixing ratio shows a lower value (0.0195 kg kg⁻¹) than the background ensemble mean (0.0199 kg kg⁻¹). This result, however, conflicts with the correlation of SLP and mixing ratio collected from the background ensemble directly (Figure 15b blue dots). The background ensemble shows a negative cross-correlation between these two variables (Figure 15b blue dots). This negative correlation suggests that if the same SLP observation is assimilated, the mixing ratio should be corrected toward a larger value compared to the background. In other words, the mixing ratio will be adjusted toward the opposite direction when PRH is used as the control variable compared to when the mixing ratio is used as the control variable. This is because the background ensemble mean saturation mixing ratio (${\overline{\mathit{w}}}^{\mathit{s}\mathit{b}}$) does not reflect the temperature change in the analysis when converting from the updated PRH to the mixing ratio according to Equation (6). To confirm this hypothesis, we conducted another experiment (EnKF_RH) where the analysis ensemble mixing ratio is generated by multiplying the updated ensemble RH with an analysis ensemble mean saturation mixing ratio (${\overline{\mathit{w}}}^{\mathit{s}}$), not the background ensemble mean saturation mixing ratio (${\overline{\mathit{w}}}^{\mathit{s}\mathit{b}}$). Figure 16b shows that with this configuration, the analysis ensemble mean mixing ratio (0.0203 kg kg⁻¹) is larger than the background ensemble mean (0.0199 kg kg⁻¹), which is consistent with the background ensemble estimated negative cross-correlation between the mixing ratio and SLP. Additional calculations shows that EnKF_RH analysis ensemble mean mixing ratio has a lower RMSE than EnKF_PRH (not shown). This is in contrast to earlier studies which revealed the advantage of using PRH for large scale DA [93,95]. In mesoscale applications, temperature and moisture fields have much larger variability in space and time. Our results suggest that the use of PRH as the moisture control variable in the EnKF DA may not update the temperature and mixing ratio consistently, which can lead to degradation of the moisture analysis.

It is also noted that the error distribution of the background ensemble PRH or RH is highly skewed and non-Gaussian. This deviation from Gaussianity is because RH has both a lower (0%) and an upper bound (100%) (Figure 16a). The non-Gaussian distribution of the background ensemble can produce a suboptimal analysis solution in the EnKF DA [96,97]. To evaluate the effect of a non-Gaussian error distribution of the background RH ensemble on the moisture analysis in the EnKF DA, we assimilate the same SLP observation used in Figure 15 using the Particle Filter (PF) DA method with the RH set as the moisture control variable (hereafter denoted as PF_RH). Because the PF method estimates the model probability densities conditioned on the observations without assuming the background ensemble error distribution [98,99], we expect that the PF DA method will produce a more realistic analysis distribution than the EnKF DA method. Note that PF_RH example presented here only assimilates one observation for demonstration purposes and no treatment on PF degeneration is implemented. The comparison between PF_RH and EnKF_RH reveals how the non-Gaussian error distribution of the background ensemble RH impacts the moisture analysis accuracy in the EnKF DA method.

Figure 17 shows the scatter plots similar to Figure 16, but with results from PF_RH. Although PF_RH and EnKF_RH have almost the same ensemble mean mixing ratio (0.0203 kg kg⁻¹), the updated ensemble RH of PF_RH has a distinct PDF shape compared to that of EnkF_RH (Figure 16a vs. Figure 17a red curves). Such a distinct shape also shows up in the converted mixing ratio (Figure 16b vs. Figure 17b red curves). In the PF DA, assimilating the SLP observation further reduces the background ensemble RH near the maximum peak of the PDF (Figure 17a), while assimilating the SLP observation using the EnKF DA decreases the background ensemble mean RH while maintaining the overall shape of the background ensemble PDF (Figure 16a). This is because the EnKF aims to correct the background ensemble mean (Equation (1)) and to reduce the variance of the background ensemble perturbations (Equation (2)) using an assumption of the underlying Gaussian error distribution. The incorrect estimate of higher moments of the analysis distribution will lead to an incorrect subsequent first guess distribution, which accumulates as the forecast cycles.

Mixing ratio was used as an alternative moisture control variable in an experimental run (EnKF_QV) to examine its impact as a mitigation strategy for the aforementioned issue. Unlike RH or PRH, which have both upper and lower bounds, the background ensemble distribution of mixing ratio only has a lower bound and is closer to a Gaussian distribution (Figure 15a vs. Figure 15b, blue curves). The impact of assimilating MURON observations on the moisture analysis using PRH and QV as control variables is compared (Figure 18). The vertical profile of RMSEs is averaged over all DA cycles, as in Figure 7. The use of mixing ratio instead of PRH as the moisture control variable (EnKF_QV) leads to an improved mixing ratio analysis compared to the control run at lower levels below 800 hPa (Figure 18). Despite the positive impact shown here, preceding studies [93,94] suggest that mixing ratio and its error have a much larger spatiotemporal variability compared to RH, which can lead to errors in the global (non-mesoscale) moisture analysis and the subsequent global forecast. The results presented here indicate that future study is needed to identify the optimal selection of moisture variables for mesoscale DA.

6. Summary and Discussion

This research conducted an OSSE to estimate the effects of potential future new ocean surface MURON wind and mass observations on the RI forecasts using typhoon Haiyan as a case study. MURON consists of a high spatiotemporal resolution observation network, providing in-situ measurements of surface mass and wind fields over the West Pacific Ocean. It is expected that MURON can provide lower-level information near the TC inner core, where the energy flux transfer from the ocean surface plays a critical role for the RI process in TCs. It is found that assimilating the MURON observation together with the conventional observations reduces the TC intensity forecast error by about 30% and 10~20% compared to the control run and the DA_CONV when averaged over all DA cycles. While assimilating conventional observations has a marginal or neutral impact on the analysis compared to the control run, assimilating the MURON observations improves the structure and magnitude of near-surface analyzed fields, including SLP, wind, temperature, moisture, latent heat flux, and PBL height. The improvement of surface analysis fields affects subsequent forecasts from the surface to the middle levels by correcting the initial vorticity and vertical mass flux. These favorable conditions for moist convection in the TC inner core improve the RI forecast and reduce TC intensity forecast errors when compared against the control run and the DA_CONV.

MURON wind and mass observations show distinct effects on the TC intensity forecast during the TD and TS stages of typhoon Haiyan. Assimilating the MURON wind observations reduces the intensity forecast error more than the MURON mass observations at the TD stage through correcting the initial vortex structure efficiently. However, the MURON wind observations cause an imbalance between the wind and pressure fields of the analysis especially, during the TS stage, which reduces its positive impact on the intensity forecast at the TS stage. Compared to the MURON wind observations, the assimilation of the MURON mass observations is less effective in reinforcing the circulation when the inner core structure of typhoon Haiyan is not well developed before the beginning of the RI. However, the assimilation of the MURON mass observations reduces intensity forecast errors more than the MURON wind observations after the onset of the RI process because the mass observations maintain a better balance between the wind and pressure fields of the analysis than wind observations.

The impact of assimilating the MURON observations on the moisture analysis is dependent on the selection of the moisture control for the EnKF update. A neutral or negative impact on the low-level moisture analysis is found when PRH is used as a moisture control variable. This problem is caused by two factors: (1) the updated temperature field is not considered in the nonlinear conversion from the updated PRH moisture control variable to the mixing ratio analysis, and (2) the EnKF can generate a suboptimal moisture analysis error distribution due to the non-Gaussian distribution of the background ensemble PRH, which happens more frequently in a humid environment. The use of mixing ratio as the control variable avoids or alleviates these issues and improves the moisture analyses using the MURON observations.

This study is the first to explore the potential impact of the future ocean surface observations from MURON on the TC intensity forecasts. This initial study is limited by the exploration of a single TC case. Further studies are planned to address how the MURON observation complements other existing observations for TC forecasts.