A Feasibility Study on the High-Resolution Regional Reanalysis over Japan Assimilating Only Conventional Observations as an Alternative to the Dynamical Downscaling

The feasibility of regional reanalysis assimilating only conventional observations was investigated as an alternative to dynamical downscaling to estimate the past three-dimensional high-resolution atmospheric fields with long-term homogeneity over about 60 years. The two types of widely applied dynamical downscaling approaches have problems. One, with a serial long-term time-integration, often fails to reproduce synoptic-scale systems and precipitation patterns. The other, with frequent reinitializations, underestimates precipitation due to insufficient spin-up. To address these problems maintaining long-term homogeneity, we proposed the regional reanalysis assimilating only the conventional observations. We examined it by paying special attention to summer precipitation, through one-month experiment before conducting a long-term reanalysis. The system was designed to assimilate surface pressure and radiosonde upper-air observations using the Japan Meteorological Agency’s nonhydrostatic model (NHM) and the local ensemble transform Kalman filter (LETKF). It covered Japan and its surrounding area with a 5-km grid spacing and East Asia with a 25-km grid spacing, applying one-way double nesting in the Japanese 55-year reanalysis (JRA-55). The regional reanalysis overcame the problems with both types of dynamical downscaling approaches. It reproduced actual synoptic-scale systems and precipitation patterns better. It also realistically described spatial variability and precipitation intensity. The 5-km grid spacing regional reanalysis reproduced frequency of heavy precipiJournal of the Meteorological Society of Japan Vol. 96, No. 6 566


Introduction
Multi-decade high-resolution atmospheric datasets with long-term consistency in quality are required for various purposes, such as the investigation of local responses to climate change, the study of past extreme mesoscale events, and the advancement of application of meteorological information in supporting decision making.
A technique widely applied to production of highresolution atmospheric data is dynamical downscaling, since attempted by Giorgi and Bates (1989) (e.g., Giorgi and Gutowski, Jr. 2015). It provides highresolution atmospheric data with physical consistency by driving a physical-based regional climate model forced by the lateral boundary conditions of the model, which are obtained from coarser atmospheric data such as global reanalyses. There are two types of dynamical downscaling approaches: one with serial time integration over a long term and the other with frequent reinitilizations. However, both types of dynamical downscaling approaches have some problems. The dynamical downscaling with serial time integration often fails to simulate the actual states of the atmosphere, including synoptic-scale systems that are appropriately described in the outer data (Kayaba et al. 2016). This fact implies that the constraints from the model dynamics with boundary conditions are insufficient to accurately determine the inner fields in the long-term simulations. In contrast, the dynamical downscaling approach with frequent reinitializations is sometimes adopted to overcome the problem with the dynamical downscaling approach with serial longterm time integration (Qian et al. 2003;Lucas-Picher et al. 2013;Kayaba et al. 2016). In this type of dynamical downscaling approach, the fields are frequently reinitialized to perform cold-start short-term time integrations. Although it successfully improves the description of synoptic-scale fields, this type of dynami-cal downscaling approach suffers from other problems of insufficient spin-up and temporal discontinuity at the joints of the downscaling cycle. Kayaba et al. (2016) completed dynamical downscaling of JRA-55 with this approach over Japan and its surrounding area with a 5-km grid spacing from 1958 to 2012 (DSJRA-55) and pointed out that DSJRA-55 underestimates the frequency of heavy precipitations. The underestimation implies that DSJRA-55 is affected by insufficient spin-up.
To address the problems with these two types of dynamical downscaling approaches, we propose a regional reanalysis to which additional constraints are introduced from the conventional observations with the systems operated for more than about 60 years such as the surface observations at observatories and the upper-air observations with radiosondes. The longterm consistency in quality of data is essential for detection of the signals of climate change and decadal variabilities. The quality of reanalysis varies accompanied by the transitions of observing systems, if all available observations are assimilated. The assimilation of only the observations available throughout the reanalysis period enables the reanalysis to be free from the transition of observing systems, maintaining long-term homogeneity (Kobayashi et al. 2014).
Several regional reanalyses, including pilot ones, have been carried out recently covering North America (Mesinger et al. 2006), Europe (Bollmeyer et al. 2015;Jermey and Renshaw 2016;Dahlgren et al. 2016), and the Arctic (Bromwich et al. 2016). Most of them used not only conventional observations but also satellite data for assimilation, not targeting pre-satellite era (Mesinger et al. 2006;Bromwich et al. 2016;Jermey and Renshaw 2016). Bollmeyer et al. (2015) and Dahlgren et al. (2016) did not use satellite or radar observations, but assimilated data including aircraft observations, the number of which has dramatically increased for the last several decades. Besides, the tation and described anomalous local fields affected by topography, such as circulations and solar radiation, better than the coarser reanalyses.
We optimized the NHM-LETKF for long-term reanalysis by sensitivity experiments. The lateral boundary perturbations that were derived from an empirical orthogonal function analysis of JRA-55 brought stable analysis, saving computational costs. The ensemble size of at least 30 was needed, because further reduction significantly degraded the analysis. The deterministic run from non-perturbed analysis was adopted as a first guess in LETKF instead of the ensemble mean of perturbed runs, enabling reasonable simulation of spatial variability in the atmosphere and precipitation intensity.
Keywords reanalysis; regional climate; data assimilation studies on regional reanalysis did not focus on the capability of the regional reanalysis assimilating the conventional observations to overcome the problems with the two types of dynamical downscaling approaches.
We are planning to conduct a regional reanalysis assimilating only the conventional observations to produce a high-resolution dataset covering Japan and its surrounding area back to 1958 with longterm homogeneity. The purpose of this study is, as the first step for the plan, to investigate the feasibility of the regional reanalysis system assimilating only the conventional observations to estimate the past three-dimensional atmospheric states with long-term consistency in quality with particular attention paid to heavy precipitations in summer, through one-month regional reanalysis experiments before conducting a long-term regional reanalysis. An outline of the paper is as follows. In Section 2, we describe the regional reanalysis system. The experimental setup is presented in Section 3. In Section 4, the performance of the regional reanalysis system is evaluated by referring to an operational analysis and observation data. Confirming the problems with the two types of dynamical downscaling approaches, we demonstrate how much the regional reanalysis overcomes the problems. The benefits of horizontal resolution enhancement to 5 km are also investigated. The sensitivities of some important configurations for designing the system for long-term regional reanalysis are discussed in Section 5. Conclusions are summarized in Section 6.
2. The regional reanalysis system The regional reanalysis system is based on the NHM-LETKF (Kunii 2014b), which is composed of time integration with the Japan Meteorological Agency's (JMA's) nonhydrostatic model (NHM: Saito et al. 2007) and data assimilation with the local ensemble transform Kalman filter (LETKF: Hunt et al. 2007), applying the one-way double nesting approach to the Japanese 55-year reanalysis (JRA-55: Kobayashi et al. 2015). The inner and outer models cover Japan and its surrounding area with a horizontal grid spacing of 5 km and East Asia with a horizontal grid spacing of 25 km, respectively (Fig. 1). The assimilated data are limited to the conventional observations, specifically surface pressures and radiosonde observations, to ensure long-term consistency in analysis quality.

Data assimilation
For the analysis, we utilize the LETKF, an ensemble Kalman filter algorithm with high efficiency in parallel computing (Miyoshi and Yamane 2007). The first guesses are corrected by referring to the background error covariances reflecting the state-dependent uncertainties that are estimated from ensemble simulations. The observations organized in one-hour time slots are assimilated every six hours with the four-dimensional The insides of the inner and outer frames are for the NHM-LETKF with a horizontal grid spacing of 5 km and 25 km, respectively. Shades indicate the elevation in the models. The dots and crosses denote the locations of the upper and surface observations for the assimilation.
ensemble Kalman filter technique (Hunt et al. 2004). The increase in number of ensemble members improves the analysis fields, because it reduces sampling errors in the background error covariances (Kunii 2014a); however, it requires more computational resources. The localization scheme enables us to create analysis fields with moderate accuracy and reasonable ensemble size (Kunii 2014a). Miyoshi and Kunii (2012) did not show large improvement in the analysis fields with the enhancement of ensemble size from 27 in their system. From the perspective of conducting long-term reanalysis, it is crucial to save computational resources in each analysis-forecast cycle. The ensemble size is set to 30 in this system, and the localization scales are 200 km in the horizontal direction and 0.4 ln p in the vertical direction. We examine the effect of further reduction of the ensemble size in Section 5.2.
The regional reanalysis is requested to reproduce spatio-temporal variability more realistically than the global reanalysis. In the general LETKF, the first guess in the data assimilation is derived from the ensemble mean of the perturbed runs (Hunt et al. 2007). In that case, however, the analysis fields become too smoothed, resulting in the loss of realistic variability particularly in the precipitation. To avoid this problem, we use the deterministic run from the analysis as the first guess instead of the ensemble mean of the perturbed runs. The validity of this treatment is shown in Section 5.3.
When the assimilated data are limited to the conventional observations, it is apparent that the assimilated data are not uniformly distributed in space, as shown in Fig. 1. The covariance inflation with a uniform multiplicative inflation factor induces overdispersions of the ensemble members over the regions of sparse observations. The method to adaptively give inflation factors corresponding to the observation distributions can avoid the overdispersions . It is also desirable that the inflation factors reflect the possible sudden variation of observation positions in case there are rejections of observations in the quality-control process. As a method to give such inflation factors, the relaxation to prior spreads method (Whitaker and Hamill 2012) is employed in this system. The relaxation factor, a tuning parameter in this scheme, is set to 0.95 referring to Whitaker and Hamill (2012).

Time integration
The time integrations are carried out to provide the first guesses, to estimate background error covarianc-es, and to simulate the variables that are not analyzed, such as precipitation. As the model for time integration, we adopt NHM, which is a regional model governed by the fully compressive equations with various physics schemes implemented. NHM is utilized for various purposes in a range of operational short-range forecasts and regional climate research (Saito 2012).
When that time integrations are performed using a regional model, lateral boundary perturbations (LBPs) improve the analysis fields Kunii 2014b). The perturbations of global ensemble forecasts are often applied to the LBPs, but such global ensemble forecasts covering 60 years are unavailable. JRA-55, to which the one-way double nesting method is applied in this study, is a deterministic reanalysis. Performing global ensemble reforecasts to generate LBPs requires huge additional computational costs. In order to perturb lateral boundary conditions of the outer NHM-LETKF, we introduce the LBPs by adding and subtracting the modes randomly chosen in each cycle from the 15 leading empirical orthogonal function (EOF) modes of the JRA-55 climatological anomalies over the domain. The perturbations are normalized to be the amplitude of the mean sea level pressure perturbation fields of 0.7 hPa, which is equivalent to the observation error assumed in the JMA's operational Meso-scale Analysis (MA: Japan Meteorological Agency 2013). The LBPs are orthogonal to each other. Their ensemble mean is zero, their time average in each member is expected to be zero, and they do not break the physical consistency within the scope of the linear relationship. The validity of the introduction of the LBPs is confirmed in Section 5.1. For the inner model, on the other hand, the initial and lateral boundary conditions are taken simply from the analysis fields and the corresponding perturbations generated by the outer NHM-LETKF.
The tops of the inner and outer models are commonly placed at a height of 21800 m with 50 layers in the generalized hybrid vertical coordinates. The elevation of the model is given from GTOPO30. The parameters representing the surface features are determined from the land-use data from the Global Land Cover Characterization and the National Land Numerical Information. The sea surface temperatures are given from MGDSST (Kurihara et al. 2006). The same physics schemes are employed in both the inner and outer models. The adopted physics schemes are the following: the 6-class bulk microphysics scheme (Ikawa and Saito 1991); the Kain-Fritsch convective scheme (Kain and Fritsch 1990); the Mellor-Yamada-Nakanishi-Niino level 3 turbulence closure model Niino 2004, 2006); the bulk method (Beljaars and Holtslag 1991) for the estimation of sensible and latent heat fluxes from the surface using the slab land surface model, consisting of four layers for soil temperature, with the heat conduction equations and three layers for soil moisture with the force restore method; the clear-sky radiation scheme (Yabu et al. 2005); and the cloud radiation scheme (Iwasaki and Kitagawa 1998;Kitagawa 2000).

Experimental setup
As stated above, we assimilated only the conventional data, specifically the surface pressure observed at the weather stations, ships and buoys, and the wind, temperature, and relative humidity at upper levels observed with radiosondes and pilot balloons qualitycontrolled for the JMA's operational analysis (Japan Meteorological Agency 2013). An example of spatial distributions of the observations is shown in Fig. 1.
We set the experimental period to August 2014. The period includes a number of heavy precipitation events occurring in Japan, enabling us to confirm the advantages of the regional reanalysis system against the dynamical downscaling approaches and the lowerresolution reanalyses in terms of the reproducibility of the heavy precipitations in summer. The experiments covering the latest period enable us to validate the system, referring to various high-quality datasets, such as the JMA's radar-estimated precipitation calibrated with raingauge data (R/A) and the JMA's operational Meso-scale Analysis data (MA). The outer model begins at 1200 UTC 1 August 2014, and ends at 0000 UTC 1 September 2014. The initial time for the inner NHM-LETKF is lagged five days to spin up the outer NHM-LETKF.
To evaluate the regional reanalysis in terms of its advantages over the dynamical downscaling, we conducted additional experiments applying two types of dynamical downscaling approaches with the NHM without any assimilation processes. The setup of the NHM is the same as it is in the regional reanalysis experiment, except for the details described below.
In the first type of downscaling approach, the integration is serially performed over the whole experimental period using the NHM with the spectral boundary coupling (Yasunaga et al. 2005) above a height of 7 km as additional constraint by the outer data instead of the assimilation of observations.
In the second type of downscaling approach, time integrations using the NHM are carried out with the fields replaced by those interpolated from the outer data after short-term time integrations. Referring to Kayaba et al. (2016), the outer NHM is initialized to JRA-55 every six hours and driven over 12 hours. The inner NHM is initialized to the three-hour timeintegrated fields with the outer NHM. The long-term data are organized with the inner 3-to 9-hour simulated fields.
Hereafter, the experimental results of the inner models with a horizontal grid spacing of 5 km of the regional reanalysis, and the two types of downscaling approaches with serial time integrations and with reinitializations are designated by RRA, DS1 and DS2, respectively.

Comparisons of the regional reanalysis with the two types of dynamical downscaling approaches
The performance of the regional reanalysis of RRA is compared with performance of the two types of dynamical downscaling experiments of DS1 and DS2, referring to the JMA's operational analysis data (MA) and the JMA's precipitation analysis data derived from radar and raingauge observations (R/A). MA is generated with the four-dimensional variational assimilation of various kinds of observations, including the advanced observations, specifically satellites, Doppler radar radial winds, R/A, the global positioning system total precipitable water vapour data, and wind profilers (Japan Meteorological Agency 2013). The performance of MA is validated by Sawada (2008, 2009). R/A is based on radar observations and corrected with raingauge observations. R/A covers the Japan islands with a grid spacing of 1 km, and its high accuracy is verified (Nagata 2011).
a. The synoptic-scale pattern of the mean sea level pressure and upper geopotential height First, we compare the departures of RRA, DS1, and DS2 from MA in order to examine the representation of synoptic-scale fields. Figure 2 shows time series of the root mean squared differences (RMSDs) of RRA, DS1, and DS2 to MA for the mean sea level pressure (MSLP) over the inner domain shown in Fig. 1. The RMSDs of RRA are smaller compared with those of DS1 and DS2, nearly a half of that of DS1 and solidly smaller than that of DS2 for most of the period. The RMSDs of DS1 include an oscillation with a period of 5 -7 days, which cannot be found in the RMSDs of RRA and DS2. This oscillation is associated with depressions passing along a frontal zone across the domain. Figure 3 shows an example of the horizontal distributions of MSLP at 0000 UTC 26 August 2014 of RRA, DS1, DS2, JRA-55 and MA. In MA, a depression is located over the Japan Sea off the western part of the main island of Japan. RRA and DS2 describe the synoptic-scale pattern. In contrast, DS1 fails to reproduce the depression, although JRA-55, the lateral boundary conditions for DS1, successfully analyzes the pattern. Figure 4 illustrates vertical profiles of biases and RMSDs of RRA, DS1, and DS2 against MA for geopotential heights over the inner-model domain during the period. Although biases of RRA are larger in the middle-upper troposphere than those of DS1 and DS2 (Fig. 4a), RRA outperforms DS1 throughout the troposphere and outperforms DS2 below the 850 hPa surface in RMSDs (Fig. 4b). Upper observations for  the assimilations are sparse in space and time, which makes RRA have small improvements compared with DS1 and causes RRA to have larger errors than DS2 in the middle-upper troposphere. In DS1, the upper layers are constrained by forcing of spectral boundary coupling toward the boundary condition of JRA-55 above the height of 7 km, which contributes to avoiding large biases and RMSDs in the upper fields. Nevertheless, the RMSDs of RRA are still smaller than those of DS1, which emphasizes the advantages of the regional reanalysis assimilating the conventional observation over the dynamical downscaling with serial time integrations.

b. Precipitation
The simulated precipitation is evaluated against R/A data as the completely independent data. The simulated precipitation data and R/A data are averaged over the JRA-55 grids, the grid spacing of which is approximately 55 km, to tolerate position errors to some extent. The bias scores and threat scores are used for the evaluations. The bias scores indicate the ability of the system to simulate the frequencies of the phenomena exceeding a given threshold, while the threat scores evaluate the ability of the system to simulate the phenomena as per spatial and temporal patterns as well. These scores are calculated as and where FO denotes the number of samples in which both of the simulated and observed values exceed a given threshold, FX is the number of samples in which only the simulated values exceed the threshold, and XO is the number of samples in which only the observed values exceed the threshold. Figure 5a shows the bias scores of RRA, DS1, and DS2. For the bias scores, all of RRA, DS1, and DS2 have small dependence on the thresholds, which suggests that the systems have the potential to represent precipitation, including intense rain. Note that the bias scores are less than one for all of the thresholds, which implies the underestimation of precipitation in the systems. The underestimation in DS2 is more apparent relative to RRA and DS1, which implies that the spinup after reinitializations is insufficient in DS2. In terms of intense precipitation (> 30 mm 6 h -1 ), the bias scores of RRA are better than those of DS1, because DS1 has more difficulty in simulating synopticscale conditions favorable for intense precipitation, which are rarely satisfied.
The threat scores are shown in Fig. 5b. The threat scores of RRA clearly prevail over those of DS1 for all thresholds, which reflects the improvement in the representation of synoptic-scale fields by the assim- ilation of the conventional observations. RRA has better threat scores than DS2, particularly for heavier precipitations, which also attributes to the better representation of precipitation frequencies.
c. Power spectral density Figure 6 shows the power spectra of zonal wind at 500 hPa averaged over the experimental period. For RRA, both of the power spectra of the analyses and the first guesses are comparable to those of MA and DS1. The comparison of the power spectra of RRA and MA indicates that RRA represents the appropriate variability.
The power spectra of DS1 are also comparable to those of MA. Since the time integration is serially performed for long term in DS1, its spectra are assumed to represent the variability of the fields sufficiently spun up by NHM with a grid spacing of 5 km. In contrast, the power spectra of DS2 are smaller than those of RRA and DS1 although they approach the power spectra of DS1 as the integration time increases. DS2 is organized with 3-to 9-hour time integrations from the reinitialized fields derived from simulation with a grid spacing of 25 km. The small power spectra of DS2 reflect the insufficient spin-up accompanied by the reinitializations, which are consistent with the result that RRA reduces the underestimation of precipitation detected in DS2 (Fig. 5a). The comparison of the power spectra, as well as that of the precipitations, demonstrates the advantages of RRA over DS2.

d. Discussion
Forcing from the boundaries only cannot fully constrain the internal disturbances, which causes the simulated states to depart from the actual trajectory, even if the reasonable boundary conditions are provided. As the departure evolves with time integrations, the inconsistency between the inner fields and the boundary conditions can also enlarge the departure. The additional constraint by the assimilation of the conventional observations only successfully prevents the simulated synoptic-scale fields from leaving the actual states while ensuring the long-term consistency in analysis quality. The results are consistent with the study of regional reanalysis by Bollmeyer et al. (2015), in which several types of observations in addition to the conventional observations are assimilated. Our results suggest that the regional reanalysis has the potential to provide high-resolution atmospheric fields with high accuracy and long-term consistency in quality.
The regional reanalysis also has advantages over another type of dynamical downscaling including frequent reinitialization process. It is certain that this type of dynamical downscaling accurately describes the synoptic-scale fields, which is consistent with several studies (Qian et al. 2003;Lucas-Picher et al. 2013;Kayaba et al. 2016); but this approach is found to underestimate spatial variability of the fields and precipitation frequencies, owing to insufficient spin-up. In the reinitialization process, the fields are replaced with those simulated by the coarser outer model, inducing the spin-up issues. In contrast, the analysis process in the regional reanalysis just adds subtle increments to the original simulated fields. Therefore, the regional reanalysis largely improves the representation of precipitations and spatial variability of the fields. Our results reveal that the regional reanalysis overcomes the problems of not only the dynamical downscaling with serial time integrations, but also the problems of dynamical downscaling with reinitializations.

Resolution impacts
The impacts of horizontal resolution enhancement are investigated by comparing the regional reanalysis with a grid spacing of 5 km (RRA) with the global reanalysis of JRA-55 with a spacing of approximately 55 km and the regional reanalysis produced as the intermediate data in the one-way double nesting with a grid spacing of 25 km (hereafter referred to RRA25).

a. Synoptic-scale fields
Time sequences of the RMSDs of RRA, RRA25, and JRA-55 to MA for MSLP are shown in Fig. 7. There are little differences among the three reanalyses. The RMSDs of MSLP depend on the synoptic-scale fields rather than on smaller-scale fields, which are resolved only in higher-resolution reanalysis. These results imply that the regional reanalysis system including the one-way double nesting successfully works to represent the synoptic-scale fields. Figure 8 shows vertical profiles of biases and RMSDs of RRA, RRA25, and JRA-55 against MA for geopotential heights over the inner model domain during the period from 1200 UTC 8 to 1800 UTC 31 in August 2014. While the biases and RMSDs of RRA and RRA25 are small and comparable to JRA-55 near the surface, the positive biases in the lower troposphere and negative biases in the upper troposphere are found in all of RRA, RRA25, and JRA-55, which is consistent with warm bias in the lower troposphere and cold bias in the upper troposphere. The biases are enlarged in RRA25, and moderated, but still remain, in RRA. This result implies that vertical heat transport by convections is insufficient in this model, particularly in the NHM with 25-km grid spacing. The convective parameterization used in the NHM, which judges occurrence of convections and determines heat transport by convections, should be more optimized. RMSDs of RRA and RRA25 are larger than those of JRA-55. One of the reasons for this is the "double penalty", which is a well-known problem that often causes higher-resolution forecasts to have a larger amount of small-scale intensity errors. Figure 9 shows the bias scores and threat scores of RRA, RRA25, and JRA-55 against the R/A data for the precipitation averaged over the JRA-55 grids. The bias scores of RRA have small dependence on the thresholds, as previously demonstrated. By contrast, the bias scores of the JRA-55 decrease as the thresholds increase. The decline of the bias scores moderates but still remains in RRA25. From the threat scores, RRA outperforms RRA25 and JRA-55 for precipitation of stronger than 30 mm 6 h -1 .

b. Precipitations
The results suggest that the horizontal resolutions of JRA-55 (about 55 km) and RRA25 (25 km) are insufficient to represent intense precipitation. The enhancement of horizontal resolutions enables the analysis to represent the precipitation better in terms of intensity.

c. Representation of local fields in anomalously cool
days During the period of 25 -31 August 2014, a surface high over the Sea of Okhotsk was dominant, and cool air accompanying low-level clouds flowed into the northeast part of Japan. This condition causes anomalously cool days in summer, with low irradiation on the windward (eastern) side of the mountainous area running north and south over the northeast part of Japan. We examined a benefit of the high-resolution reanalysis in terms of the anomalous local fields that the coarse reanalysis did not fully resolve. Figure 10 shows the surface winds of RRA, RRA25, JRA-55 and observations at 0600 UTC averaged over the period of 25 -31 August 2014. RRA represents the cool-air inflow coming from the east coast, which is blocked by the mountainous area as  observed. By contrast, RRA25 and JRA-55 simulate the inflow crossing over the mountainous area, because the topography and accompanied local circulations are not sufficiently resolved. Figure 11 shows the downward shortwave radiation at the surface simulated with RRA, RRA25, and JRA-55 and estimated from the Moderate Resolution Imaging Spectroradiometer (MODIS) observation (Saigusa et al. 2010) averaged over the period of 25 -31 August 2014. The observed shortwave radiation (Fig. 11c) illustrates the contrast between the windward and leeward sides along the mountains running north and south over the region. The contrast is more appropriately reproduced with RRA (Fig. 11a). Particularly over the southeast part of the region, RRA simulates less than 100 W m -2 of downward shortwave radiation, which indicates that RRA reproduces the low-level clouds along the windward side. In contrast, RRA25 overestimates the shortwave radiation over the east coastal area and fails to simulate the contrast between the windward and leeward sides. Although RRA25 simulates more low-level clouds along the eastern side, the amount of clouds is not sufficient. JRA-55 overestimates the shortwave radiation and tends to estimate the cloud amount more over the lee (northwest) side. These results indicate that RRA better simulates the low-level clouds. Note that even in RRA, the contrast of shortwave radiation between the eastern and western sides is insufficient, particularly in the northern part of the region. The possible reasons for this are the difficulty of simulating the low-level clouds and the insufficient model resolution to resolve the topography. Figure 12 shows the distributions of the biases of diurnal range of surface temperature against the JMA's operational surface observations during 25 -31 August. While RRA25 overestimates the diurnal range of temperature over the eastern part of the region, RRA reduces the overestimation, particularly in the southeast part of the region, where RRA simulates the larger cloud amounts and smaller downward shortwave radiation at the surface. The biases of JRA-55 are small. Note that the 2-m temperatures in JRA-55 are obtained with assimilation of surface temperature observations applying a univariate two-dimensional optimal interpolation (Kobayashi et al. 2015), while the surface temperature observations are not assimilated in RRA and RRA25. The assimilation makes the surface temperatures in JRA-55 inconsistent with the other field variables, such as solar radiation, and obscures the features induced from the model physics. Therefore, it is difficult to discuss the differences between the temperatures of the regional reanalysis and JRA-55 in terms of the physical processes.
These results indicate that enhancement of the resolution to 5 km improves the reproduction of the anomalous local fields influenced by complex terrain such as circulations, distributions of downward shortwave radiation at the surface, and temperature variation, which are important elements from the perspective of the natural variability of local climate.

Lateral boundary perturbations
The limited area model requires lateral boundary conditions that inevitably include their uncertainties. Several studies have demonstrated that the introduction of LBPs representing the uncertainties of the lateral boundary conditions avoids underestimation of ensemble spreads especially near the lateral boundaries, contributing to the improvement of the analysis Kunii 2014b). Although the perturbations derived from global ensemble forecasts are often adopted for the LBPs, they are difficult to obtain for the LBPs of the regional reanalysis system covering the period of the last 60 years. Instead, we give the lateral boundary perturbations derived from an EOF analysis of the JRA-55 climatological anomalies as explained in Section 2.2. To justify the LBPs from the EOF analysis, we conducted an additional reanalysis experiment, hereafter denoted by RRAnolbp, in which no LBPs are given to the outer NHM-LETKF. Note that the LBPs for the inner model of RRAnolbp are given from the perturbations generated by the outer NHM-LETKF of RRAnolbp, in the same way as the LBPs for the inner model in RRA described in Section 2.2. Figure 13 shows time series of the RMSDs against MA and ensemble spreads of RRA and RRAnolbp for MSLP. The RMSDs of RRAnolbp gradually become larger than those of RRA, while the spreads of RRAnolbp gradually decrease and become much smaller than the RMSDs. The RMSE-spread ratio, an important indicator to evaluate the representativity of the uncertainties by ensemble forecasts, equals one in the ideal ensemble prediction systems. The ratio is about 0.7 for RRA, but about 0.2 for RRAnolbp in the The results suggest that without LBPs for the outer model, the underestimation of the ensemble spreads gradually becomes dominant, even in the inner domain, as the number of the analysis cycles increases. The underestimation prevents the correction of fields towards the observations due to overconfidence of the first guesses, resulting in the deterioration of analysis fields. The introduction of the LBPs from the EOF analysis can maintain the magnitude of ensemble spreads, contributing to the generation of analysis fields with steady accuracy. Since the advantages of the introduction of LBPs for the outer model are expected to be more dominant in the reanalysis covering longer periods, LBPs are indispensable in the regional reanalysis system.

Ensemble size
From the perspective of the conduction of long-term regional reanalyses, it is desirable to save computational costs for each cycle. Since the ensemble size is proportional to the computational cost of the time integrations that require the most computational resources in this system, the reduction of ensemble size significantly contributes to saving costs. In contrast, to appropriately assimilate observations, sampling errors of background error covariances estimated from ensemble perturbations should be suppressed. The reduction of ensemble size can result in deterioration of the analysis due to the larger sampling errors of the background error covariances estimated from ensemble perturbations. To explore whether the system can provide a reasonable analysis even with smaller ensemble size, we performed an additional regional reanalysis experiment with an ensemble size of 10 (RRA10). Figure 14a shows time series of the RMSDs of RRA and RRA10 against MA for MSLP. The RMSDs of RRA10 are larger than those of RRA in the latter half of the experimental period, particularly around 23 August. The degradation of RRA10 to RRA is more apparent in the 6-hour forecasts (Fig. 14b), which are used as the first-guess fields in the next data assimilation cycles. Figure 15 shows the bias and threat scores for precipitations of RRA and RRA10 against R/A. While the bias scores of RRA10 are comparable to those of RRA, the threat scores are degraded in RRA10. This result suggests that while the reduction of ensemble size does not affect the averaged features of the simulated precipitation, it degrades the simulation of precipitation in individual events corresponding to the representation of the synoptic-scale fields.

First guess
The ensemble mean fields are smooth due to suppression of uncertain modes. As the result, the ensemble mean fields represent less variability than the fields in each single run. Since the LETKF generally adopts the ensemble mean of the runs of perturbed members as the first guess, the analysis fields are also smooth. From the perspective of the application of reanalysis, the analysis fields are desirable for representing realistic variability. In order to keep reasonable variability in the analysis fields, our regional reanalysis system adopted the deterministic run from the analysis as the first guess. To justify this treatment, we conducted an additional experiment in which the ensemble mean of the perturbed runs is utilized as the first guess. Hereafter, the additional experiment is designated by RRAmean. Figure 16 shows time series of RMSDs of RRA and RRAmean against MA for MSLP over the inner domain. The RMSDs of RRAmean are nearly equal to or slightly better than those of RRA. The background error covariances in both experiments are nearly equal to the ensemble mean of perturbed runs. They cannot necessarily represent uncertainties of the deterministic run, first guess of RRA, in case the departure of the deterministic run from the ensemble mean is large owing to nonlinearity of the time integration model. Besides, the ensemble mean suppresses uncertain modes in the simulated fields. Thus, it is expected to reduce the RMSDs in RRAmean compared with those in RRA. Nevertheless, the comparison does not show clear improvement in RRAmean in terms of the RMSDs of MSLP, suggesting the validity of the use of the single deterministic run from analysis as the first guess. Figure 17 shows the power spectral densities of the two experiments, RRA and RRAmean. RRA has reasonable power spectra that are comparable to those of MA as stated in Section 4.1c. In contrast, the power spectra of the analysis fields of RRAmean are substantially smaller than those of RRA. For the 6-hour de-terministic runs from the analysis fields of RRAmean, the power spectra approach those of MA although the underestimation of power spectra still remains. The power spectra of the first guess of RRAmean, which is the ensemble mean of perturbed runs, are smaller than those of the analysis of RRAmean. These results indicate that the use of ensemble mean to obtain a first guess causes the underestimation of spatial variability and the insufficiency of spin-up. However, it may contribute to the reduction of erroneous signals coming from uncertainties in the deterministic forecasts. Figure 18 shows the bias and threat scores for 6-hour accumulated precipitations simulated from the analyses produced by the two experiments against R/A. The bias scores of RRAmean are substantially smaller for all thresholds, compared with RRA, which indicates the underestimation of precipitation frequencies in RRAmean. The precipitation systems, usually accompanied by high uncertainties, are filtered out in the ensemble mean fields, which results in the apparent underestimation of precipitation frequencies.
The threat scores in RRAmean are slightly better than RRA for the range of 5 to 30 mm 6 h -1 , which can arise not only from the slight improvement of the synoptic fields but also from the reduction of the double penalty problem in RRAmean accompanied by the underestimation of precipitation frequency.
In considering various applications of the regional reanalyses it is important that the analysis fields represent the field variables with realistic variability, particularly for the amount of precipitations. The use of the deterministic run from the analysis as the first guess is justified because the spin-up issue that accompanies the excessive field smoothing is avoided.

Conclusions
We investigated the feasibility of the regional reanalysis assimilating only the conventional observations as an alternative to dynamical downscaling for generating the past three-dimensional high-resolution atmospheric fields that were free from the transitions of observing systems back to 1958. We conducted a one-month regional reanalysis experiment over August 2014 before conducting a long-term regional reanalysis. The regional reanalysis system used the NHM-LETKF applying the one-way double nesting method to JRA-55, which covers Japan and its surrounding area with a 5-km grid spacing as the inner model and the eastern part of Asia with a 25-km grid spacing as the outer model. The assimilated data are only the conventional observations from the systems operated for about 60 years, specifically surface observations made at the observatories and upper air observations made with radiosondes, for homogeneity in the reanalysis. The results were validated with the high-quality analysis data of MA, the surface observations, and the precipitation data derived from radar and raingauge observations of R/A. This study demonstrated how much the regional reanalysis assimilating only the conventional observations can reduce the problems associated with the dynamical downscaling approaches. The regional reanalysis reproduced synoptic-scale fields considerably better and contributed to simulating temporal-spatial pattern of precipitation more accurately than the dynamical downscaling with serial long-term time integration. This type of dynamical downscaling often fails to describe the actual synoptic-scale fields, although JRA-55, which was given as the lateral boundary conditions, successfully analyzed the fields. The failure of this type of the dynamical downscaling suggested that the system was ill-conditioned from the perspective of the lateral boundary value problem. Assimilating only the conventional observations in addition to forcing from the lateral boundaries successfully improves the representation of the synoptic-scale fields throughout the troposphere, particularly in the lower layer. The regional reanalysis was also revealed to have an advantage over the dynamical downscaling with frequent reinitializations, another type of dynamical downscaling approach to reproducing higherresolution atmospheric fields. The regional reanalysis represented precipitation with realistic intensity and power spectra in appropriate magnitude, while this type of dynamical downscaling substantially underestimated precipitation and power spectra of the fields due to insufficient spin-up. These results indicated that the regional reanalysis has potential to estimate highresolution atmospheric fields, particularly precipitation, over ~ 60 years while maintaining temporal homogeneity. It was independent of transitions of observing systems, and it reduced the shortages of the dynamical downscaling approaches.
The enhancement of horizontal resolution improved estimation of the frequency of heavier rainfalls, and the it improved representation of local fields of circulations, solar radiation, and diurnal temperature changes that are largely influenced by topography, which indicated the values of the regional reanalysis for studies on extreme events. As precipitation was heavier, the simulated frequency became insufficient in JRA-55. This insufficiency could be moderated with the higher-resolution regional reanalysis. The dependence of the simulated frequency on precipitation intensity was small in the regional reanalysis with a grid spacing of 5 km. The representation of local fields in anomalous cases could also be improved in the higher-resolution reanalysis. These advantages encourage the conduction of high-resolution regional reanalysis, even when the global reanalysis is provided at the equivalent resolution to the resolution of the current operational global numerical weather forecasting model.
Through some sensitivity experiments, we examined the setup of the NHM-LETKF that was optimal for the regional reanalysis system. The lateral boundary perturbations maintained the reasonable amplitude of the perturbations for the EnKF implemented in the limited area model. Even the lateral boundary perturbations generated with the simple procedure that we introduced improve the analysis fields. An ensemble size of at least 30 was required. The reduction of the ensemble size to 10 significantly degraded the analysis. Adoption of the deterministic run from the analysis field as the first guess, rather than the ensemble mean of perturbed runs as in the general LETKF, avoided the spin-up issues that accompanied taking ensemble mean, such as excessively smoothed analysis fields and underestimation of precipitation. This modification in LETKF enabled the regional reanalysis to have realistic spatial variability in the analysis fields and to simulate precipitation with reasonable intensity. That is vital from the perspective of application of regional reanalysis.
Finally, we propose future work as follows. The experimental period in this study was limited to only one month. The conduction of regional reanalysis throughout at least one year is necessary to evaluate the re-gional reanalysis system from the scope of seasonality and some extreme events. In this reanalysis system, since the ensemble method is applied, the information about uncertainties of the reanalysis can be provided through the ensemble spreads. From this perspective, the ensemble spreads should be statistically evaluated. The system is found to underestimate precipitation. Improvement of the model is expected to simulate precipitation in more realistic quantities. Assimilating the precipitation data may also contribute to reduction of the underestimation accompanied by the reduction of the spin-up problem (Koizumi et al. 2005). It is necessary to examine the impacts of assimilation of the precipitation data observed with raingauges, a kind of conventional observation. Advancement of land surface model is also required for accurate estimation of spatially detailed distributions of surface meteorological variables. These include surface temperature and diagnosis of unavailable variables in the current system, such as accumulated snow depth, which are demanded in application of regional reanalysis. To assess the potential of a longer-term regional reanalysis, like the twentieth century reanalysis (Compo et al. 2011) or the ECMWF twentieth century reanalysis (Poli et al. 2016), an experiment assimilating only surface observations should be carried out. To explore higher-quality reanalysis with reasonable computational costs, it is valuable to examine the assimilation schemes with the use of the empirical background error covariance matrix to reduce the influence of the sampling errors, which are due to small ensemble size such as the hybrid system (e.g., Penny 2014).