Competing Effects of Vegetation on Summer Temperature in North Korea

2.1. Study area

The Korean Peninsula (33–43.5 °N, 124–131.5 °E) located between China and Japan was divided into two regions along the 38th parallel after World War II. North Korea (the Democratic People’s Republic of Korea) has developed a coupled human and environmental system under different political and economic systems, compared to that of South Korea (the Republic of Korea). Under a socialist regime, the administrative division of North Korea was frequently reorganized for political purpose, and it currently consists of four special cities, which include the capital city of ‘Pyongyang’ and three special cities of ‘Nampo’, ‘Kaesong’, and ‘Rason (‘Sonbong’ in Fig. 1)’, and nine provinces (Hwanghaenam-do, Hwanghaebuk-do, Kangwon-do, Phyongannam-do, Phyonganbuk-do, Hamgyongnam-do, Hamgyongbuk-do, Ryanggang-do, and Jagang-do) as the highest-level divisions (Fig. 1). Most borders of the nine provinces (‘do’) are generally conformed to watershed boundaries which are associated with the topography of North Korea. As the major mountain ranges of the Korean Peninsula, called Baekdu-daegan, start from the highest mountainous region in the northeastern inlands of North Korea, a higher elevation is observed on the eastern side of North Korea. The lower elevated regions are the eastern coastal regions (e.g., Sinpo, Hamhung, and Wonsan) and the western regions (e.g., Pyongyang, Nampo, Hwanghaenam-do, Hwanghaebuk-do, and Phyongannam-do). The continuous mountain ranges in the eastern areas have profoundly influenced the LCLU in North Korea. For example, croplands have been broadly distributed in the western regions, where the most downstream of large rivers with alluvial plains are located.

2.2 Data

Vegetation data. We used the Normalized Difference Vegetation Index (NDVI) to detect vegetation changes in North Korea. NDVI is defined as the ratio between the near-infrared and the visible red bands of the light spectrum, and it is commonly used to indicate the amount of green biomass because it effectively detects total chlorophyll (Rouse et al. 1974). Earlier studies on monitoring deforestation and forest degradation in North Korea used the LCLU data derived from high-resolution satellite images with temporal intervals (e.g., decadal) (Kang and Choi 2014; Yu and Kim 2015; Jin et al. 2016; Choi et al. 2017; Lim et al. 2017; Dong et al. 2020). Thus, the LCLU data with high spatial resolution used in earlier studies have a relatively low temporal frequency. Rather than the higher spatial resolution LCLU data with coarse temporal resolution, the continuous time series of LCLU data with relatively coarse spatial resolution are more appropriate for matching the temporal scale of climate data in the quantitative analysis of land and atmosphere interactions (He et al. 2018). In earlier studies associated with land–atmosphere interactions, continuous NDVI data were used as an indicator of spatio-temporal vegetation changes (Lee et al. 2009, 2015; He and Lee 2016; Shull and Lee 2019). Therefore, considering the trade-off relationships between spatial and temporal resolutions, we selected NDVI time-series data derived from the Advanced Very High-Resolution Radiometer Global Inventory Monitoring and Modeling System NDVI third generation (AVHRR GIMMS NDVI3g) to identify long-term vegetation trends over the 34 years of 1982–2015, which include the 1990s when the forests in North Korea were severely destroyed. The AVHRR GIMMS NDVI3g was produced to reduce the errors and uncertainty by calibrating the AVHRR sensors, and this covers the globe with 1/12° grids (approximately 8 km) at bi-monthly intervals (15 days) from July 1981–December 2015 (Pinzon and Tucker 2014). The bi-monthly NDVI data were averaged to compose monthly time series over the 34 years of the study period (1982–2015). Using the average monthly NDVI, annual and seasonal NDVI time series were calculated for each pixel excluding the pixels with lower quality, possibly due to snow or cloud cover. The quality information of NDVI3g was represented by the percentile values in [0, 1000] for ‘good data’, [2000, 3000] for ‘spline interpolation’, and [4000, 5000] for ‘possible snow/cloud cover’ (Pinzon and Tucker 2014).

Climate data. We used meteorological observations and climate reanalysis datasets. The observation dataset of North Korea was provided by the Korea Meteorological Administration (KMA) through the Global Telecommunication System (GTS) network of the World Meteorological Organization (WMO) at 27 meteorological stations on a daily scale (see Fig. 1 and Supplementary Table S1). Using the GTS daily data, the arithmetic means of temperature and precipitation were calculated for each month of each year. Then, the average monthly temperature and precipitation were used to calculate the seasonal and annual means based on the statistical guidelines of KMA (2019). In averaging daily precipitation to a monthly mean, as suggested by KMA (2011a), the days without precipitation values were considered as non-precipitation days if at least one element among daily average, maximum, and minimum temperatures was observed; else, they would be a day with a missing value. The seasonal and annual GTS time series were used to assess the suitability of the climate reanalysis dataset during the study period (1982–2015).

For climate reanalysis, we used the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis 5th Generation-Land (ERA5-Land) dataset, which provided the climate variables related to land–atmosphere interaction with the enhanced resolution (1/10°, up to 9 km) reproduced from ERA5 reanalysis data (Muñoz-Sabater et al. 2021). The 2m temperature and total precipitation from the ERA5-Land were validated with the GTS dataset, and near-surface energy variables of surface solar radiation (SSRD), surface net solar radiation (SSR), sensible heat flux (H), latent heat flux (LE), and 2m temperature were used to investigate the land–atmosphere interactions associated with the vegetation changes in North Korea. The vertical direction of the energy flux variables in the ERA5-Land is downward in the positive values. We multiplied -1 by the original value of H and LE in each grid to represent heat energy transfers from the land to the atmosphere. Additionally, H and LE variables were used to calculate the Bowen ratio, i.e., the ratio of H to LE fluxes, to examine the dominant type of heat transfer from land to atmosphere.

2.3. Methods

Validation of climate reanalysis data. There is a limitation in acquiring reliable weather and climate data in North Korea under a tightly controlled system for generating and accessing data and information. While the observed climate data of North Korea are available through the GTS of WMO, the dataset is incomplete in statistical analysis with more than 25% missing values in 1997, 1998, 2000, 2002, 2006, and 2007 (KMA 2011a). In addition to the temporal incompletion, the GTS dataset of North Korea is relatively insufficient to cover the entire North Korea region, which is more than half of the Korean Peninsula, with only 27 stations (see Fig. 1), as compared to South Korea with 73 weather stations (KMA 2011b). To overcome the temporal and spatial limitations of the GTS dataset in North Korea, we applied climate reanalysis data with spatially and temporally continuous variables. Among the reanalysis data, we chose the ERA5-Land reanalysis data with a relatively high spatial resolution and used it to evaluate if the ERA5-Land can surrogate the GTS observations.

For the comparison of two datasets at the locations of 27 GTS stations, we used the seasonal and annual data averaged using the monthly GTS observation and the monthly ERA5-Land. Based on the seasonal and annual data of the two datasets, the means and standard deviations were calculated over the entire 27 locations during the study period of 1982–2015 and represented as a bar graph for visual comparison. To quantitatively assess the association between GTS and ERA5-Land data, we performed a parametric and non-parametric correlation analysis of Pearson’s, Spearman’s rank, and Kendall’s tau correlations. Pearson’s correlation coefficient was used as a parametric correlation analysis to represent the strength and direction of the linear relationships between two variables (Pearson 1895). As a non-parametric correlation, Spearman’s and Kendall’s rank correlation coefficients measured the strength and direction of the monotonic relationship (Spearman 1907; Kendall 1938). The three correlation coefficients were computed to evaluate the relationships between the two datasets in terms of the spatial and temporal aspects of the variances. As a spatial aspect, correlation coefficients between GTS and ERA5-Land data were calculated using the 34-year averages of seasonal and annual temperature or precipitation at the 27 locations (i.e., n=27), and the results are shown with a scatter plot of GTS and ERA5-Land data. As a temporal aspect, the correlation coefficients were estimated at each location using the 1982–2015 time-series of seasonal and annual temperature or precipitation (i.e., n=34), and then the coefficient values were represented at each location in the maps. Further, we tested the significance of the correlation coefficients for Pearson’s, Spearman’s (Best and Roberts 1975), and Kendall’s correlations (Best and Gipps 1974) at the 1% significance level.

Detection of vegetation changes. To identify the areas where vegetation has significantly changed in North Korea, we conducted a linear regression analysis using the GIMMS NDVI3g time-series from 1982 to 2015. Linear regression analysis estimated the linear relationships between the dependent variable (Y) and the independent variable (X), which was described as a linear equation, where the slope coefficient explained the rate of change in the dependent variable for every one-unit change in the independent variable (Walpole et al. 1993). At each grid point, we assigned the annual and seasonal mean of NDVI as the dependent variable and time as the independent variable, and then calculated the slope of the linear equation, which represented the NDVI trend during 1982–2015. The significance of the linear regression coefficients was tested by performing Student’s t-test at the 5% level.

Association analysis of vegetation and climate variables. As a significant change was detected during the active growing season of summer in the linear regression analysis of NDVI, we focused on investigating the association in summer (June–August, JJA) between NDVI and the near-surface energy variables, which include SSRD, SSR, H, LE, Bowen ratio, and 2m temperature. To determine the biogeophysical energy processes, a detrended correlation analysis was conducted at each grid point. First, the NDVI dataset with the finer spatial resolution was resampled to 1/10° for consistency with the ERA5-Land resolution using a bilinear interpolation method. At each grid, we then removed trends, weighted by time (year in this study) for each variable to analyze the empirical relationships between vegetation and climate variables, which could be influenced by long-term trends including global warming (Thompson et al. 2009; Lee and He 2018; He et al. 2020). After the detrended analysis, Pearson’s correlation coefficients were calculated for each grid, and the significance test was performed at the 10% level.

To quantify the magnitude of the change in temperature over vegetation, we performed a linear regression analysis using JJA NDVI as the independent variable (X) and JJA 2m temperature as the dependent variable (Y) at each grid point. The significance of the slope coefficients was tested using Student’s t-test, and that of linear regression models was tested using the F-test at the 10% level. Based on the results of linear regression, we extracted the 1982–2015 time-series of JJA NDVI and 2m temperature, which were area-averaged only over the regions with significant associations between NDVI and temperature in the western (38.0–40.0 °N, 125.1–126.5 °E) and the northern regions (40.5–41.2 °N, 125.9–127.9 °E) of North Korea. Area-averaged time-series were used to examine the vegetation and temperature associations in the western and eastern regions, respectively, and to explore the associated biogeophysical processes by performing Pearson’s correlation and linear regression analysis along with the analysis of the scatter plots.

3.1. Evaluating the ERA5-Land reanalysis with GTS observation

The 34-year means of annual and seasonal temperature and precipitation averaged over the 27 stations along with their standard deviation were visually compared between the ERA5-Land reanalysis and GTS observations (Fig. 2). Based on the GTS data, the annual average temperature was approximately 8.5 ℃, with the seasonal averages of -5.7 ℃ in winter (December through February, DJF), 8.2 ℃ in spring (March through May, MAM), 21.3 ℃ in summer (JJA), and 10.4 ℃ in autumn (September through November, SON) (Fig. 2a). Compared to the GTS data, the means from the ERA5-Land were lower by approximately 1°C for all the seasonal and annual averages. Additionally, the ranges of one standard deviation from the annual and seasonal means were similar between the two datasets. For precipitation, the annual mean from the GTS data was approximately 2.6 mm/day, with averages of 0.5 mm/day in DJF, 1.5 mm/day in MAM, 6.0 mm/day in JJA, and 2.0 mm/day in SON (Fig. 2b). The annual and seasonal means from the ERA5-Land were higher than those from the GTS by 0.4–0.8 mm/day. The ranges of one standard deviation were generally consistent between the two datasets, except for the autumn season, which had a higher deviation in the GTS than ERA5-Land. The results from the visual comparisons based on means and standard deviations supported the fact that the temperature and precipitation values from the ERA5-Land reanalysis were consistent with those from the GTS observations.

As a quantitative assessment considering the spatial aspect of variance, parametric and non-parametric correlations were estimated along with the scatter plots between the ERA5-Land and GTS data (Fig. 3 and Supplementary Fig. S1). The annual, DJF, and JJA temperatures from the ERA5-Land were linearly associated with those from the GTS over the 27 stations, although with the underestimation of the GTS (lower temperature in the ERA5-Land) (Fig. 3a, 3b, and 3c). The significantly positive associations between the ERA5-Land and the GTS temperatures were quantified as the correlation coefficients with higher than 0.95, 0.9, and 0.75 in Pearson’s, Spearman’s, and Kendall’s correlations, respectively, during all the annual and seasonal periods (all p <0.01) (Table 1a). For precipitation, statistically significant positive correlations between the ERA5-Land and the GTS were shown for all the annual and seasonal means at the 1% significance level, except for the Pearson’s r-value of MAM mean with a significant correlation at the 5% level (Table 1b). The scatter plots of precipitation showed a linear association between the ERA5-Land and the GTS with the overestimation of ERA5-Land precipitation to the GTS (more precipitation in the ERA5-Land) (Fig. 3d, 3e, and 3f). The strength of the linear relationships between the two datasets was relatively weaker during the dry season (DJF) than that in the wet season (JJA). Particularly, the highest correlation coefficients and nearly one-to-one relations between the two datasets were shown for both temperature and precipitation during the summer when vegetation and climate interactions were examined in this study.

As a quantitative assessment considering a temporal aspect of variance, the three correlation analysis methods (Pearson, Spearman, and Kendall) were conducted at each of the 27 stations using the 1982–2015 time-series of annual and seasonal temperature and precipitation. The results are shown in the maps of the Pearson’s correlation coefficients at the 27 stations (Fig. 4 and Supplementary Fig. S2) and further summarized as in the histogram (Fig. 5). The results of Spearman’s and Kendall’s correlations are shown in the supplementary materials (Supplementary Fig. S3 and S4). For temperature, there were significant positive correlations at the 1% level for the annual, DJF, and JJA means in all 27 stations (Fig. 4a, 4b, and 4c). As shown in the histogram, more than 85% of the 27 stations (more than 23 stations) had high correlation coefficients for annual and seasonal mean temperatures with Pearson’s r-values higher than 0.8, 0.8 of Spearman’s ρ, and 0.6 of Kendall’s τ (Fig. 5a). The results indicated that there was a strong association between the ERA5-Land reanalysis and GTS observation for all seasonal and annual means of temperature in nearly all of North Korea. The patterns of Pearson correlation for precipitation also showed significant positive associations at all locations (Fig. 4d, 4e and 4f), except for some stations near the northeastern coastal (i.e., Kimchaek and Sonbong) and the inland mountainous regions (i.e., Yangdok and Pungsan) for annual and JJA precipitation (Fig. 4d and 4f). At approximately 90% of the 27 stations (more than 24 stations), positive associations appeared with moderate correlation coefficients of higher than 0.5 (r), 0.5 (ρ), and 0.4 (τ) in DJF, MAM, and SON, respectively, (Fig. 5b). Compared to the drier seasons, there were relatively lower correlation coefficients (lower than 0.4) in JJA and annual means.

The means and standard deviations of temperature and precipitation between the ERA5-Land reanalysis and GTS observations were comparable (see Fig. 2), and statistically significant positive associations between the two datasets were generally quantified based on both the spatial and temporal correlations (see Table 1 and Fig. 3, 4, and 5). Based on the results of the comparison analysis, the temperature and precipitation from the ERA5-Land reanalysis were consistent with those from the GTS observations, and thus supported ERA5-Land as a suitable climate reanalysis dataset, which can substitute for incomplete observation data in North Korea.

3.2. Spatial distributions and temporal trends of vegetation

We examined the spatial patterns of annually and seasonally averaged NDVI in North Korea during the study period of 1982–2015. As shown in Fig. 6a, the annual mean values of NDVI ranged from 0.3 to 0.8 in North Korea. The higher values (0.6–0.8) were distributed mostly in the inland and northern regions, while the lower values (0.3–0.4) were in the western and near the eastern coastal regions. In the western regions, which includes the capital city of North Korea (Pyongyang), the lower NDVI values were extensively distributed even during the summer growing season, as compared to other regions (Fig. 6b). The seasonal NDVI means of DJF, MAM, and SON are shown in the supplementary materials (Supplementary Fig. S5). To detect the spatio-temporal changes in vegetation cover, the annual and seasonal NDVI trends along with their percentage changes were calculated based on linear regression analysis. Statistically significant NDVI changes primarily occurred in the western region of North Korea. In the western regions, the annual mean NDVI significantly decreased, with a trend of approximately -0.03 to -0.01 per decade (Fig. 7a), which correspond to 10–20% reduction in NDVI over the 34 years of the study period (Fig. 7b). Additionally, the seasonal means of NDVI in the western region consistently showed decreasing trends for the growing seasons throughout the spring (Supplementary Fig. S6), summer (Fig. 7c), and autumn (Supplementary Fig. S6). During JJA, the highest rate of decreasing NDVI, approximately -0.03 per decade occurred in the western region (Fig. 7c), which indicated 10–25% reduction in summer NDVI for the 34 years (Fig. 7d). The results indicated that vegetation in the western region of North Korea was intensively reduced over the thirty-four years of 1982–2015.

3.3. Associations of vegetation with near surface climates

Based on the results presented in the previous section, significantly changed vegetation was detected in the western region of North Korea with a large decrease in NDVI, and sporadically distributed, small, isolated areas of increased or decreased NDVI were detected in the northern and eastern regions. To identify how the NDVI changes can affect the near-surface climate during summer, we conducted a detrended correlation analysis of JJA NDVI with the surface energy variables of SSRD, SSR, H, LH, Bowen ratio, and 2m temperature in JJA at each grid point.

The results showed a positive correlation between NDVI and SSRD in almost all of North Korea, especially higher correlations in the northern and southern mountainous regions with an r-value higher than 0.3, which indicated a significant positive correlation at the 10% level (Fig. 8a). Contrastingly, relatively weaker positive and even negative correlations were observed in the western region (e.g., Nampo and Pyongyang, see Fig. 1 for the locations) and the eastern coastal region (e.g., Hamhung, see Fig. 1 for the location). The detrended correlation pattern of SSRD with NDVI was generally consistent with that of SSR with NDVI (Fig. 8b). The NDVI in the northwestern and southeastern regions, which showed higher r-values with both SSRD and SSR, were positively and significantly correlated with H. However, the NDVI was negatively correlated with H in the western region around Nampo (r = -0.38) (Fig. 8c). The spatial pattern of detrended correlation between LH and NDVI showed a positive association in almost all the grid points, and the statistically significant correlations were in the western and northern regions with an r-value higher than 0.3, which indicated that evapotranspiration in JJA was more active as the amount of green vegetation increased (Fig. 8d). The pattern of detrended correlation between Bowen ratio and NDVI indicated the pairs of two regions with an opposite sign of r-value (Fig. 8e), which was attributed to the results of H and LH (Fig. 8c and 8d). In the northwestern and southeastern regions, where a positive correlation for each H and LH was observed, a positive correlation of Bowen ratio with NDVI was shown. Contrastingly, a negative correlation for the Bowen ratio in the southwestern and northeastern regions, where a negative correlation for H but a positive correlation for LH, was identified. The spatial pattern of the detrended correlation analysis for temperature was generally consistent with that for the Bowen ratio. Negative correlations of temperature with NDVI were found in the western and northeastern coastal regions, especially in the western region including the significant values at the 10% level, and ranged from -0.4 to -0.3 (Fig. 8f). Contrarily, the positive associations were in the northern region, including the significant areas at the 10% level with r-values of 0.3 to 0.4, and also in inland and southeastern regions.

To estimate the change in 2m temperature over NDVI, the slope values of the linear regression between JJA NDVI and 2m temperature at each grid point were estimated (Fig. 9). The spatial pattern of linear regression was generally consistent with that of detrended correlation (see Fig. 8f), which showed negative slope values (decreasing temperature with increasing NDVI) in the western region and positive values (increasing temperature with increasing NDVI) in the northern region (Fig. 9). The significant areas identified by t-test for the slope at the 10% level are shown within the western region of 38.0–40.0 °N and 125.1–126.5 °E and the northern region of 40.5–41.2 °N and 125.9–127.9 °E, which were consistent with the significant regions tested by F-test for the linear regression models (not shown here). Specifically, as the NDVI increased by 0.1, the temperature decreased by approximately 0.4–0.6 ℃ in the western region, while it increased by approximately 0.4–0.6 ℃ in the northern region. Using the area-averaged time-series only over the grid points with significant associations between JJA NDVI and 2m temperature in the western and northern regions, Pearson’s correlation and linear regression analysis were performed, and the results are shown in Fig. 10. In the northern region, JJA temperature increased and JJA NDVI showed almost no change during 1982–2015, but for both, there was an increasing trend before the mid-1990s and no significant trend thereafter (Fig. 10a). The significant positive associations between temperature and NDVI were identified by a positive correlation coefficient (r = 0.35) and a positive slope value (0.68 ℃ per 0.1 NDVI) with significance at the 5% level (Fig. 10b). Contrastingly, in the western region, there were opposite trends between JJA NDVI and 2m temperature with a decreasing trend for NDVI and an increasing trend for 2m temperature, especially after the mid-1990s (Fig. 10c). A significant negative correlation was shown between NDVI and 2m temperature with an r-value of -0.4, and the linear regression results indicated that if NDVI increased by 0.1, then the temperature decreased by approximately 0.77 ℃ (Fig. 10d).

4.1. Competitive biogeophysical processes in North Korea

Vegetation changes including changes in vegetation density (e.g., Leaf area index, LAI), canopy structures (e.g., height and dimension), vegetation type (e.g., Plant Functional Types), and further conversion of land cover can physically affect the surface energy and moisture processes by changing the surface albedo and evapotranspiration, which leads to increasing or decreasing temperature (Foley et al. 2003; Luyssaert et al. 2014; Li et al. 2015; Zeng et al. 2017; Duveiller et al. 2018; Chen et al. 2020; Forzieri et al. 2020). While the net biogeophysical effects can be less distinct or inconsistent within the mid-latitudes due to the competing effects of albedo and evapotranspiration, the enhanced temperate forests generally induce a moderate climate by their cooling effect in summer and warming effect in winter in the areas of 35–45 °N (Li et al. 2015). Vegetation activity (e.g., LAI variation) plays a role in changing the ratio of transpiration to evapotranspiration, which closely controls the redistribution of surface energy in the temperate regions of the Northern Hemisphere (Forzieri et al. 2020). This indicates that the vegetation potentially influences temperature change through the form of latent heat flux during the summer growing season in the mid-latitude temperate zone. However, the spatial pattern of local temperature change by biogeophysical processes during the summer would be diverse due to intense competition between the two effects of albedo and evapotranspiration in the mid-latitude regions including North Korea.

The results of detrended correlation and linear regression analyses of surface energy variables with NDVI during JJA showed that the biogeophysical energy processes could occur differently between the northern regions with higher NDVI and the western regions with lower NDVI in North Korea depending on which of the two competing effects of albedo and evapotranspiration is dominant (see Fig. 8). In the northern regions, the positive associations of incoming solar energy (SSRD and SSR), sensible heat (H) and latent heat (LE) fluxes with NDVI were shown. As H increased more than LE with increasing vegetation, the increase in the ratio of H to LH (Bowen ratio) could lead to an increase in surface air temperature (2m temperature) during the summer (see Fig. 8). In an earlier study, Lee et al. (2011) identified that the intrinsic biogeophysical effect varied with latitude in the Northern Hemisphere due to the different importance among the three main factors: albedo, Bowen ratio (associated with evapotranspiration), and roughness length (almost negligible). Compared to the surrounding grasslands and bare ground, the forests to the north of 35 °N increased the temperature (warming effect), and those to the south of 35 °N decreased it (cooling effect). The associated factors were that as latitude increased, the albedo effect (amount of net shortwave absorption) was enhanced, although the effects of evapotranspiration reduced. Concerning the climatic effect of forests with altitude, the tropical mountain forests at higher altitudes induced a warming climate as the strength of the albedo effect was sustained, although that of the evapotranspiration effect decreased with increasing altitude (Zeng et al. 2021).

The northern region of North Korea, where the temperature significantly increased as the vegetation increased (see Fig. 10b), is located at latitudes higher than 35 °N (approximately 41 °N) with mountain ranges over 1000 m in elevation (see Fig. 1). Considering the geography of the northern region, the warming effect of vegetation can be dominant with enhanced vegetation resulting in increased solar radiation by decreasing albedo, increasing sensible heat flux, and thereby increasing the Bowen ratio and temperature (see Fig. 8). Thus, the warming effect could offset the cooling effect caused by the increasing latent heat flux by enhancing vegetation in the northern mountainous region. Contrastingly, the correlation results in the western region with lower NDVI, which is located at lower elevation and lower latitude, exhibited the cooling effect by increasing vegetation (see Fig. 10d). The cooling effect of vegetation was associated with the biogeophysical processes of the weak increasing or decreasing solar radiation, decreasing sensible heat flux, increasing latent heat flux, and thereby decreasing the Bowen ratio and temperature (see Fig. 8). The biogeophysical processes in the western region of North Korea were generally consistent with earlier studies in other regions across the globe, which showed that the cooling effect of vegetation was caused by decreasing the Bowen ratio along with enhanced evapotranspiration (Mahmood et al. 2014; Li et al. 2015; Zeng et al. 2017; Forzieri et al. 2020). We propose competitive biogeophysical processes in North Korea as shown in Fig. 11.

4.2. Potential impact of forest destruction in the western part of North Korea

The western region of North Korea, where a significant negative association of temperature with NDVI appeared (see Fig. 8f), experienced intensively reduced vegetation during the study period of 1982–2015. The annual trend of mean NDVI in this region was approximately -0.02 per decade, which represented a 10–20% decrease over the 34 years (see Fig. 7a). Additionally, the results of linear regression analysis for the seasonal NDVIs showed statistically significant decreasing trends during the spring, summer, and autumn (see Fig. 7b and Supplementary Fig. S6), which showed that the vegetation has continuously decreased during the growing seasons in the western region. This trend would have been associated with the forest destruction phenomenon, which is currently one of the critical environmental issues in North Korea. The forest cover change in North Korea was forced by socio-economic demands similar to other developing countries including China, Bangladesh, the Philippines, and Indonesia, which have experienced expanded deforestation driven by economic development, population density, resident income, and forest management policies of the country (Shi et al. 2017; Xu et al. 2019). In North Korea, the food and energy crisis in the mid-1990s to the early-2000s under the extreme famine, also known as the Arduous March or the March of Suffering, were the critical drivers of forest destruction. For example, local residents intensively destroyed the forest to gain food and wood fuel from the mid-1990s to the early-2000s, when food and energy shortages were experienced by internal (e.g., inefficient agricultural policies) and external forces (e.g., serious natural disasters and dissolution of the Soviet Union), which led to extreme devastation of the forest (Kong 2006; von Hippel and Hayes 2007; Stone 2012). The highest rate of devasted forests appeared in the western region (NIFS 2018; also see Fig. 7), where half of the North Korean people live (DPRK 2009). Kang and Choi (2014) stated that the forest cover fraction in the 2000s dropped to approximately half of that recorded in the 1980s over the western regions, which included the large urban areas of ‘Pyongyang’ and ‘Nampo’, and the provinces of ‘Hwanghaenam-do’, ‘Hwanghaebuk-do’ and ‘Phyongannam-do’ (see Fig. 1 for the locations). For example, the changes of forest cover fraction from the 1980s to 2000s were approximately 45–20% in ‘Pyongyang’, 35–5% in ‘Nampo’, 50–25% in ‘Hwanghaenam-do’, 60–30% in ‘Hwanghaebuk-do’, and 62–40% in ‘Phyongannam-do’. Since the 2000s, as forest development gradually reached its limit in the western region, forest cover conversion has mainly occurred in the highland areas of ‘Jagang-do’ and ‘Ryanggang-do’ in the northern mountainous regions of North Korea (Yu and Kim 2015; Dong et al. 2020). Deforestation and forest degradation were observed over the entire of North Korea, especially around suburban areas for supplying energy, building artificial parks in downtown areas, and cultivating crops (NIFS, 2018). Our findings support the fact that human-induced forest reduction has continuously appeared in the western region, which was confirmed by the significantly decreasing trend of NDVI over the recent 34 years (see Fig. 10c).

In the western region of North Korea, large-scale human-induced forest loss can induce an increase in summer temperatures through the above-mentioned biogeophysical processes. If the decrease in JJA NDVI observed in this study is maintained (-0.024/10 years, see Fig. 10c), the temperature would increase by approximately 0.18 ℃ per decade due to the negative association of temperature with vegetation. This calculation is based on the linear regression analysis of a temperature increase of 0.77 ℃ with decreasing NDVI of 0.1 during summer in the western region of North Korea (see Fig. 10d). This estimation suggests that the summer warming by forest reduction could be up to 1.5°C by the end of this century if the current rate of deforestation is not reduced. This indicates that the current forest reduction in the western region could aggravate the impacts of a warming climate on the people in the highly populated region of North Korea under future climate change. For instance, a reduction in the forests and a subsequent increase in summer temperature over the western region can enhance the frequency and strength of natural disasters (e.g., heat waves, droughts, and flash floods) and facilitate the invasion and spread of insects and parasites (KMA 2013). Consequently, approximately half of the North Korean people can be exposed to climatic disasters and climate-related health impacts. Additionally, increasing the maximum temperature can interrupt the growth of crops through the population increase and spread of exotic pests or weeds, thus reducing crop production. Considering the western region as a major breadbasket of North Korea, the enhanced warming due to forest reduction could threaten the food security of North Korea, where the agricultural policy of national self-sufficiency has been applied for food production and consumption (KMA 2013; Kim and Lee 2017). Therefore, based on our results of the forest cooling effect during the summer in the western region, we urgently propose that more effective forest management and afforestation strategies are required in North Korea, which is one of the most vulnerable regions to climate change.