Fitting to EOF1

The Imbrie model has four degrees of freedom: two for the timescales and two for insolation forcing. Imbrie and Imbrie [1980] searched for a best fit between their model output and a set of radiometrically datedδ¹⁸Orecords extending between zero to 130KY BP using a coarse grid search technique. They found optimal output for an ablation time-constant T2 = 10.6KY, accumulation time-constant T1 = 42.5KY, precessional phase φ = 16^◦, and insolation ratio of α = −2. It appears Imbrie and Imbrie [1980] define $ as the angle between vernal equinox and aphelion, rather than the more standard definition employing perihelion [e.g. Vernekar, 1972;Berger and Loutre, 1992], so that negative α and φ = 16^◦ indicates a July perihelion. This inference was checked by simulating Figure 3 in Imbrie and Imbrie [1980]. In this thesis, $ is always defined relative perihelion.

As the Imbrie model was originally tuned using data only over approximately the last 250KY, it is useful to repeat the exercise using the longer and better resolved EOF1 record. An exhaustive search of all combinations of plausible parameters can

easily be made because the Imbrie model has a small number of adjustable parame-ters. However, in anticipation of discussing models with more adjustable parameters, a simulated annealing search algorithm [e.g. Press et al., 1999] is instead used. Sim-ulated annealing is a Monte Carlo method which works in analogy with the slow cooling of a liquid from a hot and disordered state to a cool, crystallized state. If the adjustable parameters are considered molecules, and the misfit between model output and observations interpreted as energy, the cooled crystallized state corresponds to a minima in molecular energy or a local minima in the cost function.

The search domain for parameter values is restricted to one significant figure.

This gives a dramatic reduction in the number of possible values — if the search is restricted to values between one and a hundred, there are only nineteen possible values,{1,2...9,10,20...100}. Note, however, there are still an infinity of numbers with one significant figure between zero and one. Apart from making the parameter space easier to search, there are other advantages to restricting the search domain to one significant figure. The Imbrie model, and other models considered later, all represent drastic simplifications of the dynamics governing climate variability. If the behavior of the model is sensitive to the second or higher significant figures in a parameter value, it is less likely to represent a physically meaningful solution. One hopes to find a model which is robust to minor perturbations as this aids in identifying mechanisms likely to control climate. It should be kept in mind, however, that parameterizations with only one significant figure can still lead to results sensitive to minor perturbations

— they are only less likely to be sensitive. Furthermore, the climate system could itself be sensitive to minor perturbations; for instance, the glacial cycles could be chaotic in nature. The development and fitting of a simple model to observations is as much an art as a science, and there is no foolproof method for determining the adequacy of a given model or the accuracy of the parameterizations. Rather, insight into the adequacy of a given model requires careful consideration of the assumption, comparison of the results against independent observations as well as the results of other competing models.

The Imbrie model is initially fit to EOF1 over a search range with timescales T1

andT2 ranging from 1 to 120KY, a precessional phaseφbetween 10^◦ and 360^◦, and an obliquity to precession ratioαof -1 to 5. Values ofα <−1 are not permitted as these are redundant for values ofφ differing by±180^◦. As discussed in Chapter 4, only the last 650KY of EOF1 are used for fitting the models, as this is the portion which clearly

−600 −500 −400 −300 −200 −100

−600 −500 −400 −300 −200 −100

−15

−600 −500 −400 −300 −200 −100

−30

Figure 5-1: Results from the Imbrie model (red) after adjusting the free parameters to maximize the cross-correlations with EOF1 (black) over the last 650KY. EOF1 is scaled to the model output. a shows the best fit achieved when both obliquity and precession are included in the forcing, but requiring the ablation timescale to be longer than 1KY. b is the best fit when ablation timescales are permitted to be arbitrarily small, here a best fit is achieved forT₂=0.05KY.cshows the best fit when the model is forced only by obliquity variability, again using very rapid ablation.

The squared-cross-correlations between model results and EOF1 are shown at left.

Although difficult to see on this plot, the terminations in EOF1 occur over a period of roughly 10KY. At right are periodograms of the model results with vertical dashed lines indicates bands centered on 1/100, 1/70, 1/41, 1/29, and 1/23KY —- the bands of energy which exceed the 95% confidence level in EOF1.

shows 100KY variability; prior to 650KY BP, the glacial cycles have less low-frequency variability [e.g. Raymo and Nisancioglu, 2003]. To help ensure that the estimated parameters are the best global fit, rather than a localized feature in parameter space, the annealing algorithm is initiated at random locations in parameter space. For the Imbrie model, regardless of where the search is initialized, the same solution is consistently arrived at, suggesting it is globally the best fit for a one significant figure set of parameterizations.

The values which maximize the cross-correlation between EOF1 and the model results are T1 = 90KY, T2 = 10KY, φ = 60^◦, and α = −0.09 yielding a square-cross-correlation of 0.24. Relative to the original fit by Imbrie and Imbrie [1980], the reglaciation timescale derived from the fit to EOF1 is closer to observations, the obliquity variability is seen to be more important, but the squared-cross-correlation is nearly the same. The squared-cross-correlation between model and output will also be referred to as the fraction of variance described by the model. The significance of the variance described and robustness of the parameterization should be greater because the duration of EOF1 is three times that of the observations used byImbrie and Imbrie [1980]. The new fit is shown in Figure 5-1.

To gauge how sensitive the results of the Imbrie model are to the exact parame-terizations, a series of perturbation experiments are carried out. The squared-cross-correlation between model results and EOF1 is computed as a function of a single varying parameter while holding the others fixed at the optimal values determined using the simulated annealing method. Results are shown in Figure 5-2. As expected for a successful optimal fit of the model to observations, the largest cross-correlation is achieved when no perturbations are made to the model. The timescales are adjusted over a range of 10 and 200KY for T1, and 1 and 20KY for T2. The Imbrie model is most sensitive to making the accumulation timescale short (small values of T₁), dropping the squared-cross-correlation by a factor of four. The precessional phase of the insolation forcing was varied from 10^◦ to 360^◦ and the ratio of precession to obliq-uity energy was varied from zero (α = −1) to three (α = 2). Smallest correlations occur when the precession phase is changed by 180^◦, giving the equivalent of anoma-lously large Southern Hemisphere summer insolation. Overall, the cross-correlation between the Imbrie model results and EOF1 vary smoothly with changes in its pa-rameterization. Even for large changes, the Imbrie model retains some correlation with EOF1 showing that the Imbrie model is fairly robust to perturbations in the

parameterizations.

Figure 5-2: Model sensitivity to parameter perturbations. Starting from the best fit, the squared-cross-correlation between model results and EOF1 are computed over the last 650KY as a function of one of the parameters while holding the others fixed. From left to right are the T1 and T2 measured in KY, φ the precessional phase measured in degrees, and α which controls the ratio of obliquity to precession energy; minus one is no precession, zero is equal parts precession and obliquity.