The hilly and gully region of the Loess Plateau represents one of China’s most ecologically vulnerable landscapes, characterized by severe soil erosion, intensive land use, and pronounced disturbances to the structure and functionality of ecosystem services. Taking Zichang City as a case study, this research integrates grid-scale analysis with the InVEST-PLUS model and bivariate spatial autocorrelation techniques to examine the spatiotemporal dynamics and inter-relations of four critical ecosystem services—carbon storage, water yield, biodiversity, and soil retention—under varying land use intensity scenarios from 1990 to 2035. The findings indicate that (1) between 1990 and 2020, land use intensity in Zichang City steadily declined, exhibiting a spatial distribution pattern typified by central-area clustering and gradual peripheral transitions. (2) Across three development scenarios, the spatial distribution of the four ecosystem services aligned with the patterns observed in 2020, with central areas showing pronounced fluctuations, whereas peripheral regions experienced relatively minor changes. Specifically, from 1990 to 2020, the proportion of low-carbon storage areas increased by 2.89%, and high water yield areas expanded by 9.45%, while the shares of low habitat quality and low soil retention areas decreased by 5.59% and 6.25%, respectively. (3) A significant spatial autocorrelation was observed between land use intensity and the four ecosystem services, with widespread cold and hot spots reflecting dynamic spatial clustering patterns. These results offer valuable insights for optimizing land use strategies, improving ecosystem service performance, and advancing ecological conservation and sustainable development initiatives.