An exploration of the complex relationships among ecosystem services (ESs) is critical for achieving comprehensive, coordinated and sustainable human development. However, studying the optimal relationships of multiple ESs and what drives them remains an enormous challenge, especially in lake basins. Therefore, we chose the Nansi Lake Basin (NLB), the largest freshwater lake basin in North China, for ES evaluations in 2010, 2015, and 2020. The Pareto optimal method was used to screen out the Pareto optimal points (POPs) representing the optimal relationships among the three ESs. The spatial distribution characteristics and drivers of the POPs were analysed using kernel density estimation (KDE) and an optimal parameter -based geographic detector (OPGD). The results showed that the number of POPs was low, accounting for only 0.30% - 0.88% of the total samples. The distribution of POPs exhibited spatial clustering. Areas with high density values (top 30%) for ES relationships that included water yield (WY) gradually moved eastward and southward from 2010 to 2020. These ESs, which had a significant relationship in the two-dimensional (2D) analysis, significantly affected the shape of the threedimensional (3D) surface formed by the POPs. Climatic factors (temperature and precipitation), topographic factors (elevation and slope) and landscape fragmentation all greatly influenced the POP distributions of the different ES relationships during the research periods. We propose a method to explore the optimal relationships among 3D ESs and provide valuable insight for future research on multiple ES relationships in lake basins.