In this paper, we present a new, statistical framework for understanding the dynamics of polarity reversals under time-dependent forcing. The methods are thought to be generally applicable to geodynamo models in any scenario involving a time-dependent forcing. We demonstrate this framework on a simple conceptual dynamo model, the Rikitake system, which is a classical model of polarity reversals, known for its chaotic dynamics. Here, we investigate, for the first time, a non-autonomous extension of this model, in which the coupling parameter changes in time, mimicking the gradual activation or strengthening of a geodynamo. Given the non-autonomous nature of the underlying dynamics, the typical behavior of the system can only be accurately characterized using ensemble-based statistical approaches (instead of time series-based ones), in line with the so-called snapshot attractor view of dynamical systems theory. Using ensemble simulations with a “spin-up” initialization, we analyze the statistical behavior of polarity reversals, specifically, we construct the time-dependent ensemble-based cumulative distribution functions (CDFs) of chron lengths. Results show that as the system evolves and the coupling increases, longer chrons become increasingly permitted, albeit with lower probability. This behavior can be interpreted as the analogy of very rare “superchron” events, supported by paleomagnetic evidence. We find that early-stage chron distributions follow Gaussian statistics, while later stages exhibit lognormal distributions, indicating a transition in the underlying dynamical behavior as the system evolves. An intriguing and new result is the simultaneous probabilistic coexistence of a plethora of chrons with a multitude of lengths, whose ensemble characterizes the system at each point in time. These findings highlight the importance of incorporating ensemble methods in advanced geodynamo modeling.