1.

## Introduction

The warm tropical ocean acts as a source of extra latent heat and anomalously large values of equivalent potential temperature (i.e. moist entropy) in the boundary layer part of tropical cyclones (TCs), and the low pressures within TCs (Riehl, 1954; Malkus and Riehl, 1960; Ooyama, 1969). The low pressures increase the saturation equivalent potential temperature. High wind and cumulus clouds redistribute upward the extra latent heat acquired at the surface and low levels, leading to warm temperature perturbations aloft and heat loss near the low-temperature tropopause (Emanuel, 1986; Rotunno and Emanuel, 1987). Rotunno and Emanuel (1987) constructed a time-dependent, nonhydrostatic, axisymmetric numerical model with convection explicitly accounted for (hereafter referred to as the RE87 model) to demonstrate that the development of a TC may take place under nearly neutral conditions with each additional increment of latent heat put in by the ocean almost immediately redistributed upward in cumulus clouds. The RE87 model was initialized with an empirically specified model sounding neutral to the model’s cumulus clouds, mimicking the Jordan sounding (Jordan, 1958) that is nearly neutral to actual clouds. If the above idea and the RE87 model of Rotunno and Emanuel (1987) prevail in real TCs, the RE87 model should be applicable to real hurricane intensity forecasts during the time period when the axisymmetric component dominates. The RE87 model can also be used as a better dynamic and physical constraint for vortex initialization (VI) than many previous VI methods.

Earlier VI methods excluded convection so that the role of cumulus clouds as agents for transferring heat acquired at the oceanic heat source upward to the upper tropospheric heat sink is neglected. Detailed procedures differ from one method to another, mainly in the choices of which variable to specify (e.g. wind versus mass fields), what empirical formula to adopt, and which dynamic constraints to impose. Various bogus VI methods can be divided into two major groups depending on whether a wind field (Kurihara et al., 1990; Reeder et al., 1991; Van Thu and Krishnamurti, 1992; Kurihara et al., 1993) or a mass field (Fujita, 1952; Xiao et al., 2000; Zou and Xiao, 2000) is empirically specified. All other model fields are derived from the specified wind or mass field based on the nonlinear balance equation, the gradient-wind relation, the geostrophic relation, the hydrostatic equation, and the divergent mass equation.

Instead of using simplified dynamic constraints, the Pennsylvania State University–National Center for Atmospheric Research nonhydrostatic Mesoscale Model version 5 (MM5) and its full-physics (e.g. cumulus parameterization, planetary boundary layer physics, etc.) adjoint modelling system (Zou et al., 1995; Zou and Kuo, 1996) was employed for generating bogus vortices of all model variables from an empirically specified SLP (Xiao et al., 2000; Park and Zou, 2004). The 2D field of bogus sea level pressure (SLP) is specified using Fujita’s empirical formula (Fujita, 1952) and fitted as a set of prototype ‘observations’ in an MM5 four-dimensional variational data assimilation (4D-Var) framework. The Tropical Prediction Center (TPC)-observed central SLP, the radius of maximum and/or 34-kt surface wind, and the central location of the TC at the model initial time of the target TC are incorporated into Fujita’s formula. The MM5 4D-Var VI generates a dynamically and physically consistent vortex field that is not only automatically merged to the environmental flows, but also well adapted to the TC forecast model, e.g. the MM5. It has been shown that many TC structural features are seen in the 4D-Var-generated bogus vortices, including anomalous warm cores in the middle and upper troposphere, an eyewall, saturated ascent around the eye and descent or weak ascent in the eye, as well as spiral cloud and rainband features. A disadvantage of the MM5 4D-Var VI is its high computational cost and a fit to a bogus SLP.

In this study, we propose to replace the MM5 by the RE87 model and the bogus SLP by satellite-derived warm-core temperature anomalies and total precipitable water (TPW) fields. The novelty of this work involves developing a new 4D-Var VI system, making use of satellite retrieval products of warm core and TPW generated by the authors based on their past research (Tian and Zou, 2016a, 2016b; Zou and Tian, 2018), and demonstrating of a promising, improved TC intensity prediction. The article is arranged as follows. Section 2 briefly describes the nonhydrostatic axisymmetric numerical model with convection accounted for, i.e. the RE87 model and also provides a discussion of the development of the RE adjoint model and the correctness checks of the adjoint model and the gradient of a cost function. Section 3 discusses the satellite microwave retrievals of TC warm cores based on Advanced Technology Microwave Sounder (ATMS) brightness temperature observations from both NOAA-20 and Suomi National Polar-orbiting Partnership (S-NPP) satellites, microwave temperature sounding instrument (MWTS) brightness temperature observations from the Fengyun-3D (FY-3D) satellite, and TPW products obtained from the Advanced Microwave Scanning Radiometer 2 (AMSR2) onboard the Global Change Observation Mission—Water Satellite 1 (GCOM-W1) satellite. Section 4 presents a spectral analysis in cylindrical coordinates of the satellite warm-core and TPW retrievals to show the dominance of the axisymmetric components in these fields. A comprehensive 4D-Var VI assimilating the ATMS-derived warm-core temperatures and AMSR2-derived TPW is carried out for two selected TCs: Hurricane Florence (2018) over the Atlantic Ocean and Typhoon Mangkhut (2018) over the western Pacific Ocean. Results from these two numerical experiments with and without assimilating the satellite warm core and TPW retrievals are presented in Section 5. Section 6 contains concluding remarks.

2.

## Description of the RE97 model and its adjoint model

The RE87 model simulates axisymmetric and compressible flow evolutions on the f-plane. The governing equations are written in cylindrical coordinates (r, ϕ, z) following (Rotunno and Emanuel, 1987):

((1) )
$\frac{\mathit{du}}{\mathit{dt}}-\left(f+\frac{v}{r}\right)v=-{c}_{p}{\overline{\theta }}_{v}\frac{\partial \pi }{\partial r}+{D}_{u},$
((2) )
$\frac{\mathit{dv}}{\mathit{dt}}+\left(f+\frac{v}{r}\right)u={D}_{v},$
((3) )
$\frac{\mathit{dw}}{\mathit{dt}}=-{c}_{p}{\overline{\theta }}_{v}\frac{\partial \pi }{\partial r}+g\left\{\frac{\theta -\overline{\theta }}{\overline{\theta }}+0.61\left({q}_{v}-{\overline{q}}_{v}\right)-{q}_{l}\right\}+{D}_{w},$
((4) )
$\frac{\partial \pi }{\partial t}+\frac{{\overline{c}}^{2}}{{c}_{p}\overline{\rho }{\overline{\theta }}_{v}^{2}}\left\{\frac{1}{r}\frac{\partial \left(\mathit{ru}\overline{\rho }{\overline{\theta }}_{v}\right)}{\partial r}+\frac{\partial \left(w\overline{\rho }{\overline{\theta }}_{v}\right)}{\partial z}\right\}=0,$
((5) )
$\frac{d\theta }{\mathit{dt}}={M}_{\theta }+{D}_{\theta }+R,$
((6) )
$\frac{d{q}_{\mathrm{v}}}{\mathit{dt}}={M}_{{q}_{\mathrm{v}}}+{D}_{{q}_{\mathrm{v}}},$
((7) )
$\frac{d{q}_{\mathrm{l}}}{\mathit{dt}}={M}_{{q}_{\mathrm{l}}}+{D}_{{q}_{\mathrm{l}}},$
where u, v, and w are the radial, azimuthal, and vertical velocities, respectively, π is the nondimensional pressure perturbation from the initial state, θ is the potential temperature, qv is the water vapour mixing ratio, and ql is the liquid water mixing ratio. The ρ denotes the density of air–vapour mixture, c the speed of sound waves, cp the dry air specific heat under constant pressure, θv the virtual potential temperature, g the gravitational acceleration and f the Coriolis parameter. The source-sink terms associated with the phase changes of water vapour are expressed in Equations (5) to (7). R in Equation (5) represents the radiative cooling effect. ${D}_{\chi }$ ($\chi =u,v,w,\theta ,{q}_{v},{q}_{l}$) are diffusion terms. Mathematical expressions for M, D, and R are omitted here and can be found in Rotunno and Emanuel (1987).

A rigid lid is placed on the upper boundary: $w=0$ at ztop = 25 km. Also, at $z={z}_{\text{top}}$, ${\tau }_{\mathit{rz}}=0$, ${\tau }_{z\varphi }=0$ and ${F}_{z}^{\chi }=0$. A sponge layer called the ‘graveyard for old gravity waves’, (zsponge = 20 km), is set to damp out the gravity waves reflected back by the rigid lid by adding a Newtonian damping term $-\alpha \left(z\right)\left(\psi -\overline{\psi }\right)$ to the right-hand side of all prognostic equations except the one for π, where α = 0 for z ≤ zsponge and $\alpha ={\alpha }_{\mathrm{max}}$ at z = ztop. The largest horizontal distance from the TC centre is router = 1500 km. The model resolutions in the radial and vertical directions are Δr = 15 km and Δz = 1 km, respectively. The time step is Δt = 20 s. The model variables are placed on a staggered grid, whose alignment is shown in Fig. 1a. Variables including θ, π, qv and ql also sit on the grids of v.

Fig. 1.

The RE model’s staggered grid alignments in (a) horizontal and vertical coordinates. (b) Variations in $\Phi$(α) with respect to α (plotted using a logarithmic scale), used as correctness checks of the TLM of the RE model.

The forecast at time ${t}_{r}$ from an initial time ${t}_{0}$ made by the RE model can be written symbolically as

((8) )
$\mathbf{x}\left({t}_{r}\right)={\mathbf{Q}}_{r}\left(\mathbf{x}\right){\mathbf{x}}_{0},$
where the vector x includes all seven model state variables in Equations (1) to (7), and x0 represents the model initial condition $\mathbf{x}\left({t}_{0}\right)$. The forward model operator ${\mathbf{Q}}_{r}\left(\mathbf{x}\right)$ in Equation (8) is nonlinear. Linearization of the RE model gives the so-called tangent linear model (TLM), which can be written as
((9) )
$\mathbf{x}\prime \left({t}_{r}\right)={\mathbf{P}}_{r}\left(\mathbf{x}\right)\mathbf{x}{\prime }_{0}=\frac{\partial {\mathbf{Q}}_{r}}{\partial \mathbf{x}}\mathbf{x}{\prime }_{0},$
where the superscript ‘prime’ represents perturbations to the corresponding nonlinear model state variables. The correctness of the RE87 TLM can be checked to see if the following equality is satisfied:
((10) )
$\Phi \left(\alpha \right)=\frac{‖{\mathbf{Q}}_{r}\left(\mathbf{x}+\alpha \mathbf{h}\right)-{\mathbf{Q}}_{r}\left(\mathbf{x}\right)‖}{\alpha {\mathbf{P}}_{r}\mathbf{h}}=1+O\left(\alpha \right).$

The variation in $\Phi \left(\alpha \right)$ with respect to $\alpha$ is provided in Fig. 1b for ${t}_{r}={t}_{0}+1h$. The quantity Φ(α) linearly approaches unity as α decreases logarithmically. When α is less than 10−7, Φ(α) starts to gradually drift away from 1 as a result of machine round-off errors.

The adjoint model is defined as

((11a) )
${\stackrel{\mathbf{^}}{\mathbf{x}}}^{r}={\mathbf{P}}_{r}^{T}\left(\mathbf{x}\right)\stackrel{\mathbf{^}}{\mathbf{x}}\left({t}_{r}\right),$
((11b) )
where $\stackrel{\mathbf{^}}{\mathbf{x}}$ represents the adjoint variable, tR is the final time of interest, and R is the total number of time steps. The correctness of the adjoint model can be checked by the following equality:
((12) )
${\left({\mathbf{P}}_{r}{\mathbf{x}}_{0}^{\prime }\right)}^{T}\left({\mathbf{P}}_{r}{\mathbf{x}}_{0}^{\prime }\right)={\left({\mathbf{x}}_{0}^{\prime }\right)}^{T}{\mathbf{P}}_{r}^{T}\left({\mathbf{P}}_{r}{\mathbf{x}}_{0}^{\prime }\right).$

The left-hand side (LHS) of Equation (12) is an inner product of the vector $\mathbf{x}\prime \left({t}_{r}\right)$, and the right-hand side (RHS) of Equation (12) is the product of ${\mathbf{x}}_{0}^{\prime }$ with the $\stackrel{\mathbf{^}}{\mathbf{x}}\left({t}_{0}\right)$ obtained by an adjoint model integration from the ‘initial’ condition $\mathbf{x}\prime \left({t}_{r}\right)$ at ${t}_{r}$ to its ‘final’ time ${t}_{0}$. With the same setup as in Equation (10) and Fig. 1c, the LHS of Equation (12) is 1914897.7637156348, and the RHS of Equation (12) is 1914897.7637156344, achieving a more than 13-digit accuracy that is close to the machine’s precision. A detailed derivation and instruction on the development of the adjoint model can be found in Zou et al. (1997).

3.

## Satellite microwave retrievals of TC warm cores and TPW

3.1.

### Case selection

The case selected to test the 4D-Var VI system developed in this study is Hurricane Florence (2018) over the Atlantic Ocean. Figure 2 shows the best track of Hurricane Florence from 29 August to 16 September 2018 and its intensity categories during this period. Hurricane Florence experienced a strong intensification from a tropical storm to hurricane category 4 in about 36 h starting around 1530 UTC 4 September 2018. It then quickly weakened to tropical storm intensity at about 1800 UTC 6 September 2018. Figure 2 also shows the 48-h time period from 1530 UTC 4 September to 1530 UTC 6 September 2018, the period of interest in this study.

Fig. 2.

The best track of Hurricane Florence (2018) from 29 August to 16 September 2018 at 6-h intervals. The different coloured symbols indicate the intensity category. The two black arrows show the 48-h time window from 1530 UTC 4 September to 1530 UTC 6 September 2018.

The local equator crossing times (LECTs) of the S-NPP, NOAA-20, GCOM-W1 and FY-3D satellites are quite close. The first three satellites’ LECT is 1:30 pm, and that of the FY-3D satellite is 2:00 pm. Figure 3 presents a single swath of the GCOM-W1 AMSR2, the NOAA-20 ATMS, the FY-3D MWTS and the S-NPP ATMS from 1540 UTC to 1620 UTC 4 September 2018 in the western (Fig. 3a) and eastern (Fig. 3b) parts of the globe. These swaths are superimposed onto GOES-16 Advanced Baseline Imager (ABI) channel-13 (Fig. 3a) and FY-4A Advanced Geosynchronous Radiation Imager channel-12 (Fig. 3b) radiance observations at 1500 UTC 4 September 2018. During this 40-min period, the NOAA-20 ATMS and AMSR2 were observing the equator and the northern hemisphere in the western part of the globe (Fig. 3a) and the high latitudes in the eastern part of the globe (Fig. 3b). The S-NPP ATMS and MWTS2 were observing the southern hemisphere in the eastern part of the globe (Fig. 3b) and the high latitudes in the western part of the globe (Fig. 3a). The centre of Hurricane Florence at 1500 UTC 4 September 2018 was located in the middle of the NOAA-20 ATMS swath and at the left edge of the AMSR2 swath.

Fig. 3.

The swaths covered by GCOM-W1 AMSR2 (red), NOAA-20 ATMS (cyan), FY-3D MWTS (green), and S-NPP ATMS (magenta) from 1540 UTC to 1620 UTC 4 September 2018 in (a) the western and (b) eastern parts of the globe. The underlying grey-scale backgrounds show GOES-16 Advanced Baseline Imager channel-13 (10.35 μm) and FY-4A Advanced Geosynchronous Radiation Imager channel-12 (10.7 μm) radiance observations at 1500 UTC 4 September 2018 [panels (a) and (b), respectively]. The centre of Hurricane Florence at 1500 UTC is indicated by a black hurricane symbol in (a).

The AMSR2, NOAA-20 ATMS, MWTS and S-NPP ATMS swaths can all observe the same hurricane four times in a one-hour time interval twice daily. Figure 4 shows a schematic illustration of the TPW observation times and warm-core retrievals from 1541 to 1639 UTC 4 September 2018. Specifically, the AMSR2, NOAA-20 ATMS, MWTS, and S-NPP ATMS observed the centre of Hurricane Florence at 1541, 1548, 1619, and 1639 UTC 4 September 2018, respectively. Section 5 presents results from a 4D-Var VI experiment with an assimilation window from 1530 to 1639 UTC 4 September 2018 for the case of Hurricane Florence and 1450 to 1545 UTC 9 September 2018 for the case of Typhoon Mangkhut.

Fig. 4.

A schematic illustration of the TPW observation times and warm-core retrievals that are assimilated in the assimilation window from 1530 to 1639 UTC 4 September 2018. Small ticks are drawn at 1-min intervals on the horizontal axis, which corresponds to three time steps of the RE87 model integration. The AMSR2, NOAA-20 ATMS, MWTS and S-NPP ATMS observed the centre of Hurricane Florence at 1541, 1548, 1619 and 1639 UTC 4 September 2018, respectively.

3.2.

### TC warm-core retrievals from microwave temperature sounders ATMS and MWTS-2 (MWTS2)

Hurricane centres are warmer than their environments by more than several degrees, arising from adiabatic warming of the atmosphere if there is downward motion at the eye and a diabatic latent heat release associated with convective precipitation near and outside of the eyewall (Chen and Zhang, 2013). This warm anomaly located at a TC centre is called the warm core. Dolling and Barnes (2012) pointed out that a warm-core formation could signal the transformation of a tropical depression into a hurricane. Based on aircraft reconnaissance data, dropsondes and in situ observations (e.g. buoys, measuring instruments mounted on ships and aircraft, dropsondes and weather stations on islands), TC warm cores are typically found around 250 hPa, but some exist between 750 and 150 hPa depending on intensity and other factors (Hawkins and Imbembo, 1976; Halverson et al., 2006; Durden, 2013). Note that aircraft reconnaissance, dropsonde and in situ observations, i.e. conventional observations, are sparse and not available from deep-ocean areas of the world.

Compared with conventional observations, satellite microwave temperature sounders provide observations with greater spatiotemporal coverage over the global oceans. Of special interest to this study are the ATMSs currently onboard the S-NPP and NOAA-20 satellites, as well as the MWTS2 onboard the Fenyun-3D satellite. The ATMS provides measurements of atmospheric thermal emission in the microwave spectral region between 23.8 and 88.2 GHz with a total of 16 channels. Channels 1 and 2 are window channels at the two lowest frequencies of 23.8 and 31.4 GHz. The MWTS2 has only 13 sounding channels in the microwave spectral region between 50.3 and 89 GHz.

For atmospheric temperature profile retrievals, ATMS channels 5–15 and MWTS2 channels 2–12 are used for TC warm-core retrievals. The weighting functions (WFs) of these channels are vertically distributed evenly in the troposphere and lower stratosphere with the lowest WF peak at ∼850 hPa (the 52.8-GHz channel) and the highest at ∼2 hPa (the 57.29-GHz channel). The temperature-sounding channels with peak WFs below 400 hPa (i.e., ATMS channels 5–7 and MWTS2 channels 2–4) are eliminated from the regression in cloudy conditions to reduce the effect of rain contamination. The beam widths for ATMS channels 5–15 and MWTS channels 2–12 are the same (2.2°). Ninety-six ATMS scene resolution cells, symmetrically distributed on both sides of the sub-satellite path, are sampled. Since neighbouring ATMS FOVs overlap significantly, a remap algorithm is applied to convert ATMS FOVs at the original 2.2° beam width to a new 3.3° beam width (Zou and Tian, 2018).

The warm-core algorithm developed by Tian and Zou (2016) and further refined by Zou and Tian (2018) is applied to ATMS and MWTS2 brightness temperature observations of Hurricane Florence (2018). The regression coefficients are obtained separately for both clear-sky and cloudy conditions from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 reanalysis (Hoffmann et al., 2018) during the two-week training period prior to the occurrence of Hurricane Florence, i.e. 15–30 August 2018.

3.3.

### TPW retrieval from AMSR2

The GCOM-W1 satellite, operated by the Japan Aerospace Exploration Agency, was successfully launched into a polar orbit on 4 July 2012. The GCOM-W1 satellite carries the AMSR2 (Kachi et al., 2008), a conical-scanning microwave imager with 14 channels at the following seven frequencies: 6.925, 7.3, 10.65, 18.7, 23.8, 36.5 and 89.0 GHz (Kawanishi et al., 2003; Zou et al., 2014). Over oceans, AMSR2 radiance observations can be used to retrieve geophysical products such as TPW, the liquid water path, sea-surface wind and sea-surface temperatures (Wentz and Meissner, 2000). Note that the AMSR2 channels with frequencies centred at 18.7 and 10.65 GHz can encounter interference from television and radio frequency interference (TFI and RFI, respectively) signals, which introduce errors in AMSR2 geophysical retrieval products such as TPW over coastal regions. Different from RFI over land, which can be detected but is difficult to correct, Tian and Zou (2016) developed an empirical model for a TFI correction of AMSR2 data over oceans near the U.S. and Europe. False dry anomalies in the 18.7-GHz retrieved TPW in the presence of TFI near the coastal regions were successfully corrected using this model.

4.

## Azimuthal spectral analyses and the axisymmetry of Hurricane Florence

Since the RE87 model is an axisymmetric model, we select a time period when Hurricane Florence had a strong axisymmetric component. An azimuthal spectral analysis is carried out on GOES-16 ABI channel-13 (10.35 μm) radiance observations (Zou et al., 2010). Figure 5a shows a spatial distribution of radiances observed by GOES-16 ABI channel 13 at 1500 UTC 4 September 2018. A high local maximum of 10.35-μm radiance observations is seen as the eye of Hurricane Florence, embedded in a low radiance area due to convection. The amplitudes of azimuthal wavenumbers 0–4 are calculated from the ABI observations in Fig. 5a so that some insight into how much is represented by axisymmetric and asymmetric components can be gained (Fig. 5b). The axisymmetric component accounts for more than 80% near the eye (d 100 km) and about 50% between 150- and 300-km radial distances. The wavenumber one component accounts for the largest asymmetry with a maximum of slightly more than 10% at the 200-km radial distance. Figure 5c shows the temporal evolution of the wavenumber-0 spectrum calculated from GOES-16 ABI channel-13 radiance observations during 1–7 September 2018 at 3-h intervals and the observed temporal evolution of the central SLP. Hurricane Florence is characterized by an increased axisymmetry as it went through intensification during 4–6 September 2018. Due to the strong impact of convection, which has large asymmetric components, on ABI infrared channels, regions of large axisymmetry are confined to relatively small radial distances.

Fig. 5.

(a) Spatial distribution of ABI channel-13 radiance observations at 1500 UTC 5 September 2018 and (b) its radial variations in percentage spectra of the azimuthal wavenumbers 0–4 at a 25-km radial distance and in 2° azimuthal angle intervals. (c) Temporal evolution of the observed central SLP (black curve) and the wavenumber-0 spectrum (i.e. the axisymmetric component) calculated from GOES-16 ABI channel-13 radiance observations during 1–7 September 2018 at 3-hour intervals. Black dashed lines show the analysis time at 1530 UTC 4 September 2018 and the 48-hour forecast time at 1530 UTC 6 September 2018.

The TPW fields within TCs retrieved from AMSR2 observations are mostly axisymmetric. Figure 6a shows an example of the TPW retrieved from AMSR2 observations at 1624 UTC 5 September 2018. Hurricane Florence’s centre (see the cross symbol in Fig. 6a) was located in the middle of an AMSR2 swath at this time. Figure 6b shows the corresponding spectra of the TPW for azimuthal wavenumbers 0–4. The wavenumber-0 component accounts for about 90% of the total spectra from the centre to the 300-km radial distance. The spectra of wavenumbers 1 and 2 vary around 1% and those of wavenumbers 3 and 4 are below 1%. In other words, the microwave retrieval of TPW is much more axisymmetric than the ABI infrared radiance observations. This agrees with the fact that ABI observations are strongly affected by clouds that have spiral-band structures within tropical storms (Anthes, 1982).

Fig. 6.

(a) TPW (unit: kg m−2) retrieved from AMSR2 observations at 1624 UTC 5 September 2018 and (b) its radial variations in percentage spectra of the azimuthal wavenumbers 0–4. The centre for the azimuthal spectral analysis is indicated by ‘x’ in (a).

Fig. 7.

Temperature anomalies at 250 hPa (unit: K) retrieved from observations made by the NOAA-20 ATMS at 1548 UTC (top panel), the FY-3D MWTS at 1619 UTC (middle panel), and the S-NPP ATMS at 1640 UTC (bottom panel) 4 September 2018. The centre from the best track data is indicated by the ‘x’ symbol.

Fig. 8.

Radial variations in percentage spectra of the azimuthal wavenumbers 0–4 for the 250-hPa warm-core anomaly retrieved from NOAA-20 ATMS observations at 1548 UTC 4 September 2018.

5.

## 4D-Var VI experiment design and numerical results

5.1.

### Experiment design and convergence

Only the axisymmetric components of the TPW and warm cores are assimilated in the 4D-Var VI experiments because the RE87 model is an axisymmetric model. The following cost function is minimized in a 4D-Var VI experiment (Navon et al., 1992; Zou et al., 1993):

((13) )
$J\left({\mathbf{x}}_{0}\right)=\frac{1}{2}{\left({\mathbf{x}}_{0}-{\mathbf{x}}_{b}\right)}^{T}{\mathbf{B}}^{-1}\left({\mathbf{x}}_{0}-{\mathbf{x}}_{b}\right)+\frac{1}{2}\sum _{r=0}^{N}{\left({H}_{r}\left({\mathbf{x}}_{r}\right)-{\mathbf{y}}_{r}\right)}^{-1}{\mathbf{O}}_{r}^{-1}\left({H}_{r}\left({\mathbf{x}}_{r}\right)-{\mathbf{y}}_{r}\right),$
where x0 is the analysis vector on the analysis/forecast grid at time t0, xb is the background field, xr is the model state at time tr, yr represents microwave retrievals of TPW and the warm core at time tr, Hr is the observation operator for calculating TPW and the warm-core anomaly from xr, Or is the observation error covariance matrix at time tr and B is a diagonal background error covariance matrix. The first term on the right-hand side is called the background term and is denoted by Jb, and the second term measures the weighted differences between model forecasts and observation retrievals within the assimilation window [t0, tN].

Azimuthally averaged potential temperature anomalies calculated at 25-km radial-distance intervals using temperature retrievals from the NOAA-20 ATMS at 1548 UTC, the FY-3D MWTS at 1619 UTC, and the S-NPP ATMS at 1640 UTC, as well as the AMSR2-retrieved TPW at 1541 UTC 4 September 2018 are taken as ${\mathbf{y}}_{r}=\mathbf{y}\left({t}_{r}\right)$ (r = 1, …, 4), where ${t}_{1}=$1541 UTC, ${t}_{2}=$1548 UTC, ${t}_{3}=$1619 UTC and ${t}_{4}=$1639 UTC (Fig. 9). The azimuthally averaged potential temperature anomalies are positive and negative below and above an altitude h, respectively, where h varies from 13–14 km at the hurricane centre to about 11–12 km at the 500-km radial distance. The altitude range of 11 to 12 km contains the warm cores. TWP values are high within the hurricane (Fig. 9c). The magnitude of TPW dips at the eye, likely due to a downward motion at the eye or to a weaker eyewall upward motion that reduces the transport of low-level water vapour upward.

Fig. 9.

Azimuthally averaged potential temperature anomalies (unit: K) calculated at 25-km radial-distance intervals using temperature retrievals from (a) the NOAA-20 ATMS at 1548 UTC, (b) the FY-3D MWTS at 1619 UTC, and (c) the S-NPP ATMS at 1640 UTC 4 September 2018. The AMSR2-retrieved TPW at 1541 UTC 4 September 2018 (black curve, unit: kg m−2) is shown in (c).

The ECMWF’s next-generation reanalysis ERA5 at 1200 UTC 4 September 2018 after an azimuthal averaging is used as an initial condition for the RE87 model to generate the background field ${\mathbf{x}}_{b}$ at 1530 UTC 4 September 2018. Compared with the ERA-Interim reanalysis, which was introduced in 1979 and has been continuously updated since then, the ERA5 reanalysis has better spatial (31 km horizontally and 137 levels vertically) and temporal resolutions (hourly), allowing for better representations of convective updrafts, gravity waves, TCs and mesoscale-/synoptic-scale features of the atmosphere. This has led to significant improvements in the conservation of potential temperature in the stratosphere (Hoffmann et al., 2018). Figure 10 shows the background potential temperature (Fig. 10a) and water vapour mixing ratio anomalies and mixing ratio (Fig. 10b) of the 3.5-h forecast of the RE87 model initialized by the azimuthally averaged ERA5 reanalysis at 1200 UTC 4 September 2018. A warm core more than 3 K warmer than its environment is located between 6 and 11 km, with a maximum around 9 km (Fig. 10a). There is a wet core in the middle and lower troposphere with a maximum value of more than 6 g kg−1 around 2 km.

Fig. 10.

(a) Potential temperature anomaly (unit: K) and (b) water vapour mixing ratio anomaly (colour shading, unit: g kg−1) and mixing ratio (contours) of the 3.5-h forecast of the RE87 model initialized by the azimuthally averaged ERA5 reanalysis at 1200 UTC 4 September 2018.

Minimization of the cost function defined in equation (13) requires the gradient of J with respect to x0, which can be calculated by the following mathematical expression:

((14) )
${\nabla }_{{\mathbf{x}}_{0}}J={\mathbf{B}}^{-1}\left({\mathbf{x}}_{0}-{\mathbf{x}}_{b}\right)+\sum _{r=0}^{N}{\mathbf{P}}_{r}^{T}{\mathbf{H}}_{r}^{T}{\mathbf{O}}_{r}^{-1}\left({H}_{r}\left({\mathbf{x}}_{r}\right)-{\mathbf{y}}_{r}\right),$
where ${\mathbf{H}}_{r}^{T}$ is the adjoint of the observation operator, and ${\mathbf{P}}_{r}^{T}$ is the adjoint model of the RE87 model. The correctness of the gradient can be checked if the following equality is verified:
((15) )
$\psi \left(\alpha \right)\equiv \frac{J\left({\mathbf{x}}_{0}+\alpha \mathbf{h}\right)-J\left({\mathbf{x}}_{0}\right)}{\alpha {\mathbf{h}}^{T}\nabla J\left({\mathbf{x}}_{0}\right)}=1+O\left(\alpha \right),$
which is similar to Equation (10) used for the correctness check of the TLM.

Figure 11a shows the variation in |Ψ(α) − 1| with respect to the perturbation magnitude parameter $\alpha$. Not only does $\psi \left(\alpha \right)$ approach unity as the perturbation decreases from $\alpha ={10}^{-1}$ to $\alpha ={10}^{-10}$, but the variation is also linear as expected from Equation (15). A further reduction of the perturbation induces an increase in |Ψ(α) − 1| due to round-off errors. An accurate and efficient calculation of the gradient of the cost function using the adjoint operator of the RE87 model by Equation (14) is necessary for the minimization of the cost function defined in Equation (13). Figure 11b presents variations in the cost function and the norm of the gradient of the cost functions with respect to the number of iterations. The cost function decreases by two orders of magnitude in 45 iterations. This represents a significant adjustment of the model fields to fit the warm core and TWP retrievals. Since the RE87 model attempts to fit azimuthally averaged satellite retrievals, we do not expect the cost function to reduce to zero at the end of a minimization procedure. However, we do expect the norm of the gradient to be zero at the end of a minimization procedure for ${\mathbf{x}}_{0}^{*}$ to be the minimum point, i.e.

((16) )
where $\epsilon$ is a vector defining the neighbourhood of ${\mathbf{x}}_{b}$. In the 4D-Var VI experiment conducted in this study, the norm of the gradient decrease by more than three orders of magnitude in 45 iterations (Fig. 11b). The solution found at the 45th iteration, ${\mathbf{x}}_{0}^{\left(45\right)}$, is taken as the minimum: ${\mathbf{x}}_{0}^{*}\equiv {\mathbf{x}}_{0}^{\left(45\right)}$, which will serve as an initial condition for the RE87 model to make a 48-h forecast.

Fig. 11.

(a) Variations in |Ψ(α)1| with respect to the perturbation magnitude parameter α and (b) the cost function (black curve) and the norm of the gradient of the cost function (dashed curve) with respect to the number of iterations.

5.2.

### Numerical results of the analysis and forecasts

Potential temperature anomalies are also calculated by subtracting the mean temperature within the 1500 km radial distance but outside the 34-kt wind radius at each vertical level. Figure 12a shows results from the analysis of the potential temperature anomaly obtained by the 4D-Var VI experiment. Compared to the background distribution (Fig. 10a), the warm core has a maximum at about 10 km altitude that is slightly higher, also stronger in magnitude, and has a more of an oval shape than those in observations. The analysis increments in water vapour mixing ratio are small, with the largest adjustments of ∼1.6 g kg−1 located near 1 km within $±$ 100-km radial distances (Fig. 12b). The TPW values calculated from the 4D-Var VI analysis are greater than the background TWP by about 3 kg m−2 at all radial distances with $±$ 500 km (Fig. 12b). In other words, the 4D-Var VI generated an analysis giving a higher warm core and a stronger wet core than those from the background field.

Fig. 12.

(a) Analysis of the potential temperature anomaly (unit: K) and (b) analysis increments of the water vapour mixing ratio (colour shading, unit: g kg−1) and the background (black curve), and analysis of TPW (blue curve, unit: kg m−2) at 1550 UTC 4 September 2018.

Although a realistic forecast is not expected from the RE87 model since its simulated hurricane does not move as the real Hurricane Florence hurricane did, it is of interest to show the differences between the forecasts made by the RE87 model with two different initial conditions. Figure 13 shows the temporal evolutions of the minimum SLP of the 48-h forecasts for Hurricane Florence made by the RE87 model that is initialized by the 4D-Var VI analysis vortex and the background field (referred to as the control experiment hereafter) at 1530 UTC 4 September 2018. Hurricane Florence experienced a rapid deepening from 1530 UTC 4 September to 1830 UTC 5 September and a weakening from 1830 UTC 5 September to 1530 UTC 6 September 2018. The 4D-Var VI experiment captures reasonably well the variation in intensity during the 24-h period but not the control experiment. If the bogus SLP used in Zou and Xiao (2000) and Park and Zou (2004) is assimilated, the intensity forecast has similar deepening as the 4D-Var VI experiment and the best track during the first 6 h, but gradually weakened while the best track continued deepening in the next 12 h. The bogus SLP experiment performs better than the control but not as good as the 4D-Var VI assimilating satellite retrievals. There are large structural differences in potential temperature anomaly forecasts, liquid water mixing ratios and vertical velocity fields.

Fig. 13.

Temporal evolutions of the minimum sea level pressure (MSLP) forecasts for Hurricane Florence by the RE87 model initialized by the 4D-Var VI vortex (blue curve) and the background field (black curve) at 1530 UTC 4 September 2018. The best track records are also provided (red curve). The results from assimilating bogus SLP are shown with the magenta curve.

Figure 14 presents the axisymmetric structures of the 24-h forecast potential temperature anomalies from the 4D-Var VI experiment, control experiments and observation retrievals at the same time. These cross-sections allow us to assess the capability of RE87 model prediction in simulating the inner-core structures of Florence. The model forecasted warm cores from the 4D-Var VI experiment are located in the upper troposphere within a ± 150 km radial distance from the centre. Warm anomalies greater than 6 K are can be found at ∼10 km in the 24-h forecast. Figure 14b also shows well-defined warm-core anomalies at the centre of the storm between 5 and 12 km only significantly weaker than that seen in the 4D-Var VI experiment (Fig. 14a). When comparing with the warm core structures from the ATMS observations near the 24-h forecast time (Fig. 14c), the warm core structure and intensity in the upper troposphere from the 4D-Var VI experiment is more comparable to the observed features than that of the control experiment. Figure 15 provides a comparison of radial variations of liquid water path (LWP) between the 24-h RE87 model forecast from the 4D-Var VI experiment and the AMSR2 observations valid at 1624 UTC 5 September 2018. The 24-h forecast of LWP compared quite favourably with the AMSR2 retrieval in terms of the maximum values near the storm centre and the decreasing rate of LWP as the radial distance increases. The 24-h forecasted LWP values are slightly smaller than those observed outside the radial distances of 50 km. The LWP values outside the 150 km radial distance from the control experiment are similar with those from the 4D-Var VI, while those near the centre are significantly smaller.

Fig. 14.

Cross-sections of the potential temperature anomalies (unit: K) of the 24-h forecasts from the (a) 4D-Var VI and (b) control experiment. (c) Azimuthally averaged potential temperature anomalies of S-NPP ATMS at 1621 UTC September 5, 2018.

Fig. 15.

Radial variations of LWP from the 24-h RE87 model forecast in 4D-Var VI experiment (solid curve), control experiment (dotted curve), and AMSR2 observations (dashed curve) at 1624 UTC September 6, 2018.

Figure 16 shows two vertical bands of large values of liquid water mixing ratio and vertical velocity located at the eyewall and a minimum in the centre of the storm at the 24-h forecast time in the results of 4D-Var VI and control experiment, respectively. There is a slight increase in the scale of the eye with height. The simulated Florence has an eye varying from less than 50 km near the surface to slightly more than 50 at 12 km. The largest vertical velocities are seen in the 10–km altitude range. In the cross section resulted from the control experiment, nonetheless, the two maximum bands of cloud liquid water mixing ratio at the eyewall are very weak. The largest vertical velocity is found at the centre of the simulated storm and not at the eyewall as seen in the 4D-Var VI experiment. The lack of well-defined vertical and radial structures of the warm core and clouds in the control experiment results from a poor specification of the initial vortex and a potential spin up of the model vortex due to an inconsistency between the ERA-derived background initial vortex and the RE87 model’s dynamics and physics.

Fig. 16.

Cross-sections of liquid water mixing ratio (colour shading, unit: g kg−1) and vertical wind velocity (contours, unit: m s−1) of the 24-h forecasts from the (a) 4D-Var VI and (b) control experiment.

A similar experiment is conducted for Typhoon Mangkhut (2018) (Fig. 17). The analysis time is 1450 UTC 9 September 2018, which was the time when Typhoon Mangkhut started to intensify to its peak strength (dashed line in Fig. 17b). The observations were available at 1454 UTC from NOAA-20 ATMS, 1521 UTC from AMSR2 and 1545 UTC from S-NPP ATMS on September 9, 2018 and are assimilated to generate an initial vortex at 1450 UTC 9 September 2018. The temporal evolution of MSLP in the subsequent 48-h forecast period from the control (black curve) and 4D-Var VI (blue curve) experiments are shown in Fig. 18, which can be compared with the best track records (red curve). Similar with the results in the case of Hurricane Florence, the forecast with the 4D-Var VI initial conditions produced an intensity change closer to the best track than the control forecast with a steady evolution instead of intensification.

Fig. 17.

(a) The best track of Typhoon Mangkhut (2018) with the intensities marked in different shapes and colours during the time period from 5 to 17 September 2018 at a 6-h interval. (b) The mean sea level pressure (black) and 34-kt wind radius (red) from 5 to 17 September 2018. The analysis time of 1450 UTC is marked with the vertical dashed line.

Fig. 18.

The evolution of the minimum sea level pressure in the forecast made with RE87 model in the case of Typhoon Mangkhut taking the 4D-Var initialized vortex (blue curve) as initial conditions and the control experiment (black curve). The MSLP in Best Track records are shown with the red curve.

6.

## Summary

Advances in satellite remote sensing allow the TC structures of many variables (e.g. warm cores, TPW, total column ozone) important in numerical weather prediction to be well resolved. However, the TC vortex initialization does not fully incorporate these satellite retrieval products. The purpose of this study is to present a 4D-Var VI scheme to improve hurricane initialization by assimilating satellite retrieval products. The proposed 4D-Var VI scheme allows the axisymmetric components contained in the fields retrievable from satellite remote sensing observations to be assimilated into a nonhydrostatic axisymmetric numerical model in which convection is accounted for, i.e. the RE87 model. The axisymmetric hurricane prediction model serves as a strong constraint to spin up other fields not readily assimilated. The 4D-Var VI scheme has been developed and tested on the initialization and the 48-h prediction of Hurricane Florence (2018) and Typhoon Mangkhut (2018) during their deepening stages. We show that the 4D-Var VI scheme assimilating satellite microwave-retrieved TPW and warm-core retrievals is effective in generating a consistent three-dimensional structure of the initial vortex without having to use extremely simple dynamics to derive other variables from the assimilated variables. By assimilating the azimuthally averaged TWP and warm-core anomaly in 4D space with the hurricane forecast model and its adjoint model that carries the information forward and backward in time, a model state can be generated that contains a well-defined warm core in the upper levels near the centre, a strong upward motion and high cloud liquid water content at the eyewall, and a wet core in the low levels near the centre of the initial and forecasted vortices. Because the initial model state obtained by the 4D-Var VI was well adapted to the forecast model, dramatic improvements are observed in the intensity forecasts as well as in the description of the inner structures of the predicted storm during a 24-h forecasting period. Due to the use of a short time window (slightly longer than one hour), the proposed 4D-Var VI scheme, though carried out in 4D space, is computationally much cheaper than a 4D-Var experiment involving a much more complex nonaxisymmetric model.

Applications of the proposed 4D-Var scheme to a sufficiently large statistical sample of hurricanes are still needed to test its general promise to improve hurricane intensity forecasts. In this regard, assessing the sensitivity of the 4D-Var results to the selection of the assimilated satellite retrieval products and the storm axisymmetry will first be done. A final intended application is to incorporate the 4D-Var VI system developed in this study into realistic models for vortex initialization, such as the Hurricane Weather Research and Forecasting system or the Advanced Research Weather Research and Forecasting model, or both.