Targeted Calibration to Adjust Stability Biases in Complex Dynamical System Models

TL;DR

The paper tackles the problem of stabilizing or destabilizing equilibrium states in non-differentiable, high-dimensional climate models where conventional gradient-based tuning is infeasible. It introduces a targeted calibration method that builds a local VAR(1) model of observables, perturbs parameters along slow trajectories, and updates parameters to maximize the largest eigenvalue of $A_o + A_{op}\left(\mathbb{I}_{d_o} \otimes \vec{p}\right)$ while preserving the mean observables $\bar{o}$. The method is demonstrated on a simple double-well system and a physically plausible five-box AMOC model, showing efficient, polynomial-scaling convergence and revealing how parameter changes can affect stability and hysteresis without altering observed means. This approach offers a practical pathway to identify stability biases and worst-case tipping scenarios in Earth system components, with broad applicability to EMICs and other complex dynamical systems.

Abstract

Models of complex dynamical systems like the Earth's climate often involve large numbers of uncertain parameters. Comprehensive exploration of the parameter space is typically prohibitive due to excessive computational costs, and systematic gradient-based parameter optimization is not feasible because such models are typically not differentiable. This is especially problematic in cases where the models describe highly nonlinear and possibly abrupt dynamics, where sensitivity to parameter changes is high. Components of Earth's climate system, such as the North Atlantic Overturning Circulation or the polar ice sheets, are at risk of undergoing critical transitions in response to anthropogenic climate change. Concerns have been raised that these Earth system components are too stable in state-of-the-art models. Here, we introduce a method for efficient, systematic, and objective calibration of dynamical complex system models, targeted at adjusting system stability. Given a number of physical or observational constraints, our method moves the system in a direction where the system loses or gains stability, guided by indicators of critical slowing down. In contrast to a brute force approach, where the computational cost exponentially increases with the number of parameters, our method scales polynomially and thus evades the curse of dimensionality. We successfully apply our method to a conceptual double-fold bifurcation model and a physically plausible reduced-order model of the global ocean circulation. Our method can efficiently adjust stability biases in a range of complex system models and help reveal potentially hidden instabilities and resulting state transitions in such models. These results have important implications, e.g., for Earth system models and ongoing efforts to improve their representation of key multistable Earth system components.

Figures (5)

• Figure 1: Main steps of our calibration method to automatically update parameter values that change the stability of a simulated system in a targeted way. o are observables constructed from the output of the complex model, p are the model parameters to be calibrated (in general, both are vectors). A vector-autoregressive model (VAR) is fitted to the data with the aim of varying p such that it maximizes the largest eigenvalue of matrix A (red), while keeping the time mean of o constant. The nature of vector and matrix multiplications are not specified here; for more detail on each step see Sect. \ref{['sec:Recipe']}.
• Figure 2: Application of the destabilization method to the double-well system. a) Evolution of the parameters $p_1$,...,$p_4$, and (b) evolution of the Jacobian $\lambda$ and the mean value $\bar{x}$ of $x$ during the iterative destabilization of the double-well system; note the small y-axis range for the latter. c) The right-hand side (RHS) of Eq. \ref{['eq: Double well ODE']} for the parameters before (red) and after (black) the destabilization. By design, observable $x$ is constrained to stay at $x=1$ (grey dotted line). The original stable state (negative slope of the red line) merges with an unstable one in a saddle-node bifurcation when the system is fully destabilized.
• Figure 3: Schematic illustration of the five-box ocean model (reproduced from Fig. 1a in Wood2019). The five boxes represent major ocean basins, and arrows represent freshwater fluxes between them. The variable q (purple arrows) is the property we consider as the 'observable' in this study. Freshwater hosing is applied by modifying the green surface fluxes $F_i$. Parameters $\gamma$, $\eta$, $K_N$, $K_S$, and $K_{IP}$ are the five parameters we calibrate.
• Figure 4: Destabilization of the five-box AMOC model with (a) the evolution of the parameters $p_1$,...,$p_5$, as well as (b) the Jacobian $\lambda$ and the mean value $\bar{x}$ of $x$ during the iterative destabilization process. (c) and (d): Hysteresis in the five-box AMOC model Wood2019 before (red) and after (black) destabilization, when performing the hosing experiment as described in Wood2019 to both systems. (c) Closeup of the hysteresis emphasizing the initial AMOC collapses due to the hosing. (d) Hysteresis curve over the full range of hosing. Note that the horizontal range of hysteresis is strongly increased by the parameter change. (e) 10,000 years of stationary time series when running the model with each parameter set, and without hosing (H=0).
• Figure 5: Computational cost of our calibration method to adjust stability for (a) the double well system, (b, black) the box model system, and (b, red) the optimized box model system. Both figures show how the length of the trajectory that is needed in a single iteration step in order to fulfill the accuracy condition (see Sect. \ref{['sec: Analysing the computational cost of the method']}), depends on the number of parameters under consideration. The small transparent data points correspond to (averaged) measurements for a given fixed parameter subset whereas the large data points display the average value taken over all parameter subsets, containing the same number of parameters. By fitting a linear function to the averaged data points in the log-log plot for each system respectively, we infer a polynomial dependence of the trajectory length on the number of parameters. We excluded the data points corresponding to one parameter being varied from the linear fits due to its unique significance for the generation process of the shown data.