GSMaP RIKEN Nowcast: Global Precipitation Nowcasting with Data Assimilation

Since January 2016, RIKEN has run an extrapolation-based nowcasting system of global precipitation in real time. Although our previous paper reported the effectiveness of using data assimilation in a limited verification period, the long-term stability of the forecast accuracy through different seasons has not been investigated. In addition, the algorithm was updated seven times between January 2016 and March 2018. Therefore, this paper aims to examine how motion vectors can be derived more accurately, and how data assimilation can stably constrain an advection-diffusion model for extrapolation for the long-term operation. The Japan Aerospace Exploration Agency’s Global Satellite Mapping of Precipitation (GSMaP) near-real-time product is the only input to the nowcasting system. The motion vectors of precipitation areas are computed by a cross-correlation method, and the Local Ensemble Transform Kalman Filter is used to generate a smooth, complete set of motion vectors. Precipitation areas are extrapolated in time up to 12 hours ahead, and the product, called GSMaP RIKEN Nowcast, is disseminated on a webpage in real time. Most of the algorithmic updates involved improving the estimation of the motion vectors, and the forecast accuracy was gradually and consistently improved by these updates. In particular, the threat scores increased the most at approximately 40°S and 40°N. A performance decrease in the northern hemisphere winter was also reduced by reducing noise in advection. The time series of the ensemble spread demonstrated that an increase in the number of available motion vectors by a system update led to a decrease in the ensemble spread, and vice versa. Corresponding author: Shigenori Otsuka, RIKEN Center for Computational Science, 7-1-26 Minatojima-minami-machi, Chuo-ku, Kobe, Hyogo 650-0047, Japan E-mail: shigenori.otsuka@riken.jp J-stage Advance Published Date: 3 September 2019 Journal of the Meteorological Society of Japan Vol. 97, No. 6 1100


Introduction
The short-term prediction of precipitation is of great interest in various fields, such as disaster prevention, agriculture, the economy, and daily life. Ground-based precipitation observations, such as weather radar and gauge networks, provide accurate and timely information. However, sparsely observed areas, such as developing countries and ocean basins, suffer from water-related disasters partially due to the lack of information. Satellite observations are a powerful tool over these areas. Satellite-borne precipitation radars, such as the Tropical Rainfall Measuring Mission (TRMM) Precipitation Radar (PR) (Kozu et al. 2001;Iguchi et al. 2009) and the Global Precipitation Measurement (GPM) Dual-frequency Precipitation Radar (DPR) (Kojima et al. 2012;Seto and Iguchi 2015), are useful for global quantitative precipitation estimation; however, their spatiotemporal coverage is limited. Therefore, passive microwave observations and infrared observations are the main sources of near-realtime precipitation information, and passive microwave rainfall algorithms have been improved using TRMM PR and GPM DPR.
Using microwave observations from the GPM Micro wave Imager (GMI) on the GPM core satellite and other microwave sensors of the GPM constellation partners, the following near-real-time satellite-based products are distributed: the Japan Aerospace Exploration Agency (JAXA) Global Satellite Mapping of Precipitation (GSMaP; Kubota et al. 2007;) and the National Aeronautics and Space Administration (NASA) Integrated Multi-satellite Retrievals for GPM (IMERG; Huffman et al. 2018).
For real applications, short-term predictions are also useful. Ground-based radar precipitation nowcasting has a long history (e.g., Rinehart and Garvey 1978;Dixon and Wiener 1993;Bowler et al. 2006), whereas satellite-based global precipitation nowcasting is a relatively new field. Recently, Otsuka et al. (2016a) developed a global precipitation nowcasting system using GSMaP, and the system has run at RIKEN in real time since January 2016 (GSMaP RIKEN Nowcast, hereinafter referred to as RNC). The predictions are made open to the public with the permission of the Japan Meteorological Agency (JMA), as required by Japanese law. Although Otsuka et al. (2016a) presented initial verification results for a particular month, the long-term stability of RNC in different seasons has not yet been investigated. Our real-time operations spanning more than two years have revealed its advantages and limitations on longer time scales. The knowledge gained from the real-time operations motivated us to update the system to improve its accuracy and computational performance.
This paper presents the verification results for the period from January 2016 to March 2018 including the impact of each algorithm update. The key questions are as follows: 1) how motion vectors for extrapolation can be derived more accurately, and 2) how an advection-diffusion model for extrapolation can be stably constrained by data assimilation. The first question will contribute to the better implementation of motion tracking techniques, whereas the latter will lead to the better handling of uncertainties in data assimilation. The paper is organized as follows. Section 2 introduces the system design, Section 3 presents the verification results, Section 4 provides a discussion, and Section 5 summarizes the paper. System details are described in the Appendix.

Overview
The workflow of the system is as follows (Fig. 1): 1. Data transfer of GSMaP from JAXA to RIKEN 2. One-hour ensemble forecast of motion vectors u from time t -1 to t : u t-1 a, m ® u t f , m , where superscripts f , a, and m represent the forecast, analysis, and ensemble member, respectively 3. Motion vector computation (u o t ) using GSMaP Near-Real-Time (NRT) images at time t -1 and t (denoted NRT t-1 and NRT t , respectively) 4. Ensemble Kalman filter using u t f , m and u o t to obtain analysis ensemble-mean motion vectors u t a and analysis ensemble u t a, m 5. Deterministic 12-hour forecast of precipitation from time t to t + 12 (RNC t t+12 ) using u t a and NRT t as the initial condition 6. Visualization 7. Visual inspection by a licensed forecaster Keywords nowcast; precipitation; data assimilation Citation Otsuka, S., S. Kotsuki, M. Ohhigashi, and T. Miyoshi, 2019: GSMaP RIKEN Nowcast: Global precipitation nowcasting with data assimilation. J. Meteor. Soc. Japan, 97, 1099-1117, doi:10.2151/jmsj.2019 8. Dissemination at RIKEN's webpage and JAXA's webpage Running an operational system helps identifying key problems for future research. This includes an end-to-end workflow, such as data transfer, computational cost, visualization, and dissemination to end users. For this purpose, GSMaP RNC forecasts are made readily available to the public in real time at https://weather.riken.jp/ as a part of our research activities. See Appendix A for more information.

Data and algorithm
Here, we provide a brief summary of the system. See Otsuka et al. (2016a, b) for details. The input to the system is the GSMaP NRT algorithm version 6 Kachi et al. 2011) rain rate at a resolution of 0.1° × 0.1° from 60°S to 60°N. NRT is available every hour, after four hours from the end of the one-hour observation time window. The motion vectors of the NRT rain fields are computed by a variant of the Tracking Radar Echoes by Correlation (TREC) method (Rinehart and Garvey 1978).
To obtain a smooth, complete set of motion vectors from noisy TREC-based motion vectors, the Local Ensemble Transform Kalman Filter (LETKF; Hunt et al. 2007;Miyoshi and Yamane 2007) is used with 20 ensemble members. The covariance localization function is a Gaussian with a length scale of 100 km, while the covariance inflation method is a multiplicative inflation with a factor of 1.01 simultaneously with the relaxation to prior spread (Whitaker and Hamill 2012) so that the ensemble spread does not change in time (Otsuka et al. 2016a).
The prediction model is an advection equation with a divergence damping term, while the time integration method is the second-order Adams-Bashforth scheme. The spatial discretization method is the fifth-order Weighted Essentially Non-Oscillatory scheme (Liu et al. 1994) for rain, and the first-order upwind scheme for the motion vectors. An equirectangular projection is used throughout the system, so that the input data are simply considered plain bitmap images, while the motion vectors are treated in the units of pixels per hour, or 0.1° h −1 .
A TREC-based method is used to obtain the motion vectors; this type of motion-tracking technique is also known as cross-correlation or optical flow. The principle of this technique is to identify a spatial shift in a local pattern (so-called template) between two or more consecutive time steps by performing spatial correlation under the assumption that the pattern shape does not change in time. Of all possible location shifts, the one with the highest correlation is selected as the estimated motion. This procedure is called template matching, and the search space is called the cross-correlation surface. In this system, several changes are made to the original TREC: fractional motion, circular template, treatment of missing pixels, and quality control (Otsuka et al. 2016a, b).
To improve the accuracy of motion vectors for extrapolation, the system had seven major updates to its algorithms and parameters during the period from January 2016 to March 2018 (Table 1). These updates aimed to improve the quality of the motion vectors (V4, V10, V23), increase the number of available motion vectors (V4, V11, V23, V25), and stabilize the system (V1, V4a, V23, V25). The key goals of each update are as follows: • In the version presented in Otsuka et al. (2016a), the advection occasionally exhibited numerical noise.
To prevent this, V1 shortened the time interval of advection. • V1 and earlier versions did not estimate motion vectors over weak precipitation areas; this degraded the performance over those areas due to the lack of information. V4 lowered the threshold of rainy pixels in TREC to compute the motion vectors of weak precipitation. • NRT initially assigns large negative values to missing pixels, such as land covered by snow (Seto et al. 2005 The details of each update are provided in Appendix B, and the impact of these updates is discussed in Section 3.

Results
First, Sections 3.1 -3.2 describe the overall performance of the system. Next, Section 3.3 provides quantitative comparisons between different versions at the time of each update; this is primarily aimed to evaluate how motion vectors for extrapolation can be derived more accurately. Finally, Section 3.4 provides analyses of how an advection-diffusion equation can be stably constrained by data assimilation.
The accuracy of the RNC precipitation forecast is verified against the GSMaP standard product (moving vector with Kalman filter, MVK) version 6 because MVK is considered more reliable than NRT; MVK uses more satellite data and considers the deformation of rain areas and the temporal change of rain rate in a more sophisticated manner ). The accuracy of the tracking algorithm is also verified using NRT, because the motion vectors are designed to represent the motion of the NRT precipitation pattern. The accuracy of the RNC is compared with that of the NRT-based Eulerian persistence forecast (PER), which is considered the baseline. Although the NRT algorithm is occasionally updated, both RNC and PER commonly include the effects of these updates; thus, the difference between RNC and PER is supposed to represent the effect of space-time extrapolation. This paper uses the operational version of RNC; creating a reanalysis dataset with the latest algorithm is beyond the scope of this paper. Figure 2 presents an example of a nowcast initialized at 0000 UTC 23 March 2017 by the latest version (V25). The color shading represents hit (green), miss (red), and false alarm (blue) of a rain rate greater than 0.1 mm h −1 with respect to MVK. At a forecast time (FT) of 0 h, the initial condition is identical to NRT, that is, most of the rainy pixels greater than 0.1 mm h −1 are green in Fig. 2a. Red and blue pixels also exist due to the differences between NRT and MVK. As the FT increases, blue and red pixels increase. At FT = 12 h (Fig. 2g), red areas appear to the east of blue areas at higher latitudes, indicating that the eastward motion of extratropical cyclones tends to be underestimated. This slow-bias problem is discussed further in later sections.

Visual analysis of spatial patterns
A comparison between FT = 6 and 12 h indicates that rapid changes in the rain areas are one of the reasons for the underestimation of eastward motion (Figs. 2e, g); the rain area south of Australia quickly forms a new rain area to the east, while another T-bone rain area in the southern Pacific displays a rapid shift to the southeast. These changes appear to be linearly unpredictable.

Statistical verification
Before investigating the temporal changes in the forecast accuracy, time-averaged performance is examined as a function of forecast time to provide fundamental information on the forecast skill. Figure  3 presents the mean threat scores (TS) against the MVK for threshold values of 0.1, 1, and 5 mm h −1 as a function of FT for RNC and PER, as well as RNC TS minus PER TS (ΔTS), computed for the northern hemisphere (20 -60°N), the tropics (20°S -20°N), the southern hemisphere (20 -60°S), and the entire computational domain. In all verification regions, TS is the highest with a threshold of 0.1 mm h −1 and lowest with a threshold of 5 mm h −1 . The TS at FT = 0 h is smaller than unity, which reflects the difference between NRT and MVK. In particular, the TS for 5 mm h −1 at FT = 0 h is less than 0.6, while at the 0.1 and 1 mm h −1 thresholds, TSs are in the range of 0.71 to 0.78; this signifies that the differences between NRT and MVK precipitation are less noticeable for lower rain thresholds.
The threat score difference, ΔTS, is the largest for 0.1 mm h −1 and the smallest for 5 mm h −1 . Although the TS is the highest over 20°S -20°N (Fig. 3b), the ΔTS is the smallest over 20°S -20°N, indicating that linear extrapolation in the tropics is not as advantageous against PER as in the extratropics, presumably due to the slower motion of large-scale precipitation systems, such as the Madden-Julian oscillation. It should be noted that the lifecycle of subgrid-scale convective cells does not affect the accuracy of spacetime extrapolation on a spatiotemporal scale of 0.1° × 0.1° and 1 hour, even if convective systems play a dominant role in the tropics. The analysis regions 20 -60°N and 20 -60°S have similar performance to each other (Figs. 3a, c). The ΔTS for 0.1 mm h −1 (red dashed line) maximizes at FT = 5 h over 20 -60°N and 20 -60°S, whereas that over 20°S -20°N maximizes at FT = 3 h. These results indicate that the advantage of advection in weak rain areas lasts up to 12 hours, and lasts longer in the extratropics than in the tropics. The ΔTS for 1 mm h −1 (green dashed line) maximizes at FT = 3 -4 h, and ΔTS for 5 mm h −1 (blue dashed line) always maximizes at FT = 1 h. The prediction skill for intense rain areas quickly decreases as FT increases. The general tendency is similar to the results in Fig. 5 of Otsuka et al. (2016a).
Next, the temporal change of the forecast skill in response to the system updates is examined. Figure 4 presents the time series of the monthly-mean global TS and ΔTS. The performance of RNC for 0.1 mm h −1 improved during this period for the entire forecast time, whereas performance of PER did not exhibit this trend (Fig. 4a). The time series of ΔTS for 0.1 mm h −1 (Fig. 4d) has an annual cycle with two negative spikes around February and December 2016; details are discussed in Fig. 6. The time series of ΔTS for 0.1 mm h −1 also indicates that the advantage of space-time extrapolation rapidly increases between FT = 1 and 2 h.
The performance of RNC for 1 mm h −1 (Fig. 4b) is more stable than that for 0.1 mm h −1 during this experimental period. The difference between the nowcast and persistence (ΔTS in Fig. 4e) displays a similar temporal change as in ΔTS for 0.1 mm h −1 (Fig. 4d); however, the amplitude is much smaller. The advantage of nowcasting maximizes at approximately FT = 2 -3 h for 1 mm h −1 , which is consistent with the green dashed line in Fig. 3d. The performance of RNC for 5 mm h −1 is also stable during the experimental period (Fig. 4c), and the difference between RNC and PER for 5 mm h −1 is almost identical between FT = 1 -3 h (Fig. 4f).
Although the scores for 1 and 5 mm h −1 imply that the changes in the system do not significantly improve the forecast performance, the TSs computed with respect to NRT do improve during this experimental period (Fig. 5). Figure 5 illustrates the performance at FT = 1 h to highlight a time at which advection is successful as a forecasting approach; longer forecasts are affected more strongly by the evolution of precipitation systems. Therefore, the changes in the system help better capture the motion of NRT precipitation. However, uncertainties in both NRT and MVK have a larger impact on the TSs than the improvements in the motion vectors, as improving forecast performance also requires improvements in the quality of the input.
As the precipitation rate increases, the discrepancy between NRT and MVK increases (Fig. 5c). Figure 6 is similar to Figs. 4d -f, but is verified in each region. The ΔTS for 0.1 mm h −1 over 20 -60°S (Fig. 6g) displays the greatest improvement in performance; the ΔTS increases by approximately 0.02 -0.03 during this period. Another uniqueness in Fig. 6g is an increase in the ΔTS with the forecast time. This may be related to the faster "apparent" zonal motion of synoptic weather systems in the units of pixels per hour over 40 -60°S on the equirectangular projection. Faster apparent motion is more advantageous to the space-time extrapolation than Eulerian persistence. In the northern hemisphere, in contrast, the north Pacific storm track region lies at slightly lower latitudes than the storm track region in the southern hemisphere (e.g., Catto et al. 2012), and this effect is weakened. Apparent zonal motion at 35°N on the equirectangular projection is 14 % smaller than that at 45°S if the actual zonal motion is identical.
The ΔTS for 0.1 mm h −1 over 20 -60°N displays a large decrease around February and December 2016 (Fig. 6a), whereas other regions or other thresholds do not exhibit such a decrease (Figs. 6b -i). This is caused by numerical noise of advection along the boundaries between valid and missing pixels over Eurasia (Appendix B.3). The annual cycle of decreased performance appears because land with snow cover leads to missing pixels in GSMaP (Seto et al. 2005). The problem was addressed by the update on 20 January 2017 (Appendix B.3), and the boreal winter in 2017 -18 did not exhibit a decrease in performance.

Impact of each system update
Rerunning each version of the operational system for the entire period is unrealistic; instead, it is normal practice to run old and new systems in parallel, (i.e., parallel run) at the time of each system update. This process ensures that the new system runs properly in preparation for the actual operational update in place. Figure 7 presents the TS improvements of RNC against MVK version 6 as a function of forecast time and latitude due to the following system updates: V1 to V4 (Figs. 7a -   provided due to the misconfiguration of the covariance inflation (Appendix B.8). It should be noted that the impact of each parameter change or new method on forecast accuracy is assessed before each system update; however, for simplicity, differences between the major versions are presented in this paper.

a. V1 to V4
This update displayed the largest improvement during the period from January 2016 to March 2018 (Figs. 7a -c), and the most significant improvements appeared at 40°S and 40°N at FT = 5 h for 0.1 mm h −1 . These changes are consistent with the algorithmic update; computing additional motion vectors in the weak rain areas expands the "observed" areas ( Fig. 8a), leading to the highest improvement in the weak rain areas (Appendix B.2). Figure 9a demonstrates that this update accelerates eastward motion in the extratropics and westward motion in the tropics, which is effective in reducing the slow bias of the motion.

b. V4a to V10
The performance at FT = 1 -3 h generally improved at all latitudes, which is consistent with the update of the motion vector quality control. As in the previous update, the forecasts improved the most in the extratropics at approximately FT = 6 h for 0.1 and 1 mm h −1 (Figs. 7d -f). In the tropics, the TS improved up to FT = 3 h; however, some degradation appeared over a longer forecast time. This may be due to the nature of tropical convective systems; precisely tracking an existing convective system does not necessarily improve longer forecasts beyond the life cycle of that convection. As a result of stricter quality control, the number of available motion vectors decreased (Fig.  8b). Figure 9b indicates that the motion vectors that passed quality control did not change the mean value before and after the update.

c. V10 to V11
The TS slightly improved at the highest latitudes for 0.1 and 1 mm h −1 (Figs. 7g, h), which is consistent   Table 2. The thresholds are 0.5 (left), 1 (center), and 5 mm h −1 (right). Data without hatching are statistically significant with a 95 % confidence interval estimated from 10,000 bootstrap simulations.
with the increase of available motion vectors near the lateral boundaries (Fig. 8c, Appendix B.5). The motion vectors that existed before this update outside the boundary regions were not changed by this update (Fig. 9c). For 5 mm h −1 , the TS difference was not always positive near 60°S (Fig. 7i).

d. V11 to V23
The forecasts improved the most in the extratropics at a longer forecast time (Figs. 7j -l). The number of available motion vectors also increased (Fig. 8d), and the eastward component of the motion vectors became stronger (Fig. 9d). The accelerated eastward motion appeared to reduce the slow bias; however, the TS for 0.1 mm h −1 slightly degraded in the tropics, where the number of available motion vectors decreased (Fig.  8d). Again, the improvements are more evident for a weaker precipitation rate. This is consistent with the tuning of the spatial scale at high latitudes, and the performance change in the tropics can be considered a side effect of the changes in the quality control methods (Appendix B.6).  Table 2. e. V23 to V25 The forecasts improved the most between 40 -60°S, and a secondary peak appeared near 30 -40°N. The eastward component of the motion vectors at these latitudes became stronger (Fig. 9e), further reducing the slow bias of the extratropical cyclones (Fig. 2). These latitudes correspond to the population of cold and warm fronts (e.g., Catto et al. 2012). This change is consistent with the introduction of a latitude-dependent search space in TREC (Appendix B.7). The improvement is largest for 0.1 mm h −1 and smallest for 5 mm h −1 . No difference appears in the tropics.

Analysis on the assimilation system
The ensemble Kalman filter finds an optimum solution based on the forecast ensemble covariance and observation error covariance. The observation error covariance is given a priori, whereas the forecast ensemble covariance is controlled by the internal dynamics and the covariance inflation method. In a numerical weather prediction (NWP) system, it is common that internal dynamics plays a central role. However, the current system adopts an advectiondiffusion equation, and the internal dynamics lacks mechanisms of error growth that are essential in the real atmosphere. Therefore, the ensemble spread of motion vectors in the analysis cycle must be carefully examined.
First, an instantaneous distribution of the ensemble spread is examined to understand the uncertainties that the advection-diffusion model represents. Figure  10 presents the horizontal distributions of rain, the zonal component of the analysis motion vectors u a , the ensemble spread of the zonal component of the first guess motion vectors, and the corresponding analysis increments u au f at 0000 UTC 25 December 2017. An elongated area of large ensemble spread appears near 130 -150°E, 20 -30°N (Fig. 10c), which corresponds to a "frontal" region to the southeast of a rain area centered at 147°E, 40°N (Figs. 10a, b). This also coincides with a large increment in Fig. 10d. Uncertainty in the location of the shear line  Table 2. produces a large ensemble spread, and "observations" of motion vectors effectively fix the location based on the background error covariance structure. It should be noted that this frontal structure does not necessarily correspond to the frontal system in the atmospheric dynamics; the only thing the nowcasting system handles is the apparent motion of rain patterns. An area of large ensemble spread to the north of the rain area is produced by the covariance inflation method (Fig. 10c). In ordinary NWP systems, the forecast error grows following the internal dynamics; however, the current extrapolation system lacks such internal dynamics and requires stronger covariance inflation. If not inflated artificially, the ensemble spread continuously decreases due to the dissipative nature of the advection-diffusion equation. If a new precipitation system appears in the high-spread area, the motion vectors of the system are quickly assimilated. More precisely, two consecutive images of a new precipitation system are required to compute the motion vectors of that precipitation system. A new precipitation system is usually small when it initially appears and moves in the direction of pre-existing motion vectors.
As an operational system, the ensemble spread must stay within a reasonable range, because the assimilation system will eventually blow up if the spread continues to increase or decrease. Figure 11 presents the time series of the ensemble spread of analysis motion vectors. The ensemble spread was unstable in V1 (black line), whereas it was mostly stabilized from V4 (red line). The increase in the number of available motion vectors in V4 (Fig. 8a) appeared to decrease the spread from V1 to V4. A small spread decrease in July 2016 appeared to correspond to a bug fix in LETKF. Another decrease in the meridional component from December 2016 to February 2017 was caused by the system misconfiguration. An increase in spread from V4a to V10 (green line) is consistent with the decrease in the number of available motion vectors (Fig. 8b) due to stricter quality control criteria in V10. A spread decrease from V10 to V11 (blue line) is consistent with the increase in the number of available motion vectors near the south/north lateral boundaries (Fig. 8c). From V11 to V23, the spread was slightly decreased, and from V23 to V25, the spread did not significantly change.

Discussion
The forecast error of space-time extrapolation stems from multiple sources: tracking error, advection error, and the error from the Lagrangian persistence assumption. The first two errors are the main themes of this paper, whereas the last error is left for future work. The system updates had clear contributions to the increased number of motion vectors (Fig. 8), to the increased quality of motion vectors (Fig. 9), and to the improved accuracy of forecast by advection (Fig. 7). The system strongly relies on the input data because the model consists of the advection and diffusion terms only; to constrain the behavior of the system, information must be effectively extracted from noisy observations. As mentioned in Section 3, the eastward motion of extratropical cyclones tended to be underestimated. This was partly caused by the suboptimality of parameters, such as the patch diameter and search space in TREC at high latitudes. This problem was mitigated by latitudinally varying these parameters (Table 1, 9,. However, there remains a slow bias of motion vectors in the case of extratropical cyclones (Fig. 2); further improvements in both motion tracking and data assimilation are thus required. Another problem in detecting motion vectors is the spatiotemporal continuity of this input data. For example, small red dots appear and disappear suddenly to the west of Australia in Figs. 2c, 2e, and 2g. In this case, it is difficult to track the motion of these areas. These sudden changes may be due to the satellite-to-satellite differences between sensors used in GSMaP, as well as the algorithm used to produce the mosaic product. In the present study, data assimilation provides reasonable estimates over these areas. However, cloud tracking by geostationary meteorological satellites may be less influenced by the discontinuity of data.
In principle, space-time extrapolation cannot capture the lifecycle of extratropical cyclones. For example, a comma-shaped rain area in the northwestern Pacific (Fig. 2b) changes to a linear rain band (Fig. 2h) in 12 hours. NWP can predict the evolution of precipitation systems within the limit of predictability based on the physics, and machine learning techniques may also be able to predict the evolution based on past data.

Summary
This paper presented an overview of GSMaP RNC, a GSMaP-based global precipitation nowcasting system using data assimilation. The main foci were on 1) how to better estimate motion vectors, and 2) how to stably constrain the advection-diffusion model with data assimilation.
The system has been stably running in real time since January 2016 for three and a half years as of July 2019, and seven major system updates have been applied (Table 1). Most of the updates involved better estimating the motion vectors. The forecasts are disseminated via RIKEN's and JAXA's webpages every hour in real time under a weather forecasting license.
Verification against GSMaP MVK demonstrated that the space-time extrapolation functioned properly . A comparison with Eulerian persistence forecasts revealed that the advantage of space-time extrapolation gradually increased from January 2016 to March 2018 as a result of seven major system updates Fig. 11. Time series of spatially-averaged analysis ensemble spread of (a) zonal and (b) meridional components of motion vectors from January 2016 to April 2018. V1 (black), V4 -V4a (red), V10 (green), V11 (blue), V23 (purple), and V25 (cyan). . In particular, the updates ameliorated the decreased winter performance in the northern hemisphere. The long-term operation revealed this seasonal performance drop and led to the improvement of the system. Comparisons of TSs before and after each system update demonstrated that the scores improved the most at around 40°S and 40°N (Fig. 7). Light precipitation areas improved significantly, whereas heavy precipitation areas did not display as much improvement.
Analysis on the ensemble spread demonstrated that the increase and decrease of the number of available motion vectors changed the ensemble spread between each system update (Fig. 11). The motion vectors were effectively modified over the regions in which the ensemble spread increased, such as areas near the frontal systems (Fig. 10). Better quantifying the uncertainties of the motion vectors in TREC estimations and in the first guess are the next steps for further improving the stability of the system as well as the forecast accuracy.
The space-time extrapolation method has an apparent limit of forecast accuracy, particularly for strong rain areas (Fig. 3). It is natural to run an NWP system for longer forecasts beyond this limit; for example, Kotsuki et al. (2017) developed a global NWP system that assimilates GSMaP to better predict precipitation. In regional weather prediction, blending extrapolation and NWP is widely practiced to take advantage of each system (e.g., Golding 1998;Sun et al. 2014). Blending is also beneficial in global precipitation forecasting (Kotsuki et al. 2019).
Probabilistic forecasting appears to be an essential tool in practical precipitation nowcasting, and ensemble-based methods have been widely used (e.g., Germann and Zawadzki 2004;Bowler et al. 2006). Our system uses an ensemble Kalman filter for motion vector construction (Fig. 1), and the "observed" motion vectors are effectively assimilated using a background error covariance structure (Fig. 10). Although the current GSMaP RNC runs a deterministic precipitation forecast using ensemble mean motion vectors, the system can be easily extended to an ensemble forecast.
Another direction for future work may be the use of data-driven approaches in nowcasting. One such approach is machine learning; for example, Han et al. (2017) used a support vector machine for precipitation nowcasting. In addition, as a quantitative precipitation estimation algorithm, Hirose et al. (2019) adopted the random forest algorithm. Another example is the analog forecasting (Lorenz 1969), which has applications in precipitation nowcasting (e.g., Obled et al. 2002). With the help of increased computing power and data size, these existing techniques have the potential to contribute to precipitation nowcasting. mm h −1 to R ³ 0 mm h −1 . This increased the number of available motion vectors by approximately four times (Fig. 8a).
• The quality control criterion of TREC motion vectors was changed. TREC estimates motion vectors by template matching, in which similarity is measured by the Pearson correlation coefficient, r. Coefficient r becomes much smaller than 1 if TREC fails to find a similar pattern. In this update, the criterion of acceptance was changed from r > 0 to r > 1 ´ 10 −5 . Stricter quality control was thus possible as a result of the increased number of TREC motion vectors.

B.3 Update on 20 January 2017 (V4a)
• Missing pixels are initially assigned by large negative values in NRT. However, sudden changes in pixel values between valid rain rates and missing values lead to numerical noise in advection. To suppress this noise, the missing values in NRT are replaced by 0 mm h −1 at the beginning of a 12-hour forecast.

B.4 Updates on 10 February 2017 (V10)
Although data assimilation reduces noise in the motion vector field, the quality of the "observed" motion is still crucial. Thus, V10 introduced the following three changes: • The quality control criterion of the correlation coefficient in TREC was changed from r > 1 ´ 10 −5 to r > 0.2. • The quality of a motion vector depends not only on r, but also on the clarity of the peak of r on a cross-correlation surface in TREC.  (Fig. 9d). The following three changes were made to introduce the latitude-dependent d CC . • The criterion for the variable range of r in TREC was changed from Δr ³ 0.4 to a normalized version Δr N ³ 0.8 (Appendix C.1). The quality control criterion of the correlation coefficient was also changed to r > 0.6. • Rain fields were smoothed only for TREC computation by a box-mean method with a box size of 0.1 d CC ´ 0.1 d CC (Appendix C.3). • The minimum ratio of valid pixels within the template in TREC was changed from 5 % to 1.25 %.

B.7 Updates on 1 March 2018 (V25)
• Although V23 improved forecast accuracy at high latitudes, there remained room for latitude-dependent optimization. Cross-correlation was performed within a finite search space S defined by | (δi, δ j ) | £ L, where (δi, δ j ) represents the offset in the zonal and meridional directions, and L = 0.6°. To take into account the length of the latitudinal circle, S was changed to | (δi cos θ, δ j ) | £ L, where θ is the latitude of the analysis point. This accelerated the eastward component of the motion vectors at high latitudes (Fig. 9e) and slightly increased the number of available motion vectors at high latitudes in the southern hemisphere (indistinguishable in Fig. 8e). • The time interval of advection Δt was changed from 60 s to 30 s for numerical stability. This change was due to the introduction of latitude-dependent S described above; the maximum zonal motion at high latitudes increased.

B.8 Other system changes
There were several additional changes to the system. The computational performance was improved several times during this period. A bug fix of LETKF applied at 0100 UTC 5 July 2016 resulted in a jump in the time series of the ensemble spread.
MPI parallelization was introduced at 0400 UTC 16 October 2016. Another bug fix in LETKF was applied at 2200 UTC 6 December 2016. Covariance inflation for the north-south motion variable in LETKF was unintentionally not applied in the assimilation cycles from 1200 UTC 27 December 2016 to 2100 UTC 8 February 2017. This was caused by a system misconfiguration during the system update at 1200 UTC 27 December 2016. , ι (1) where d CC denotes the diameter of the circular template for TREC, t denotes time, and subscript n denotes the time index. If min (r (i, d CC , S, t n-1 , t n-1 )) < 0, Δr is used instead of Δr N . In the case of Fig. C1f, Δr = 0.27 and Δr N = 1.4. From 29 November 2017 onwards, motion vectors with Δr N < 0.8 were rejected. Figure C2a presents an example of erroneous motion vectors. Most of the motion vectors represent the motion of a strong rain area moving to the right (shaded red), whereas the edge of the motion vector field represents the motion of weak rain area moving to the left (blue). The discontinuity between the higher rain rate at the center and the surrounding lower rain rate produces this ring-shaped discontinuity in the motion vector field. This pattern frequently appears in real data due to noise in the rain field. To avoid this "ring effect", the following procedure was introduced: if all the rainy pixels at time t -1 with rain rate R greater than threshold R 0 are located between the prescribed distances d 0 /2 and d CC /2 from the center of the template, the motion at that grid point is not computed. Diameter d 0 controls the ring effect, where d 0 ranges between 0 and d CC . Using this method, the appearance frequency of this pattern decreases (Fig.  The left column displays the input fields ( f (x)), while the right column displays the corresponding cross-correlation surfaces (red) and "predicted" cross-correlation surfaces (black) as functions of the location shift. Note that the x axis in (e) is different from that in (a) and (c). The template size is the same as the computational domain.

C.3 Smoothing rain field in TREC
Since 29 November 2017, rain fields have been spatially smoothed during the computation of cross-correlation. Raw rain fields include structures on scales much smaller than the template size, d CC . If the template size d CC increases, it becomes increasingly difficult to determine the global maximum on a cross-correlation surface due to multiple local maxima that appear as a result of the scale gap. For example, small-scale features may move differently from largescale features. Similarly, some portions of the template may move differently from other portions. To address this problem, rain fields are smoothed by two-dimensional box averaging; the box size is specified as 0.1d CC ´ 0.1d CC . This modification was introduced in association with the latitude-dependent template size.