The analysis of the Earth system and interactions among its spheres is increasingly important to improve the understanding of global environmental change. In this regard, Earth observation (EO) is a valuable tool for monitoring of long term changes over the land surface and its features. Although investigations commonly study environmental change by means of a single EO-based land surface variable, a joint exploitation of multivariate land surface variables covering several spheres is still rarely performed. In this regard, we present a novel methodological framework for both, the automated processing of multisource time series to generate a unified multivariate feature space, as well as the application of statistical time series analysis techniques to quantify land surface change and driving variables. In particular, we unify multivariate time series over the last two decades including vegetation greenness, surface water area, snow cover area, and climatic, as well as hydrological variables. Furthermore, the statistical time series analyses include quantification of trends, changes in seasonality, and evaluation of drivers using the recently proposed causal discovery algorithm Peter and Clark Momentary Conditional Independence (PCMCI). We demonstrate the functionality of our methodological framework using Indo-Gangetic river basins in South Asia as a case study. The time series analyses reveal increasing trends in vegetation greenness being largely dependent on water availability, decreasing trends in snow cover area being mostly negatively coupled to temperature, and trends of surface water area to be spatially heterogeneous and linked to various driving variables. Overall, the obtained results highlight the value and suitability of this methodological framework with respect to global climate change research, enabling multivariate time series preparation, derivation of detailed information on significant trends and seasonality, as well as detection of causal links with minimal user intervention. This study is the first to use multivariate time series including several EO-based variables to analyze land surface dynamics over the last two decades using the causal discovery algorithm PCMCI.