Exploring the complex relationships and drivers of ecosystem services across different geomorphological types in the Beijing-Tianjin-Hebei region, China (2000–2018)

•Trade-off scores were higher in mountainous areas than that in relatively flat areas.•The strength of association between drivers and ESs shows high spatial heterogeneity.•When diverse drivers involved in, the direction of correlations may be different.•The complex relationships mainly caused by shared and dominant drivers of ESs. Ecosystem services (ESs), as important indicators for ecosystem conditions, are always interacted complexly, however, the underlying mechanisms remain not fully understood. The Beijing-Tianjin-Hebei (BTH) region is a typical example of this complexity and is one of the most populous regions globally. In this study, we take the BTH region as the study area and try to identify the spatiotemporally complex relationships and dominant drivers of ESs, by utilizing multiple methods: Pearson correlation analysis, partial correlation analysis and the GeoDetector method. Three critically important ESs: soil erosion, water yield, and net primary production (NPP) are simulated for the period of 2000 to 2018. The results show that the Pearson correlation coefficient indicated a synergistic relationship between soil erosion and water yield and between water yield and NPP and a trade-off relationship between soil erosion and NPP, while the partial correlation analysis using the normalized difference vegetation index (NDVI) and precipitation as controlling factors showed that the relationships between soil erosion and water yield and between water yield and NPP were transformed to trade-offs. We confirmed that the drivers (especially the shared drivers and dominant drivers of each ecosystem service) played important roles in the direction of complex relationships; when different factors were involved in the relationships, the direction of positive or negative relationships may be diverse. Furthermore, the slope (q = 0.32), precipitation (q = 0.42) and NDVI (q = 0.72) were the dominant factors for soil erosion, water yield and NPP, respectively. The power of determinants for ESs showed high spatial heterogeneity in different geomorphological types due to differences in each type of inner characteristic; for example, the q value of slope for soil erosion decreased with increasing mountain relief, but the q value of NDVI increased. In addition, the trade-off degree among the three ESs was spatially heterogeneous; in mountainous areas the trade-off scores were much higher than those in relatively flat areas. Our general approaches and results provi