The major cause of sea-level change during ice ages is the exchange of water between ice and ocean and the planet’s dynamic response to the changing surface load. Inversion of ∼1,000 observations for the past 35,000 y from localities far from former ice margins has provided new constraints on the fluctuation of ice volume in this interval. Key results are: (i) a rapid final fall in global sea level of ∼40 m in <2,000 y at the onset of the glacial maximum ∼30,000 y before present (30 ka BP); (ii) a slow fall to −134 m from 29 to 21 ka BP with a maximum grounded ice volume of ∼52 × 10⁶ km³ greater than today; (iii) after an initial short duration rapid rise and a short interval of near-constant sea level, the main phase of deglaciation occurred from ∼16.5 ka BP to ∼8.2 ka BP at an average rate of rise of 12 m⋅ka⁻¹ punctuated by periods of greater, particularly at 14.5–14.0 ka BP at ≥40 mm⋅y⁻¹ (MWP-1A), and lesser, from 12.5 to 11.5 ka BP (Younger Dryas), rates; (iv) no evidence for a global MWP-1B event at ∼11.3 ka BP; and (v) a progressive decrease in the rate of rise from 8.2 ka to ∼2.5 ka BP, after which ocean volumes remained nearly constant until the renewed sea-level rise at 100–150 y ago, with no evidence of oscillations exceeding ∼15–20 cm in time intervals ≥200 y from 6 to 0.15 ka BP. Significance Several areas of earth science require knowledge of the fluctuations in sea level and ice volume through glacial cycles. These include understanding past ice sheets and providing boundary conditions for paleoclimate models, calibrating marine-sediment isotopic records, and providing the background signal for evaluating anthropogenic contributions to sea level. From ∼1,000 observations of sea level, allowing for isostatic and tectonic contributions, we have quantified the rise and fall in global ocean and ice volumes for the past 35,000 years. Of particular note is that during the ∼6,000 y up to the start of the recent rise ∼100−150 y ago, there is no evidence for global oscillations in sea level on time scales exceeding ∼200 y duration or 15−20 cm amplitude.

Summary: This paper presents a new derivative-free search method for finding models of acceptable data fit in a multidimensional parameter space. It falls into the same class of method as simulated annealing and genetic algorithms, which are commonly used for global optimization problems. The objective here is to find an ensemble of models that preferentially sample the good data-fitting regions of parameter space, rather than seeking a single optimal model. (A related paper deals with the quantitative appraisal of the ensemble.)The new search algorithm makes use of the geometrical constructs known as Voronoi cells to derive the search in parameter space. These are nearest neighbour regions defined under a suitable distance norm. The algorithm is conceptually simple, requires just two ‘tuning parameters’, and makes use of only the rank of a data fit criterion rather than the numerical value. In this way all difficulties associated with the scaling of a data misfit function are avoided, and any combination of data fit criteria can be used. It is also shown how Voronoi cells can be used to enhance any existing direct search algorithm, by intermittently replacing the forward modelling calculations with nearest neighbour calculations.The new direct search algorithm is illustrated with an application to a synthetic problem involving the inversion of receiver functions for crustal seismic structure. This is known to be a non-linear problem, where linearized inversion techniques suffer from a strong dependence on the starting solution. It is shown that the new algorithm produces a sophisticated type of ‘self-adaptive’ search behaviour, which to our knowledge has not been demonstrated in any previous technique of this kind.

Monte Carlo inversion techniques were first used by Earth scientists more than 30 years ago. Since that time they have been applied to a wide range of problems, from the inversion of free oscillation data for whole Earth seismic structure to studies at the meter‐scale lengths encountered in exploration seismology. This paper traces the development and application of Monte Carlo methods for inverse problems in the Earth sciences and in particular geophysics. The major developments in theory and application are traced from the earliest work of the Russian school and the pioneering studies in the west byPress[1968]to modern importance sampling and ensemble inference methods. The paper is divided into two parts. The first is a literature review, and the second is a summary of Monte Carlo techniques that are currently popular in geophysics. These include simulated annealing, genetic algorithms, and other importance sampling approaches. The objective is to act as both an introduction for newcomers to the field and a comprehensive reference source for researchers already familiar with Monte Carlo inversion. It is our hope that the paper will serve as a timely summary of an expanding and versatile methodology and also encourage applications to new areas of the Earth sciences.

Summary: Monte Carlo direct search methods, such as genetic algorithms, simulated annealing, etc., are often used to explore a finite-dimensional parameter space. They require the solving of the forward problem many times, that is, making predictions of observables from an earth model. The resulting ensemble of earth models represents all ‘information’ collected in the search process. Search techniques have been the subject of much study in geophysics; less attention is given to the appraisal of the ensemble. Often inferences are based on only a small subset of the ensemble, and sometimes a single member.This paper presents a new approach to the appraisal problem. To our knowledge this is the first time the general case has been addressed, that is, how to infer information from a complete ensemble, previously generated by any search method. The essence of the new approach is to use the information in the available ensemble to guide a resampling of the parameter space. This requires no further solving of the forward problem, but from the new ‘resampled’ ensemble we are able to obtain measures of resolution and trade-off in the model parameters, or any combinations of them.The new ensemble inference algorithm is illustrated on a highly non-linear wave-form inversion problem. It is shown how the computation time and memory requirements scale with the dimension of the parameter space and size of the ensemble. The method is highly parallel, and may easily be distributed across several computers. Since little is assumed about the initial ensemble of earth models, the technique is applicable to a wide variety of situations. For example, it may be applied to perform ‘error analysis’ using the ensemble generated by a genetic algorithm, or any other direct search method.

Summary. The reversible jump algorithm is a statistical method for Bayesian inference with a variable number of unknowns. Here, we apply this method to the sei