AbstractHighly sensitive laser strainmeters at Moxa Geodynamic Observatory (MGO) measure motions of the upper Earth's crust. Since the mountain overburden of the laser strainmeters installed in the gallery of the observatory is relatively low, the recorded time series are strongly influenced by local meteorological phenomena. To estimate the nonlinear effect of the meteorological variables on strain measurements in a non‐stationary environment, advanced methods capable of learning the nonlinearity and discovering causal relationships in the non‐stationary multivariate tectonic‐climate time series are needed. Methods for causal inference generally perform well in identifying linear causal relationships but often struggle to retrieve complex nonlinear causal structures prevalent in real‐world systems. This work presents a novel model invariance‐based causal discovery (CDMI) method that utilizes deep networks to model nonlinearity in a multivariate time series system. We propose to use the theoretically well‐established Knockoffs framework to generate in‐distribution, uncorrelated copies of the original data as interventional variables and test the model invariance for causal discovery. To deal with the non‐stationary behavior of the tectonic‐climate time series recorded at the MGO, we propose a regime identification approach that we apply before causal analysis to generate segments of time series that possess locally consistent statistical properties. First, we evaluate our method on synthetically generated time series by comparing it to other causal analysis methods. We then investigate the hypothesized effect of meteorological variables on strain measurements. Our approach outperforms other causality methods and provides meaningful insights into tectonic‐climate causal interactions.