Quantification of the interannual variability of the nationwide electric power supply from photovoltaic systems in Japan

Information on the variation in photovoltaic (PV) power generation is essential for resource assessment. This work investigated the interannual variability of the nationwide electric power supply from PV systems in Japan. Objectives of this study were twofold: one was the quantification of the annual variability of the nationwide PV power supply. The other was identifying the causes of the variability. However, the time span of available observation data on the PV power supply is inadequate to evaluate its variability, as PV systems have been rapidly installed in recent years. We used simulation to bypass this limitation. Due to the lack of available information for modeling, a hybrid modeling approach, combining a parametric model, and estimating parameters by fitting the model to observations, was employed. Nationwide PV power supply simulations were performed using historical weather data for 30 years, from 1991–2020. The long-term simulation data enabled us to quantify the interannual variability of the nationwide PV power generation. The annual variability measured with the range from the minimum to the maximum was approximately 9% of the mean. The variability for each month was less than 30% of the monthly mean for every month except for July when it was approximately 40%. An increasing trend in the annual mean PV power supply was observed over the 30 years, with an increase of 0.16% per year of the mean over the whole period. We found that the variations in sea surface temperature (SST) in the Tropics are factors contributing to the variability of nationwide PV power supply. Specifically, the variation in SST in the tropical Indian Ocean is one of the possible driving factors of the annual variability. The framework proposed in this study can provide valuable information for assessing solar energy resources on an interannual scale.


Abstract
Information on the variation in photovoltaic (PV) power generation is essential for resource assessment. This work investigated the interannual variability of the nationwide electric power supply from PV systems in Japan. Objectives of this study were twofold: one was the quantification of the annual variability of the nationwide PV power supply. The other was identifying the causes of the variability. However, the time span of available observation data on the PV power supply is inadequate to evaluate its variability, as PV systems have been rapidly installed in recent years. We used simulation to bypass this limitation. Due to the lack of available information for modeling, a hybrid modeling approach, combining a parametric model, and estimating parameters by fitting the model to observations, was employed. Nationwide PV power supply simulations were performed using historical weather data for 30 years, from 1991-2020. The long-term simulation data enabled us to quantify the interannual variability of the nationwide PV power generation. The annual variability measured with the range from the minimum to the maximum was approximately 9% of the mean. The variability for each month was less than 30% of the monthly mean for every month except for July when it was approximately 40%. An increasing trend in the annual mean PV power supply was observed over the 30 years, with an increase of 0.16% per year of the mean over the whole period. We found that the variations in sea surface temperature (SST) in the Tropics are factors contributing to the variability of nationwide PV power supply. Specifically, the variation in SST in the tropical Indian Ocean is one of the possible driving factors of the annual variability. The framework proposed in this study can provide valuable information for assessing solar energy resources on an interannual scale.

Introduction
Photovoltaic (PV) power generation systems are expected to be key in reducing carbon dioxide emissions and mitigating global warming (IEA 2021a). The installation of PV power-generating systems is constantly increasing in many countries, and this trend is expected to continue further in the future (IEA 2021b). This situation is common in Japan. The installed PV capacity in Japan has increased rapidly since 2012 and became the primary source of renewable electricity generation in 2018. Installation of PV systems keeps increasing though the growing ratio is slowing down (IEA 2021c). Power generation by PV systems is characterized by the wide variability of its time scale from seconds to more than decades as the solar energy received at the surface fluctuates due to variations in weather and climate. Information on the characteristics of the variation in power generation by PV systems is useful for the operation of power grid systems, assessment of renewable resources, decision-making for energy policy, and investment in solar power projects. Information on interannual variability and trend is essential for project financing, as it helps evaluate the resource risk posed by solar energy (McMahan et al 2013). Surface solar energy data over several decades are required to understand the characteristics of PV power supply variation (McMahan et al 2013). When the trend is significant in the analysis period, the climatology of solar energy changes depending on the choice of the analysis period for the assessment, affecting the decision-making process for new solar power projects. Variations also influence the estimation of climatology. However, if the analysis period is chosen appropriately to filter the periodic variations, their effect on resource assessment could be minimized.
The current study aimed to understand the power supply variability from PV power-generating systems (PV power supply) on the interannual scale ranging approximately from a year to one or two decades, spanning from a regional to a national scale. This study was conducted from the view point of climatology and had two objectives. The first was the quantification of the nationwide interannual variability of solar energy in Japan. The second was an investigation of the mechanism underlying how the variability is generated. It is known that the climate system in Japan is affected by climatic variations in remote areas, like the Tropics, Midlatitude, and Arctic (e.g., JMA 2023a). Thus, there is the possibility that natural variations in solar energy in Japan are influenced by climate variations in remote regions Ohba et al (2022). investigated extreme cases of variations in the PV power supply in an area in Japan from the viewpoint of climatology, focusing on the variation on the scale of a few days. Ohmura (2009) showed the trend of solar energy, which is a longer-term variation over a few decades. Currently, there is a lack of understanding of the variability of solar energy on the interannual scale. Our study can bridge this gap and provide a new understanding of variations in solar energy on the interannual scale.
It is challenging to obtain data on the PV power supply from earlier than 10 years, as most PV powergenerating systems were installed in the past decade. The limitation in PV power generation data availability can changed in the period of the reanalysis. Changes in the amounts of assimilated information and observation instruments affect the representation of the reanalysis product (McCarty et al 2016, Randles et al 2017). Therefore, evaluating the representation of meteorological variables from reanalysis data products is necessary. In some situations, data from the reanalysis data product might be incorrect for use in resource assessment.
Data from ground observations include all atmospheric radiative effects, like clouds aerosols, water vapor, and other gaseous compositions (Stoffel 2013). Nonetheless, it has several limitations. Missing observation data are unavoidable. Discontinuity of data possibly happens because of several reasons, like accidents, replacement of observation equipment, and displacement of the observation site. Observation stations are generally located far apart, and obtaining observation data at a real PV power plant is often challenging. Additionally, strict quality control measures are essential to ensure the quality of observation data.
Additionally, two important factors should be considered while assessing long-term surface solar energy: global dimming and brightening (Müller et al 2014). The other is volcanic eruptions. Generally speaking, the dimming phase, a period of a gradual decrease in surface solar irradiance, occurred in the 1960s and continued up to the 1990s (Wild 2009). A brightening phase is followed by a gradual increase in surface solar irradiance. This sequence of events comprises different features among different regions worldwide. The trend of global solar irradiance on the surface over Japan changed from the dimming to the brightening phase between the late 1980s and early 1990s (Ohmura 2009). Haywood et al (2011) reported that decreases in atmospheric aerosols and cloud cover are the major causes of the dimming-brightening event. This trend should be evaluated separately from interannual variations when long-term data are analyzed for resource assessment. Concerning volcanic eruptions, Molineaux and Ineichen (1996), Ohmura (2009), andTanaka et al (2016) have suggested that aerosol particles emitted from volcanic eruptions affected solar irradiance for at least 3 years. As the effects of volcanic eruptions are irregular, when significant effects of the volcanic eruption occur, it should be treated as an unusual year during resource assessment analysis.
Considering the above literature review, the framework of this work was established. The dataset used in this work contains past data on the regionally aggregated PV power supply for approximately 5 years in 10 areas of Japan; however, these datasets do not include sufficient information for modeling using the physical approach. Therefore, a hybrid parameter estimation approach was adopted by fitting the parametric model to data for a limited period. First, 10 regional models were constructed; then, a nationwide model was constructed by aggregating all the regional models. Numerical simulation using a national-scale predictive model was performed to provide long-term nationwide aggregated PV power supply data. Ground-based weather observation data in Japan for 1991-2020 (30 years) are available. These data enable us to evaluate the variability due to the real natural variation and trend in the present climate. The eruption of Mt. Pinatubo in the Philippines in June 1991, one of the largest volcanic eruptions of the 20th century, provided an opportunity for evaluating the irregular impacts of volcanic eruptions on surface solar energy.
The remainder of this paper is structured as follows: section 2 explains the datasets of PV power supply and weather used for the analysis. Section 3 describes the modeling of the hourly PV power supply at the regional and national scales and methods of variability analyses. The constructed predictive models are validated in section 4. Section 5 explains the numerical simulation process of the constructed model and investigates the quantitative characteristics of the variation in PV power supply at the national scale and the relationship between the variation in the PV power supply and climate variations in the Tropics. Section 6 provides the conclusions of this work.

Data
2.1. Regionally aggregated electric power supply of photovoltaic power-generating systems Japan has 10 power transmission and distribution business operators (TDOs; figure 1(a)), whose control areas do not overlap. The TDOs have published hourly regionally aggregated PV power supply data in MW since 1st April 2016. These data do not include PV power supply for self-consumption. Data on the curtailment of PV power supply is available and it is added to the PV power supply to provide the total generated power from PV power-generating systems. The data period ranged from April 1, 2016, to December 31, 2020. The nationwide hourly PV power supply was obtained by aggregating all corresponding regional hourly PV power supplies.

Weather observation data
Hourly weather data published by the JMA were used (JMA 2023b). Five variables from the JMA dataset, namely global horizontal surface solar irradiation, air temperature, pressure, snowfall rate, and snow-cover depth, were used. For our work, data covering a long duration is required; hence, observation stations were selected based on the following criteria: (1) The period of data is 30 years and covers the period from 1991-2020; (2) The number of missing data in every month is <20 days; (3) The observation station is not located in a small island isolated from the grid system. Based on these criteria, 38 stations were selected ( figure 1(b)).
The JMA conducts quality control of each datum and provides a quality flag (JMA 2021), which is used to judge its reliability. There exist four data quality categories. Only the best and second-best data in terms of quality were used for analyses because data with higher quality are computed from sufficient observations to maintain reliability. The hourly data are computed using 10 min solar irradiance data in an hour. The best quality is defined as that including all six valid 10 min solar irradiance data, which are used to compute the hourly data. The second-best quality is defined as including five valid 10 min data and one missing or invalid data point, and the five valid 10 min data are used to compute the hourly data. The frequency of snowfall varies in different areas; therefore, only areas with snowfall for more than 50 days in a year were considered snowfall areas. Three northern areas were chosen as snowfall areas (table 1).

Climate indexes of the sea surface temperature in the Tropics
The JMA monitors the SST in several tropical areas and publishes information as indexes (JMA 2023c). Three indexes: NINO.3, NINO-WEST, and IOBW were used. NINO.3, NINO-WEST, and IOBW quantify the variations in the SST in the eastern tropical Pacific Ocean, western tropical Pacific Ocean, and tropical Indian Ocean, respectively. Definitions of indexes are explained in appendix B. Each index is monthly data, and data from 1991 to 2020 were used. Yearly data of the index is computed by averaging monthly data over a year.

Modeling for regionally aggregated photovoltaic power supply
The prediction model of the regionally aggregated PV power supply was built using the hybrid approach and has the structure of the model chain (figure 2). The model chain that converts solar radiative flux to electricity power comprises three major processes: (1) the decomposition of the global solar irradiance into direct and diffuse components; (2) the computation of irradiance incident on the surface of the panel of the array (RPOA) with or without the effect of snow cover; (3) the conversion of the RPOA to the PV power production. The prediction model includes eight parameters listed in table 2. Details of modeling are provided in appendix C.

Modeling of the nationwide aggregated photovoltaic power supply
After establishing the model for the regionally aggregated PV power supply of each area P a , where a (a = 1, K, 10) is an index of each TDO, the model for the nationwide aggregated PV power supply, P n was constructed. The national-scale model is an aggregation of all 10 regional models: If any area contains missing data for a specific date and hour, the nationwide aggregated data corresponding to that date and hour were assigned as missing data.

Estimation of parameters of the predictive model
In the predictive model, weather observation data are used as explanatory variables, and PV power supply data are used as the target variable (table 3). The data period for training the model spans approximately five years, from April 2016 to December 2020. The parameters remain constant throughout the training period, implying that the estimated parameters are representative of the entire duration, and any changes in parameters are not considered.
In the initial step of parameter estimation, the nominal capacities of the installed PV power plant systems were estimated. It was assumed that the installed capacity of PV systems would increase linearly for the training period. This assumption is based on the fact that the gradual increase of the maximum PV power supply for the training period is observed in every TDO (the results are omitted for the sake of simplicity). The maximum PV Contribution of a virtual PV power plant k (k = 1, K, s) to the aggregated power supply of an area.  power supply and its time of occurrence each year were detected. A linear regression model was used for the maximum PV power supply as a function of date and hour t:

( )
The model in equation (2) was fitted to a time series of the maximum PV power supply to estimate slope k and intercept l. Then, the normalized regionally aggregated PV power supply from the observation was computed using equation (C1).
The parameters to be inferred are those shown in table 2, excluding the azimuth angle of the panel. ψ was set to zero for all areas, which implies that every panel in an area is directed toward the south. The parameters were inferred by fitting the normalized regionally aggregated PV power supply model formulated as equation (C2) to the observations. The target function E is a form of the root-mean-square error: where P p,t and P o,t are the prediction and observation at t (t = 1, K, T), respectively. The estimation of parameters was terminated when the minimization of the target function converged. The iteration method was used to seek the minimum E. To implement the minimization procedure, the 'optim' function in the stats package of R software was used (R Core Team 2020).

Evaluation of the performance of the predictive model
The k-folding cross-validation method was used for the validation (Hastie et al 2009). Available data were used by dividing them by k. One part of data was used to estimate the prediction error, and the remaining were mixed and used to estimate the parameters during the training process. Different combinations of k parts of data were used for the estimation of parameters and validation. Then, the model was validated k times. In this study, the k value was set to five. The choice of data for k-folding cross-validation was a random choice of date and hours in the period of parameter estimation. Furthermore, the date and hours of training and validation data at each validation round were common to all areas. The aim of selecting data was to evaluate the prediction qualities of both regional and nationwide aggregated PV power supply through analyses. The prediction of the nationwide PV power supply P n, cv, t for the validation was obtained by aggregating the predictions of regional PV power supply P a, cv, t over all 10 areas a (a = 1, K, 10) at time t during the same validation round cv. ( ) The metrics for the validation are the mean-bias error MBE, the mean absolute error MAE, the root-meansquare error RMSE, and Pearson's correlation coefficient COR: where e t is the difference between predicted data x m and observed data x o at time t: e t = x m,t −x o,t . A variable with an overbar denotes the mean of the variable of T samples; σ m and σ o represent the standard deviation (SD) of T samples for prediction and observation, respectively. Additionally, normalized MBE, MAE, and RMSE were also used, which are defined as follows: The means of k values of a metric are obtained during the k-folding cross-validation process.
where M v represents a metric for the validation at the vth round of the k-folding cross-validation.
3.5. Quantification of the variability of the photovoltaic power supply Annual and monthly means were estimated to quantify the variability of the regional and nationwide PV power supply. The monthly mean was computed using the daily mean, obtained by averaging the five hourly PV electric supply values from 10:00 to 15:00 each day. If data for more than two-thirds of the daily means in a month were missing, the monthly mean of that month was considered missing data. The number of valid daily means is summarized in table A1. The annual mean was computed using the monthly means of each year. If the mean of a month in a year was unavailable, then the annual mean of the year was deemed missing data.
Linear regression fitting was used to quantify the trend of a time series. The trend was evaluated based on the slope of the linear regression model. Variability was quantified using SD, the range from the 5th to 95th percentile, which covers the middle 90 percentile of the distribution (MP90), and the range from minimum (Min) to maximum (Max), which was estimated as Max−Min.
A boxplot diagram was used to visualize the distribution of data. The box in the boxplot diagram was drawn from the first to the third quartile with a horizontal line at the median. Lower whiskers were extended to the lowest values −1.5 of the interquartile range (IQR) from the bottom of the box. Upper whiskers were extended to the highest values lesser than or equal to +1.5 of the IQR from the top of the box. Data beyond the ends of the whiskers were plotted as points, which were considered possible outliers.
The principal component analysis (PCA) is a method that is used to reduce the dimension of a multivariate dataset and obtain a new set of variables (Hastie et al 2009). The new variable is called principal components (PC). The PCA analysis is performed using the singular vector decomposition method using the statistics packages of R software (R Core Team 2020).
An eigenvector, e, and an eigenvalue, λ, of a covariance matrix [S] of the considered multivariate data x satisfy the equation.
The eigenvalue quantifies the contribution of each principal component to the total variance. The principal component with the highest eigenvalue is the first principal component (PC1), the second-highest eigenvalue is the second principal component (PC2), etc.
A periodogram, computed based on Fourier transformation, was used to estimate the spectral density. Statistical packages of R software were used for this analysis (R Core Team 2020). Periodogram I(ω) of a frequency ω is defined as the square of Fourier coefficient for ω of time-series X = {X 1 , K, X T }: Figure 3 shows the scatter diagrams between the predicted hourly regionally aggregated PV power supply and the observation of a TDO. These scatter diagrams consist of results for every validation round of the k-folding cross-validation process and do not include data used to estimate parameters. The metrics are summarized in table 4. The MBE of all areas is negative. MAE and RMSE were high, as the magnitudes of the means were high. The normalized metrics were used to compare the prediction performances of different areas. The normalized metrics indicated that the prediction quality depends on the number of observation stations in an area. There are seven stations in Kyushu and six in Tokyo and Tohoku (table 1). The models of these areas yielded a better prediction quality using normalized metrics. Okinawa contains only one station, and therefore, the associated model exhibited a lower prediction quality. However, although six observation stations exist in Hokkaido, this area is an exception with respect to the relationship between the number of observation stations and the prediction performance. Pearson's correlation coefficient was > 0.9 for each area; therefore, these regional models can reproduce realistic variations in PV power supply. Figure 4 shows the comparison of the prediction and observation of hourly nationwide aggregated PV power supply. Table 5 shows the metrics for the validation. Because MBE was negative, the national-scale model underestimated the PV power supply. As shown in figure 4(a), the national-scale model showed good prediction performance over the entire range in predicting hourly PV power supply, which yielded a high COR. Although the MBE was negative, in some cases, the predictive model yielded overestimated predictions with high values.

Results of modeling validation
The time series of the monthly mean of nationwide aggregated PV power supply from the prediction shows that the model has a good correlation with the observations (figure 4(b)). The increasing trend and variation in the monthly mean are well captured by the model. The verification of the annual mean also demonstrates good prediction quality. However, the error in the annual mean varies across different years ( figure 4(c)). The annual mean of 2018 has the highest negative error among the five years analyzed. The monthly mean of 2018 exhibits negative errors in most months, but the fluctuations in the monthly mean are synchronized. Therefore, a significant portion of the error may arise from discrepancies in estimating the nominal capacities.  5. Evaluation of nationwide variability in photovoltaic power supply 5.1. Simulation of regional and nationwide aggregated photovoltaic power supply The simulation provided data on the long-term regional and nationwide aggregated PV power supply. It aimed to provide data on the characteristics of the variation based on actual weather conditions. A numerical simulation was conducted using the national-scale PV power supply model, expressed in equation (1). The weather observation data for 30 years, from 1991-2020, were used as inputs for the model (table 3). The nominal capacities of regional models in equation (C1) was fixed at the mean computed by averaging throughout the training period, as shown in table 1; it implies that changes in installed capacity and configurations of the power systems were not considered in the simulation. The parameters in the model were the same as those determined  in section 4. These settings of the simulation are good for the quantification of the variability of the PV power supply due to the natural variation of weather and climate because the variation in simulated PV power supply was caused only by variations in weather and climate in the past 30 years. It should be noted that the simulation does not aim to reproduce the historical PV power supply. An ensemble model consisting of k models constructed using k-folding cross-validation was used for the simulation. The model outcome was the ensemble mean, obtained by simply averaging outcomes from k models. The ensemble model was used to reduce the influences of overfitting to the chosen data to estimate parameters, which were also used for the simulation. Figure 5 shows the time series of the annual mean of the simulated nationwide aggregated PV power supply showing an evident increasing trend during the period. The boxplot shows four potential outliers in data: 1991-1993 and 1998. The first 3 years were assigned as unusual years while considering the following two factors: First, the influence analysis indicated that data from 1993 exerts the highest significant influence on the robustness of the linear model, and data from 1991 also exerts a considerable influence (see appendix D). Second, these 3 years were most likely strongly influenced by the eruption of Mt. Pinatubo in June 1991 (Ohmura 2009). The linear regression model was fitted to the time series of the annual mean using the equation:

Quantification of interannual variation of the nationwide aggregated photovoltaic power supply
where Y represents the annual mean and X the year; the slope a and intercept b were estimated. The model for the entire period was named L1, and the model excluding the unusual years was named L2. Table 4 shows the results of the simulation. L1 exhibited an increase of 57.29 MW y −1 (mean increase ratio = 0.30%). L2 exhibited an increase of 30.66 MW y −1 (mean increase ratio = 0.16%). The difference in slopes between L1 and L2 was deemed significant according to the null hypothesis test with a 95% confidence level, according to Dupont and Plummer (1990).
Additionally, a variation of the linear regression model was created by adding a term that represents the effects of outliers, which was named L3: where X and Y are defined the same way as in equation (15). Z is a dummy variable distinguishing the unusual and usual years. When X 1993, Z = 1, otherwise Z = 0. The slope a, intercept b, and the effects of outliers c were estimated. The estimated slope was 30.27 MW y −1 , approximately the same as that of the linear regression model L2 (table 6). The increasing ratio to the mean was approximately 0.16%. The estimated c was −1500 MW, and the ratio to the mean was approximately −7.9%.
The detrend annual mean, which is the residual of data from the linear regression model L3 was analyzed to investigate the variability on a shorter time scale than that of the trend (figure 5). Three variability metrics are summarized in table 7. The distribution of residuals is skewed negatively because of several high negative residuals. Therefore, SD may not be suitable for the evaluation of variability. MP90 could better serve the assessment because the quantile range does not depend on the shape of a distribution. The range covers residuals in the period to facilitate the evaluation of extreme cases. Thus, the variation in the annual mean due to the natural variation in weather and climate, from which the trend was excluded, can be estimated at approximately 9% of the mean of the annual mean.
The periodicity of the annual mean was analyzed using a time series of residuals (figure 6(a)). The power spectrum density peaks at the frequency of 0.23 year −1 , corresponding to about 4.3 years cycle. The lower limit of its 95% confidence interval exceeded the red noise level. This supports that a 4-or 5-year cycle is dominant in the variation.
Analyses of annual means of simulated data for each TDO were conducted using the same analysis procedure described above to investigate the homogeneity and heterogeneity of nationwide variability in Japan. The results of these analyses are presented in appendix E.
To find structures of the variation in the annual mean, the PCA is performed. Let us consider vector X, which consists of a PV power supply from 10 TDOs each year. The vector in the n-th year is represented by X n where n = 1, K, 30. The PCA is performed for {X n }. PC1 and PC2 account for approximately 49% and 16% of the total variance, respectively, a total of approximately 65% of the total variance. Eigenvectors of PC1 correspond to the situation when the PV power supply of all TDOs except Okinawa and Hokkaido varies simultaneously with the same sign ( figure 7(a)). PC2 corresponds to the situation when there are two groups of TDOs in Japan and the PV power supply of two groups fluctuates with opposite signs: one group consists of three TDOs: Kyushu, Chugoku, and Shikoku, while the other group consists of Okinawa, Hokuriku, Tohoku, and Hokkaido ( figure 7(b)). However, Kansai, Chubu, and Tokyo have lower eigenvector values.
Time-series of PC1 and PC2 is shown in figure 8. The time-series of PC1 has similar characteristics to the time-series of the residual of the annual mean of the nationwide aggregated PV power supply ( figure 8(d)). Pearson's correlation coefficient between those two time-series is +0.97. The analysis of the power spectrum for PC1 shows that its dominant period is approximately 4-5 years, which is similar to the results of the nationwide aggregated PV power supply (figure 6(c)).
The time-series of PC2 shows a slower variation than those of PC1 ( figure 8(e)). The correlation coefficient between PC2 and the nationwide aggregated PV power supply is −0.01. PC2 corresponds to the situation when two groups whose variations have opposite signs. Then, fluctuations of the two groups are canceled out. Thus, PC2 does not strongly affect the nationwide aggregated variation. PC2 does not have a confident result of the periodicity analysis (figure 6(d)). However, the result indicates that one of the dominant periods of PC2 can be estimated to be a few decades. The length of data is not long enough to obtain confident results of the periodicity of a longer period than a decade. It can be concluded that the variation in the annual mean of the nationwide aggregated PV power supply is mainly controlled by the variation mode of PC1.

Seasonality of the variability of the nationwide aggregated photovoltaic power supply
The variability in the monthly mean of each month was quantified using the range of MP90 by following the same method as the variability estimation in the annual mean. Figure 9(a) shows the median, the mean of monthly mean, and the variability quantified with MP90 in each month. Figure 9(b) also shows variability, where the median for the monthly mean is placed as a horizontal line at zero. The variability for July was observed to be the highest. April and May also showed high variability. The variability during winter tended to be low. The months in which the monthly mean is higher tended to exhibit higher variability. Thus, the range of MP90 of a month was divided by the monthly mean to reduce the relationship between the magnitude of the monthly mean and variability ( figure 9(c)). The variability of July was still the highest after the modification. The ratio of MP90 in July to the mean was approximately 40%. However, the ratio of other months was lower than 30%. The result for each TOD is shown in appendix E in the same manner as in figure 9(a). The PCA was performed for the monthly mean of each month with the same procedure of the analysis of the annual mean. PC1 in every month has a similar structure of the eigenvector as that of the annual mean shown in figure 9(a), in which most TDOs fluctuate simultaneously with the same sign ( figure F1). Additionally, PC1 at each month has a strong correlation with the monthly mean (table 8), indicating that the variability of each month is controlled by the variation mode of PC1. The contribution of PC1 to the total variance in July is the Horizontal bars represent the confidence interval of 95% level at each frequency. The curve represents the periodogram of the red noise. When the curve of the red noise is lower than the lowest limit of the confidence interval at a frequency, the periodogram at the frequency is statistically confident. Each time series is standardized using its mean and standard deviation in the spectrum analysis. highest in the year, which possibly causes the variability of the nationwide aggregated PV power supply in July to be the highest in the year.

Connection between the variation in the photovoltaic power supply in Japan and variations in sea surface temperature in the Tropics
To investigate the mechanism of occurrence of the variability of the nationwide aggregated PV power supply in Japan, the relationship between the variation in the PV power supply and climate indexes was analyzed. The correlation coefficient between PC1 of the annual means and each yearly climate index is shown in table 9. Timeseries of each yearly climate index is shown in figure 8. The IOBW has a weak negative correlation with PC1 at the 95% confidence level. The IOBW also negatively correlates with the nationwide aggregated PV power supply (table 9). Therefore, this result suggests that the variation in SST in the tropical Indian Ocean is one of the drives of the nationwide variation in solar energy in Japan. When the SST in the Indian Ocean is higher (lower) than usual, the solar energy in most of Japan tends to be lower (higher) than usual. The index which has the strongest correlation with PC2 is NINO3 though the correlation is weak at a confidence level of 90%, indicating that when the SST in the tropical Pacific Ocean is higher (lower), the PV power supply in western Japan tends to be lower (higher) than usual, and that in eastern Japan tends to be higher (lower).
The correlation coefficient between PC1 of the monthly means and each of the monthly climate indexes in each month shows different features (table 8) and changes month by month. The IOBW tends to have a negative relationship with PC1 in the whole year except in September and October. The highest negative correlation coefficient was observed in July. The correlation coefficient between PC1 and NINO.WEST in each month tends to have almost the opposite sign to that between PC1 and IOBW. The correlation coefficient between PC1 and NINO.3 tends to become stronger in winter; however, the tendency of the seasonal change of the correlation coefficient is unclear in summer, in a usual year, the rainy season in Japan ends in July. The onset of the rainy season in a usual year in Japan starts in the southern part of Japan in May, and it migrates northward and eastward gradually. The withdrawal of the rainy season also starts in southern Japan and migrates northward and eastward gradually; then, the rainy season ends at the end of July over whole areas in Japan in a usual year. The occurrence of the highest variability of the nationwide aggregated PV power supply matches the withdrawal of the rainy season. Takaya et al (2020) reported that the variation in the SST in the Indian Ocean affects the activity   of the frontal zone in the rainy season, where the weather condition tends to be cloudy and rainy. The highest variability of solar energy in July seems to be related to the activity of the front.

Conclusions
Major conclusions were arrived at based on the quantification of the variability of the nationwide aggregated PV power supply: (1) The time series of the annual mean showed an increasing trend over the 30 years. The increasing ratio of the mean of the annual mean was estimated at 0.16% per year.
(2) The variability of the annual mean was quantified using three metrics. The range from the minimum to the maximum of the annual mean did not exceed 10% of the mean. Table 8. Correlation coefficient between the first principal component (PC1) and monthly mean of the nationwide aggregated PV power supply or each climate index in each month. The marks * and ** indicate that the null hypothesis is rejected at the 10 and 5% significance levels, respectively.  Table 9. Pearson's correlation coefficient between the annual mean or the first (PC1) and second (PC2) principal components of the nationwide aggregated PV power supply and each of the climate indexes. The marks * and ** indicate that the null hypothesis is rejected at the 10 and 5% significance levels, respectively. (3) Variability also showed seasonal changes. Variability was quantified using the range from the 5th to 95th percentile (MP90), and it was <30% of each monthly mean, except for July when it was about 40%.
(4) The dominant period of variation in the annual mean was estimated to be 4 or 5 years.
(5) Irregular effects related to the volcanic eruption of Mt. Pinatubo were estimated. A reduction of 7.9% to the mean of the annual mean was observed correspondingly.
Variations in the SST in the Tropics are one of the possible driving factors of the variability of the annual mean of the nationwide aggregated PV power supply in Japan. The variation in the SST in the Indian Ocean is one of influential factors of the variability in the nationwide solar energy over Japan. The higher (lower) SST in the Indian Ocean causes the lower (higher) PV of the nationwide power supply in Japan. There is a seasonal change in the effects of the variations in the SST on monthly solar energy in Japan. The SST in the Indian Ocean weakly correlates with the monthly solar energy in January and July, but a clear correlation is not observed in other months.
Though the climate variations originating from variations in the SST in the Tropics are investigated in this work, climate variations in other areas, such as the Midlatitude and Arctic regions, also affect the climate in Japan. Therefore, to understand the variability of solar energy more deeply and completely, more complicated climate systems should be investigated in future works. Data for 30 years is insufficient to investigate climatescale phenomena whose time scale exceeds a few decades. Therefore, data for longer periods are required to obtain robust results. Using reanalysis data or data derived from climate simulations is an option to get longterm data. However, validation of such data is necessary to judge their suitability for the investigation.
First, the regionally aggregated PV power supply in area a, which is represented as P a , is normalized by the nominal capacity of the PV power plants in the area P a,m to obtain the normalized regionally aggregated PV power supply P′ a , as follows: FSC. In the PV snow coverage model, the snow sliding on the panel surface and the reduction ratio of the snow per hour due to sliding were considered. The occurrence of sliding at a site k is governed by the following inequality: where m is a parameter used to characterize the relationship, and its unit is W/m 2 /°C, and = RPOA RPOA k a q j Rdir Rdiff r , , , , .
) When the equation is satisfied, snow sliding occurs. If the occurrence of sliding is recognized based on equation (7), the reduction ratio of snow cover on the POA is estimated from the quantity of sliding of snow (QSS) in an hour, which is provided by the following equation: where v is the sliding coefficient. FSC is computed using the following procedure: To simplify expressions, let FSC(t) express the FSC at time t, and QSS (t) the QSS at time t.
(3) If the occurrence of snow sliding is recognized, FSC(t+1) is computed using the following relationship: 1 1 in the case of 1 . In addition, 1 0 in the case of 1 (4) Snow cover becomes zero after 24 h of snowfall, even when the snow remains on the PV panel surface, implying that the valid period of snowfall is 24 h. If a snowfall occurs within 24 h after the first snowfall, the valid period of the effect of snow cover restarts from the time of the end of the second snowfall. The valid period is introduced to include other effects, such as the removal of snow cover due to melting and wind.

Appendix D. Influence analysis and outlier test
Influential analysis and outlier tests for a univariate regression model were performed according to a previous study (Weisberg 2005). The influential analysis aims to study the robustness of a model against perturbations in predictors. Cook's distance was used to evaluate the influence of a data point on the model parameter. Scenarios in which Cook's distance was high suggest a substantial influence on both the model and fitted values. The 4/(n −k−1) is an empirical threshold of Cook's distance to consider the possibility of the influence, where n is the sample size and k is the number of explanatory variables (Fox 2015). When the Cook's distance exceeds 1, data are most likely influential on the robustness of the model. The outlier test aimed to detect data points that do not follow the same model as the rest of data. The Bonferroni outlier test was performed for each data point. The test statistic in the test was the studentized residual. The high absolute value of the studentized residual indicates the possibility of the outlier. The Cook's distance for every annual mean is shown in figure D1(a). The annual mean of 1993 had the highest Cook's distance, and that of 1991 the second highest. These two data were significantly higher than the rest and exceeded the threshold of 0.1429 (≈ 4/28). The outlier test concluded that there were no outliers with a 5% significance level ( figure D1(b)). The p value of 1993 was about 0.9. However, the studentized residual of 1993 was the highest, and that of 1991 was the third highest. Figure E1 shows the time series of the annual mean of each TDO in the same manner as in figure 5(a). Moreover, table E1 summarizes the results of regression analyses for 10 TDOs. Notable points of the results are as follows: an increasing trend is detected in four TDOs: Kansai, Chubu, Tokyo, and Tohoku. Effects of unusual years range approximately from 1 to 10% to the mean. Effects of unusual years of Kyushu and Tokyo are higher than 9% in the ratio to the mean. TDOs in western Japan, including Kyushu, Chugoku, and Shikoku, tend to have high effects of unusual years. Figure E2 shows the time series of the residual of the annual mean in the same manner as in figure 5(c). Table E2 shows the quantification of the variability of the annual mean for each TDO. The variability measured with MP90 ranges from 7 to 13% in the percentile to the mean. The variabilities of Okinawa and Tokyo are the first and second highest among TODs, respectively. The range of each TDO exceeds 10% of the mean and is larger than that of the nationwide annual mean, implying that the aggregation of TDOs causes a decrease in the variability of solar energy. Figure E3 shows the variability in the monthly mean for each TDO. There are differences and similarities in the seasonality of the variability among TDOs. All TDOs except Okinawa and Hokkaido have two peaks: one in April or May, and the other in August. The variabilities of those TDOs are the highest in July except in Kyushu, where the highest variability is observed in June. In August, Okinawa has only one peak, and the variability reaches its maximum in May. Hokkaido has a peak in April, and its maximum variability also occurs in April. Seven TDOs had the highest variability in July.