## 2 Materials and Methods
### Scale-Balancing Method  
The use of different spatial scales can affect \(S_{t}\) and \(D_{t}\) as described above. Here, we need to determine an appropriate spatial scale to evaluate the risk and inequality of risk of IAPS for protected areas using PCA that uses \(S_{t}\) and \(D_{t}\) across different spatial scales. Some previous studies showed that the appropriate spatial scale was useful for SDMs such as Maxent [9; 17; 18]. Hence, we have aimed to propose a simple method to solve this specific problem. PCA is the most commonly used multivariate statistical analysis method to determine the principal components and select the appropriate variables accounting for the total variation (all scales [20; 21]). PCA could explain most of the variation between different scales and help us select the most likely spatial scale based on relationships between two or more feature sets of source data (scale data: 2.5, 5.0, 10.0, and 30.0 arcminutes) in this study [9; 17; 19]. Hence, we combined the datasets of the risk of IAPS for protected areas across all scales into a single analysis.  
First, we used \(S_{t}\) and \(D_{t}\) across different scales (2.5, 5.0, 10.0, and 30.0 arcminutes) as the PCA variables for the current concentration, low gas concentration, and high gas concentration, respectively. PCA was used to estimate the correlation matrix and percentage of variance of \(S_{t}\) and \(D_{t}\), respectively. We used the \(S_{t}\) and \(D_{t}\) values of 2.5, 5.0, 10.0, and 30.0 arcminutes as the inputs of the PCA, respectively. We conducted PCA based on the correlation among these four resolutions (i.e., 2.5, 5.0, 10.0, and 30.0 arcminutes). We extracted the loadings of the first principal component (PC1) explaining more than 60% of the variance in \(S_{t}\) and \(D_{t}\). The PC1 was significant (\(p<0.001\); Monte Carlo test) and accounted for more than a cumulative 60% of \(S_{t}\) and \(D_{t}\) across all spatial scales, respectively, for the current concentration, low gas concentration, and high gas concentration.  
Second, we extracted the \(S_{t}\) and \(D_{t}\) scores of protected areas in the first principal component (PC1), respectively. PC1 was regarded as a single index representing the risk of IAPS for protected areas across all scales. Third, we separately assessed the relationship between IAPS risk and PC1 in protected areas at all scales (\(S_{t}\) and \(D_{t}\)) using simple linear regression analysis. This analysis was used to explore the bias of IAPS risk for protectedareas between a certain scale and all scales. Fourth, we performed a simple linear regression analysis on the relationship between \(S_{t}\) and \(D_{t}\) (the results from all the scales and PC1). We aimed to explore the impact of scale effects on the relationship between \(S_{t}\) and \(D_{t}\). Finally, we chose the appropriate spatial scales used for the risk of IAPS for protected areas under climate change based on two criteria: (1) the highest representativeness, accounting for the highest amounts of \(S_{t}\) and \(D_{t}\) at all scales (and the smallest bias of risk of IAPS for protected areas between a certain scale and all scales), and (2) the smallest impact of scale effects on the relationship between \(S_{t}\) and \(D_{t}\).