**Solar Irradiance Forecasts by Mesoscale Numerical Weather Prediction Models with Di**ff**erent Horizontal Resolutions**

**Hideaki Ohtake 1,2,\*, Fumichika Uno 1,2, Takashi Oozeki 1, Syugo Hayashi 2, Junshi Ito 2,3, Akihiro Hashimoto 2, Hiromasa Yoshimura <sup>2</sup> and Yoshinori Yamada <sup>2</sup>**


Received: 7 March 2019; Accepted: 5 April 2019; Published: 9 April 2019

**Abstract:** This study examines the performance of radiation processes (shortwave and longwave radiations) using numerical weather prediction models (NWPs). NWP were calculated using four different horizontal resolutions (5, 2 and 1 km, and 500 m). Validation results on solar irradiance simulations with a horizontal resolution of 500 m indicated positive biases for direct normal irradiance dominate for the period from 09 JST (Japan Standard Time) to 15 JST. On the other hand, after 15 JST, negative biases were found. For diffused irradiance, weak negative biases were found. Validation results on upward longwave radiation found systematic negative biases of surface temperature (corresponding to approximately −2 K for summer and approximately −1 K for winter). Downward longwave radiation tended to be weak negative biases during both summer and winter. Frequency of solar irradiance suggested that the frequency of rapid variations of solar irradiance (ramp rates) from the NWP were less than those observed. Generally, GHI distributions between the four different horizontal resolutions resembled each other, although horizontal resolutions also became finer.

**Keywords:** solar irradiance forecasts; numerical weather prediction model; different horizontal resolution; forecast errors; validation; ramp rates

#### **1. Introduction**

Photovoltaic (PV) power forecasts are based on solar irradiance forecasts generated by a numerical weather prediction model (NWP) and/or machine learning predictions using data from a NWP. However, considerable errors are often included in the solar irradiance forecasts from NWPs (Schiermeier [1]). Previous studies (Ohtake et al. [2,3]) have investigated the relation between day-ahead global horizontal irradiance (GHI) forecasts and cloud types monitored in Japan Meteorological Agency's (JMA) operational model (a mesoscale model called "MSM"). These studies indicated that when certain types of clouds (i.e., stratus, altocumulus and cirrus) were monitored, GHI forecasts tended to have large one day-ahead forecast errors.

Inter-comparisons of solar irradiance and simulated cloud from NWPs with different horizontal resolutions have been performed (e.g., Mathiesen and Kleissl [4]; Perez et al. [5]; Lorenz et al. [6]), but no operational solar irradiance forecasts based on mesoscale NWPs with horizontal grid spacing finer than 2 km have been conducted. The potential usefulness of solar irradiance forecasts with a horizontal resolution finer than 1 km grid spacing is a subject of interest.

Gregory and Rikus [7] validated GHI, direct normal irradiance (DNI), and diffused irradiance (DIF) forecasts using the Australian Community Climate and Earth-System Simulator (ACCESS) operational model with a horizontal resolution of 0.11◦. DNI forecasts were found to be overestimated and DIF forecasts were found to be under-estimated, depending on a two-stream approximation asymmetry factor in a radiation process in their NWP. Sosa-Tinoco et al. [8] compared solar irradiance forecast performance using Weather Research and Forecasting model (WRF) with different horizontal resolution models (27 and 9 km mesh models) for Mexico and suggested that a higher horizontal resolution model has little benefit; however, this is partly because their target region had not drastic topography change. Charabi et al. [9] performed a case study for Oman, which was characterized by high potential of solar energy (Sunbelt) and investigated the usefulness of a higher horizontal resolution NWP considering 40, 7 and 2.8 km horizontal resolutions for assessing PV energy. Charabi et al. [9] described that the downscaled NWP model gives much more accurate estimation of the monthly as well as the annual average solar radiation. Mathiesen et al. [10] validated hourly-averaged solar irradiance forecasts for marine clouds over the southern California coast using an NWP with a 1.3 km horizontal resolution. Similarly, Lin et al. [11] compared cloud forecasts generated by multiple nested WRF simulations with horizontal resolutions ranging from 20 to 0.8 km.

Solar irradiance forecasts using higher horizontal resolutions have a higher calculation cost; thus if we choose to use higher horizontal resolution models in the future, the costs and benefits will need to be carefully evaluated.

A previous study by Ito et al. [12] reported validation results on brightness, temperature, and surface precipitation from the JMA's non-hydrostatic model (JMA-NHM), with horizontal resolutions between 5 km and 0.5 km. They did not find a significant improvement in performance between horizontal resolutions of 2 km and 0.5 km. This study validates the same products as in Ito et al. [12]'s study, but the radiation processes (shortwave and longwave radiations) are verified using surface and satellite monitoring data.

For electric power operators, forecasting short-range variability, or ramp events, of solar irradiance is important. Ramp events frequently cause rapid increases and/or decreases of PV power generation. Wellby and Engerer [13] categorize city-scale (defined as a region size of approximately 30 km) ramp events according to meso-scale and synoptic-scale meteorological phenomena (e.g., cold fronts, thunderstorms, cloud bands). If the NWP can accurately predict these meteorological phenomena, then the ramp events can also be predicted.

The remainder of this paper is organized as follows: in Section 2, surface and satellite monitoring datasets for model validations are introduced. In Section 3, model setup and model performance metrics are described. Validation results for the outputs of radiation processes are described in Section 4. In Section 5, short-term variability (ramp events) of solar irradiance is analyzed, and Section 6 summarizes the findings.

#### **2. Observation data**

#### *2.1. Global Solar Irradiance Data*

JMA observes GHI at 48 stations on the Japanese islands. GHI, DNI, DIF, upward longwave radiation (ULW), and downward longwave radiation (DLW) were observed at the JMA Aerological Observatory (Tsukuba (Tateno) station; 36◦3.4 N, 140◦7.5 E) and included in the Baseline Surface Radiation Network (BSRN) integrated database. This data has a time resolution of one second and has been verified through the JMA quality check (Gilgen et al. [14]; Ohmura et al. [15], McArthur [16]; BSRN web site [17]). The location of Tsukuba station is shown in Figure 1. Photos of instruments that measure GHI and longwave radiation at Tsukuba station are shown in Figure 2.

**Figure 1.** Topography of dx500m and the location of the Tsukuba station (yellow square). The yellow dashed rectangle indicates the Tokyo electric power company (TEPCO) area.

For GHI, DNI, and DIF observations, Kipp & Zonen CMP-21, Kipp & Zonen CHP-1 and Kipp & Zonen CMP-22 were used. Kipp & Zonen CGR-4 monitors DLW. Detailed descriptions of GHI monitoring instruments (pyranometers and pyrheliometers) in JMA stations can be found in Ohtake et al. [3]. For ULW and DLW, Kipp & Zonen CGR-4 and/or CG-4 were used.

#### *2.2. Satellite-Derived Solar Irradiance*

A geostationary satellite, Himawari-8, was launched on 7 July 2015 and is currently in operational use (Meteorological Satellite Center of JMA [18]). Detailed specifications of Himawari-8 have been reported in Bessho et al. [19]. Satellite-derived GHI datasets (including DNI and DIF), commonly referred to as "AMATERASS", is estimated based on the algorithm described in Takenaka et al. [20]. In this algorithm, cloud properties (cloud optical thickness, cloud particle size, and cloud top temperature) that influence the GHI were taken into consideration (e.g., Nakajima and Nakajima [21]; Kawamoto et al. [22]). However, aerosol optical depth (AOD) information was not sufficiently included in this algorithm. The data had a temporal resolution of 2.5 min and a horizontal resolution of 1 km for the Japanese islands and was provided by the Solar Radiation Consortium [23]. AMATERASS datasets

from Himawari-8 were also used to validate spatial distributions in GHI model simulations. Validation of the satellite-derived GHI data with ground-based data resulted in a mean bias error (MBE) within the range of 20–30 W m−<sup>2</sup> under all-sky conditions (10–15 W m−<sup>2</sup> under clear-sky conditions) and a root mean square error (RMSE) of approximately 80 W m−<sup>2</sup> (Damiani et al. [24]).

#### **3. Numerical simulations**

#### *3.1. Model*

The JMA developed the JMA-NHM (Saito et al. [25,26]; JMA web site [27]). In this study, the JMA-NHM was used to produce numerical simulations with horizontal resolutions of 5 km, 2 km, 1 km, and 500 m. Each simulation was referred to as "dx05km", "dx02km", "dx01km" and "dx500m", respectively. The experimental period in this study was the summer season from 1 July to 31 August 2015 (Exp-Summer) and the winter season from 12 January to 18 February 2016 (Exp-Winter).

Model settings for each numerical experiment are summarized in Table 1. The horizontal grid numbers (x, y) for each resolution were 220 × 180, 550 × 450, 1100 × 900, and 2200 × 1800, respectively. The number of vertical layers was 60. Model time steps were 15, 10, 5, and 3 s for each horizontal resolution, respectively.

**Table 1.** Model set-up for numerical experiments with four different horizontal resolutions (dx05km, dx02km, dx01km, and dx500m).


The three-class bulk ice (3-ICE) microphysical scheme (Lin et al. [28]; Cotton et al. [29]) was used in these numerical experiments. The Mellor-Yamada-Nananishi-Niino 2.5 (MYNN2.5) level scheme (Nakanishi and Niino [30]) was used for the planetary boundary layer process. For radiation transfer processes of the JMA-NHM, downward shortwave processes for water clouds were based on Slingo [31], while processes for ice clouds were based on Ebert and Curry [32], respectively. Longwave radiation processes were based on Hu and Stamnes [33] for water clouds and Ou and Liou [34] for ice clouds. A partial condensation scheme (Sommeria and Deardorff [35]) was introduced to treat sub-grid scale condensation because small scale clouds could not be resolved. The radiative transfer processes were calculated at 60 s intervals for every two model grids in the x and y directions to reduce calculation costs. The finer horizontal resolutions (e.g., <2 km), which are often referred as convection-allowing models, are so fine that the parameterization of cumulus is thought to be no longer required. Cumulus parametrization schemes were not used in this study.

The target area was the Tokyo electric power company (TEPCO), which is located in the center of the Japanese islands. The topography of dx500m is shown in Figure 1. Model initialization and boundary conditions were given by the JMA meso-analysis data provided every 3 h. Initialization times were 00, 06, and 18 UTC (three times a day). Twelve hour time integrations were performed during Exp-Summer in 2015 and Exp-Winter in 2016. To evaluate GHI simulations with Himawari-8 satellite data, simulation data was output every 2.5 min (see Section 2.2).

#### *3.2. Performance Metrics*

To validate forecasted GHI values, we used several evaluation parameters; MBE, mean absolute error (MAE), and RMSE. They are defined by the following equations:

$$\text{MBE} = \frac{1}{N} \sum\_{i=1}^{N} \left( \text{FCST}\_i - \text{OBS}\_i \right) \tag{1}$$

$$\text{MAE} = \frac{1}{N} \sum\_{i=1}^{N} |\text{FCST}\_i - \text{OBS}\_i| \tag{2}$$

$$\text{RMSE} = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} \left( FCST\_i - OBS\_i \right)^2} \tag{3}$$

Forecasts (FCST) and observations (OBS) are forecasted GHI and GHI surface-observed values, respectively. *N* is the sample number. Standard deviation (SD) is defined by the following equation:

$$\text{SD} = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} (x\_i - \overline{x})^2} \tag{4}$$

where *xi* and *x* represent sample data and the mean value of the sample data.

#### **4. Validation results**

#### *4.1. Validation of GHI, DNI and DIF*

First, solar radiation simulations for dx500m were compared to the BSRN observations at the Tsukuba station (The location of Tsukuba is shown in Figure 1). Figure 3a,b show time cross sections of GHI, DNI, and DIF values for the Exp-Summer of 2015.

**Figure 3.** Time cross sections (each day from July to August (horizontal axis; Month)-hours in a day (vertical axis; Hour)) for GHI (top), DNI (middle), and DIF (bottom) for the Exp-Summer of 2015 at Tsukuba station. Left panels show GHI observations, center panels show simulations by dx500m with 00 UTC (09 JST) initialization time and right panels show the difference for dx500m (dx500m–observations). Color scales are shown at the bottom of each column.

The horizontal axis represents each day of the two month period and the vertical axis represents local hours, ranging from 08 JST (Japan Standard Time, JST = UTC + 9 h) to 22 JST. Temporal variations of GHI simulations from dx500m were like that of GHI observations (Figure 3a,b). Based on the difference between the GHI observations and the simulations, negative biases were generally found during the daytime (from 12 JST to 15 JST) for the latter half of July and the first half of August. In both the first half of July and the latter half of August, positive biases were also observed.

Figure 3d–f,g–i show time cross sections of DNI and DIF simulations (like Figure 3a–c). Positive biases for DNI dominate for the period from 09 JST to 15 JST. On the other hand, after 15 JST, negative biases were found. For DIF, weak negative biases were found. Similar results to this validation were identified for each of the four different horizontal resolutions (not shown).

Generally, for Exp-Winter, GHI, DNI and DIF biases were weaker than those in summer (Figure 4). Similar positive biases of GHI and DNI values were found during the daytime in the middle of February. GHI and DNI values in the daytime were lower than those in summer because of the lower altitude of the sun. Large simulation errors were found in GHI, DNI, and DIF values in both Exp-summer and Exp-winter, although the frequency was low.

**Figure 4.** Same as Figure 3, but for the Exp-Winter of 2016.

Figure 5 shows hexbin plots for GHI, DNI, and DIF for the Exp-Summer and the Exp-Winter (e.g., Davy et al. [36]). The hexbin plots for GHI resemble a one to one line, suggesting that GHI simulations were close to GHI observations. Hexbin plots of DNI values for clear sky conditions (around 800 W m<sup>−</sup>2) and optically thick conditions (under 200 W m−2) were close to the observations. However, the middle level (from 300 to 500 W m<sup>−</sup>2) of DNI (corresponding to optically thin conditions) exhibits larger simulation errors. DIF values did not exceed 500 W m<sup>−</sup>2. DIF simulation errors tended to be large.

Results for the Exp-Winter in 2016 were also like those of the Exp-Summer in 2015 (Figure 5b). However, the frequency of GHI and DNI simulations under clear-sky conditions was higher around the one to one line when compared to the Exp-Summer. Under clear sky conditions, radiation processes in the JMA-NHM were validated. Figure 6 shows time series for GHI, DNI, and DIF simulations on 11 February 2016 at the Tsukuba station. On this date, the aerosol optical depth (AOD) was low enough (AOD < 0.2). In the case of clear sky conditions, GHI, DNI, and DIF simulations are almost consistent with surface observations, meaning that radiation processes under clear sky conditions are accurate.

**Figure 5.** Hexbin plots of GHI (top), DNI (middle), and DIF (bottom) for dx500m in (**a**) the Exp-Summer of 2015 and (**b**) the Exp-Winter of 2016 at Tsukuba station.

**Figure 6.** Comparison of (**a**) GHI, (**b**) DNI, and (**c**) DIF time series under clear sky conditions (AOD~0.2) at Tsukuba station on 11 February 2016 for dx05km (blue), dx02km (light blue), dx01km (green), and dx500m (yellow) with 00 UTC (09 JST) initialization time, respectively. Data was plotted every 2.5 min. Black circles indicate observations (OBS) and gray circles indicate extra-terrestrial solar irradiance (EXT).

Tables 2 and 3 show the three metrics (MBE, MAE, and RMSE) for simulation errors for GHI, DNI, and DIF, for Exp-Summer in 2015 and Exp-Winter in 2016. The simulations of 00 UTC initialization times are shown. From these metrics, no significant improvement was found as horizontal resolution increases. Validation results for dx05km, dx02km, and dx01km were similar to those of the dx500m.

**Table 2.** Forecast errors metrics (MBE, MAE, and RMSE) for GHI (left), DNI (center), and DIF (right) for Exp-Summer (July (upper tables) and August (lower tables)) in 2015.The forecasts with 00 UTC (09 JST) initialization time are shown. Model set-up for numerical experiments with four different horizontal resolutions (dx05km, dx02km, dx01km, and dx500m).


**Table 3.** Same as Table 2, but for the Exp-Winter (January (upper tables) and February (lower tables)) in 2016.


#### *4.2. A Case Study of Cumulus Cloud Simulation*

Short range variability of solar irradiance is one of significant issues for energy management fields. Cumulus clouds (i.e., thunderstorms) often tend to develop rapidly and cause swift decrease of GHI values (i.e., ramp events, see Section 5). Wellby and Engerer [13]) identified weather types (thunderstorm, cold front, cloud band etc.) for city-scale negative ramp events in Australia. In this subsection, a case study of cumulus clouds simulation was performed. Cloud and GHI simulations with finer resolution simulation were investigated.

Figure 7a shows images from Himawari-8 Band-3 at 14 JST on 30 July 2015. The spatial resolution for the visible band is 1 km. Developed cumulus clouds over Tokyo can be seen in the satellite image (indicated by the yellow arrows in Figure 7a). Cloud regions (noted by the dashed ellipse region "A" in Figure 7a) prevail over the northern part of the TEPCO area. Figure 7b–e show comparisons of cloud field simulations (with close-up views for the TEPCO area) with 00 UTC initialization times for each of the four different horizontal resolutions (5 h ahead clouds simulations were calculated in the radiation process). Cloud distributions from dx05km simulations tend to be mosaic. The higher the horizontal resolution in the NWPs, the finer the cloud structures are. Locations of cumulus clouds on

the northern part of the TEPCO region prevailed toward the northeastern region when compared to the satellite observations. The target cumulus clouds (shown by the yellow arrow) were not developed enough and the location of the target simulated clouds were prevailed to further west of Tokyo even in the finest resolution model, dx500m.

**Figure 7.** Comparison of the cloud horizontal distributions between the satellite observations and simulations. (**a**) Satellite visible image (Band-3) for the Tokyo electric power company (TEPCO) area at 14 JST on 30 July 2015. Cloud forecasts (5 h ahead forecasts of all cloud amounts used in radiation process) obtained from (**b**) dx05km, (**c**) 02km, (**d**) dx01km, and (**e**) dx500m with 00 UTC (09 JST) initialization time are shown, respectively. The yellow arrow in (**a**) indicates the location of the target cumulus clouds.

Figure 8 shows the validation results of GHI simulations performed at different horizontal resolutions using satellite-derived GHI observation provided by the Solar Radiation Consortium (see Section 2.2). Relatively low GHI regions (see location "A" in Figure 8a) found in the northern part of the TEPCO area were reproduced in the simulation results. As noted above, the false low GHI regions also prevailed in the northeastern region of the simulated GHI. In the location indicated by the red arrow in Figure 8a, a massive low GHI region (under approximately 200 W m−2) comprised of the target cumulus clouds was observed over Tokyo and identified from the satellite data. Simulations using different horizontal resolutions have not reproduced this low GHI region (Figure 8b–e). GHI regions in each simulation were located further west of Tokyo than what was observed in the satellite data. Generally, GHI distributions between the four different horizontal resolutions resembled each other, although finer GHI distributions were simulated as horizontal resolutions also became finer.

To validate the simulated GHI values using the satellite-derived GHI ones, Figure 9 shows frequency distributions of both the satellite-derived GHI and the simulated GHI with the four different horizontal resolutions between July and August in 2015. The target area analyzed in this study was the area around TEPCO (shown in Figures 7 and 8). Here, satellite-derived GHI data collected every 30 min were analyzed over the period of 1 July to 3 August. After 4 August, satellite-derived GHI data collected every 10 min were used in this analysis. GHI values lower than 10.0 W m−<sup>2</sup> were excluded. For the four different horizontal resolution models, 2.5 min time interval data were used.

**Figure 8.** Comparison of the GHI horizontal distributions between the satellite observations and simulations. (**a**) Satellite-derived GHI distributions for the TEPCO area at 14 JST on 30 July 2015. GHI forecasts (5 h ahead forecast) obtained from (**b**) dx05km, (**c**) dx02km, (**d**) dx01km, and (**e**) dx500m with 00 UTC (09 JST) initialization time. The red arrow in (**a**) indicates the location of the target cumulus clouds.

**Figure 9.** Frequency of (**a**) satellite-derived GHI values for TEPCO area during the two summer months of 2015. Frequency of GHI forecasts (3 h ahead forecast) obtained from (**b**) dx05km, (**c**) dx02km, (**d**) dx01km, and (**e**) dx500m with 00 UTC (09 JST) initialization time. "Num", "MEAN", "SD", and "MAX" in each panel represent data number, mean values, standard deviation and maximum values, respectively.

From this comparison, satellite-derived GHI shows a mean GHI ("MEAN") of 385.4 W m−<sup>2</sup> and a standard deviation ("SD") of 300.3 W m−<sup>2</sup> during the two summer months. From the dx05km through dx500m simulations, MEAN values were 418.9 W m−2, 424.6 W m−2, 426.2 W m−2, and 428.0 W m−2, respectively. Mean values in each simulation were generally larger than the satellite observations. SD values for the dx05km through dx500m simulations were 299.7 W m<sup>−</sup>2, 299.5 W m−2, 299.4 W m−2, and 299.4 W m−2, respectively. SD values in each simulation were significantly closer to satellite observations. Maximum values of GHI in each simulation tended to be little lower than satellite observations. The maximum value of GHI ("MAX") identified via satellite is 1100.9 W m<sup>−</sup>2. MAX values in each of the simulations were 1036.0 W m−2, 1044.2 W m−2, 1048.1 W m−2, and 1048.3 W m<sup>−</sup>2, respectively. Aerosol information (e.g., optical depth of aerosol) has not been considered in the current satellite algorithms for surface GHI value (see Section 2.2). For further verification of satellite-derived GHI values using surface GHI observations, see the previous publication by Damiani et al. [24].

#### *4.3. Longwave Radiation*

To discuss the error in terms of clouds and GHI simulations, validations regarding longwave radiation at the surface were also performed. Figure 10 shows time cross sections of ULW for both the Exp-Summer in 2015 and the Exp-Winter in 2016 at the Tsukuba station. Negative biases of ULW during the period from early morning to full daytime were significant in both summer and winter (Figure 10c,f). These results suggest that ULW (surface temperature) simulations could be under-estimated in the land surface processes of the JMA-NHM. During the mid-night (after sunset) in summer, the negative biases were significantly reduced. In case of clear sky conditions, negative biases of ULW for day time were relatively large in comparison to biases during cloudy conditions (not shown). For the Exp-Winter, positive biases were also found after 15 JST; however, daily variation of the biases was not observed. Systematic negative biases of ULW through both the Exp-Summer in 2015 and the Exp-Winter in 2016 were observed, although the validation was performed by a one-site validation in this analysis.

**Figure 10.** Time cross-sections of (**a**) ULW observations at the Tsukuba station (OBS), (**b**) dx500m forecasts (MDL) with 00 UTC (09 JST) initialization time, and (**c**) the difference (= MDL − OBS) for the Exp-Summer in 2015 (top panels) and the Exp-Winter in 2016 (bottom panels).

Figure 11 shows the validation results for DLW during the Exp-Summer in 2015 and the Exp-Winter in 2016. Generally, DLW biases tended to be weak negative biases during both summer and winter. It can be speculated that the clouds amounts simulated in the JMA-NHM were lower than they would be in a realistic atmosphere. Quantitatively, values for DLW negative biases were smaller than values for ULW biases (Figures 10 and 11).

**Figure 11.** Same as Figure 10, but for DLW.

Tables 4 and 5 show the three metrics (MBE, MAE, and RMSE) for quantifying forecast errors of both ULW and DLW for the four different horizontal resolutions for the Exp-Summer in 2015 and the Exp-Winter in 2016. These results suggested that there was little difference in the simulation errors. ULW is given by the Stefan-Boltzmann law, ~σT4. Here, is emissivity (0.95 for ground-surface), T is temperature, and <sup>σ</sup> = 5.67 <sup>×</sup> 10−<sup>8</sup> (W m−<sup>2</sup> K−4) is the Stefan-Boltzmann constant. With summer (T = 20 ◦C) and winter (T = 0 ◦C) air conditions, 1 K temperature difference would correspond to approximately 5.4 W m−<sup>2</sup> and 4.5 W m−<sup>2</sup> of ULW difference, respectively. The range of MBE values of ULW for the Exp-Summer in 2015 (−9.9 to <sup>−</sup>12.1 W m−<sup>2</sup> for July and <sup>−</sup>12.2 to <sup>−</sup>12.4 W m−<sup>2</sup> for August) corresponded to −1.8 to −2.2 K and approximately −2.3 K, respectively. For the Exp-Winter in 2016, MBE values of ULW for January (−1.1 to <sup>−</sup>3.6 W m<sup>−</sup>2) and February (−3.9 to <sup>−</sup>4.1 W m<sup>−</sup>2) corresponded to −0.2 to −0.8 K and approximately −0.9 K, respectively.


**Table 4.** Forecast error metrics (MBE, MAE, and RMSE) for ULW and DLW with the four different horizontal resolutions for the Exp-Summer (July (upper tables) and August (lower tables)) in 2015. The forecast uses a 00 UTC (09 JST) initialization time.


**Table 5.** Same as Table 4, but for the Exp-Winter season (January (upper tables) and February (lower tables)) in 2016.

Previous studies (Hara [37]; Kusabiraki [38]) reported negative temperature biases when using JMA-NHM, an operational regional model with a horizontal grid spacing of 5 km, for summer, and suggested that there were research tasks on estimations of underground thermal transport and evaporation efficiency at the surface. Weverberg et al. [39] performed an intercomparison of NWPs on radiation processes and suggested that all models in the study had too warm near-surface temperature biases and that a warmer surface temperature ought to have a larger upwelling of longwave radiation. For DLW biases, Mocrette [40] suggested that the ECMWF formulation using ice cloud optical properties from Ebert and Curry [32] and this algorithm estimated an effective particle diameter from simply diagnosed temperature, which failed to capture the full effect of clouds on the DLW.

#### **5. Ramp rate of solar irradiance**

For the electrical power and renewable energy management fields, solar irradiance (or PV power generation) ramp events are an important research subject. In this subsection, we confirmed a possibility to reproduce rapid variation of GHI, DNI, and DIF based on a higher horizontal resolution model (dx500m) and BSRN data from the Tsukuba station.

Figure 12 shows a comparison of frequency of ramp rates for GHI, DNI, and DIF values, for both observations and simulations (dx500m) for July and August in Exp-Summer of 2015. Here, ramp rate, RR, is defined by the following equation:

$$\text{RR} = \frac{\left(\mathbf{x}(t) - \mathbf{x}(t-1)\right)}{dt} \tag{5}$$

*x*(*t*) and *dt* are solar irradiance and time interval data, respectively. *dt* was set to 5 and 10 min in this paper. MEAN and SD (Equation (4)) values in legends in each figure were calculated.

Although the maximum rates during 5 min intervals (positive ramp) of GHI, DNI, and DIF observations were 936.0 W m−2, 878.0 W m−2, and 592.0 W m−2, respective positive ramp rates for simulations were 768.3 W m<sup>−</sup>2, 606.6 W m−2, 391.7 W m−2. Similar tendencies for negative ramp rates during 5 min of GHI, DNI, and DIF observations were also found (−938.0 W m−2, <sup>−</sup>844.0 W m−2, <sup>−</sup>487.0 W m−<sup>2</sup> for GHI, DNI, and DIF observation and <sup>−</sup>482.9 W m<sup>−</sup>2, <sup>−</sup>353.5 W m<sup>−</sup>2, <sup>−</sup>250.5 W m−<sup>2</sup> for those simulations, respectively). SD values of ramp rates for GHI, DNI, and DIF simulations tended to be smaller than SD values for observations. Generally, maximum and minimum ramp rates for simulations tended to be smaller than ramp rates for observations.

Ramp rates during 10 min intervals were also investigated (not shown). From SD values for simulations and observations, distributions with 5 min intervals tended to be broader than that of 10 min intervals.

Ramp rates for the January and February (Exp-Winter in 2016) were also investigated (Figure 13). In winter season, weather condition tend to be stable because Tsukuba station locates on the leeward

side of mountains in a winter monsoon. SD values for GHI, DNI, and DIF for observations (simulations) were 77.2 W m−<sup>2</sup> (49.6 W m−2), 150.5 W m−<sup>2</sup> (43.8 W m−2), 36.8 W m−<sup>2</sup> (18.2 W m−2), respectively. The SD values for GHI, DNI, and DIF for both of observations and simulations were smaller than the summer season.

**Figure 12.** Frequency of solar radiation ramp rate (per 5 min) for (**a**) observations and (**b**) dx500m forecasts in July and August of 2015 at the Tsukuba station. Results for GHI (top), DNI (middle), and DIF (bottom). Num", "MEAN", "SD", "MAX", and "MIN" in each panel mean data number, mean values, standard deviation, maximum values and minimum values, respectively.

**Figure 13.** Same as Figure 12, frequency of solar radiation ramp rate in the winter season of 2016.

NWPs with different horizontal resolutions did not reproduce the developed cumulus clouds that cause ramp events of solar irradiance simulation (see Section 4.2). Above results suggested that time variations in GHI, DNI, and DIF simulations tended to be smoother than in observations.

#### **6. Summary**

In this study, solar irradiance simulations were created using NWPs (JMA-NHM) for four different horizontal resolutions (5 km, 2 km, 1 km, and 500 m). This study examines the performance of radiation processes (shortwave and longwave radiations) using NWPs. Results were validated using surface-observed solar irradiance datasets and satellite-derived observations. The validation period was during the summer of 2015 and the winter of 2016.

GHI, DNI, and DIF simulations were validated using BSRN datasets for the Tsukuba station in Japan. Validation results showed that GHI simulations and DNI values for both clear sky conditions (around 800 W m<sup>−</sup>2) and optically thick conditions (under 200 W m−2) tend to be close to observations. DIF simulation errors tend to be large.

In addition, a case study of developed cumulus clouds on 30 July 2015 was performed. The cumulus clouds were prevailed to further west of Tokyo compared with satellite observations even in the finest resolution model, dx500m. Differences in cloud distribution relating to horizontal resolutions were small. Statistical validations using satellite-derived solar radiation found that differences of frequency of GHI values between simulations and observations was small in the TEPCO area.

Validation results on ULW found systematic negative biases of surface temperature (corresponding to approximately - 2 K for summer and approximately - 1 K for winter). The negative temperature biases could not be improved even in the finer model, dx500m. DLW biases generally tended to be weak negative biases in the summer and winter seasons.

For electric energy utility, rapid variations in solar irradiance (called "ramp events") were also investigated, using higher horizontal resolution models. From the comparison of ramp rates for solar irradiance simulations between observations and simulations, frequency of ramp rates of GHI, DNI, and DIF tended to be narrower than those of observations.

From our numerical experiments on the four different horizontal resolutions of NWPs, it was found that differences between horizontal resolutions for clouds and solar irradiance distributions were not large. In addition to enhancing the resolution of the NWP, further improvements to cloud and/or radiation schemes will be required to realize more accurate cloud distributions and/or solar irradiance forecasts. These biases could be caused by our approximate radiative calculations. As a future work, the use of a three dimensional radiative transfer is worth attempting to improve accuracy of the radiation calculations.

**Author Contributions:** Data curation, S.H., J.I., A.H., and Y.Y.; Formal analysis, H.O.; Investigation, H.O.; Methodology, H.O.; Supervision, T.O.; Writing—original draft, H.O.; Writing—review & editing, H.O., F.U., T.O., S.H., J.I., A.H., H.Y., and Y.Y.

**Funding:** This study was supported by the Core Research for Evolutional Science and Technology (CREST) project, Grant Number JPMJCR15K1, provided by the Japan Science and Technology Agency.

**Acknowledgments:** We are grateful to Osamu Ijima of the Aerological Observatory of the Japan Meteorological Agency (JMA) for providing the BSRN datasets and giving us useful comments regarding radiation measurements. We are also grateful the Solar Radiation Consortium for providing the satellite-derived GHI values. This study was carried out in collaboration with the Meteorological Research Institute (MRI) of the JMA and simulations were performed using the FUJITSU PRIMEHPC FX100 super computer at their institute. Thanks to Generic Mapping Tools (GMT) for the software that was used to draw the figures in this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. Schiermeier, Q. And Now for the Energy Forecast. Germany Works to Predict Wind and Solar Power Generation. Nature.com, 2016. Available online: http://www.nature.com/polopoly\_fs/1.20251!/menu/main/ topColumns/topLeftColumn/pdf/535212a.pdf (accessed on 6 March 2019).


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
