Estimating Probability Distributions from Complex Models with Bifurcations: The Case of Ocean Circulation Collapse

Studying the uncertainty in computationally expensive models has required the development of specialized methods, including alternative sampling techniques and response surface approaches. However, existing techniques for response surface development break down when the model being studied exhibits discontinuities or bifurcations. One uncertain variable that exhibits this behavior is the thermohaline circulation (THC) as modeled in three dimensional general circulation models. This is a critical uncertainty for climate change policy studies. We investigate the development of a response surface for studying uncertainty in THC using the Deterministic Equivalent Modeling Method, a stochastic technique using expansions in orthogonal polynomials. We show that this approach is unable to reasonably approximate the model response. We demonstrate an alternative representation that accurately simulates the model’s response, using a basis function with properties similar to the model’s response over the uncertain parameter space. This indicates useful directions for future methodological improvements.

section_content: Estimating probability distributions of uncertain model outputs has long been a challenge for models requiring large amounts of computation time. A variety of methods have been developed for this problem, including specialized sampling methods [1] and constructing response surface approximation methods [2, 3]. One obstacle to using most response surface methods occurs when the model response exhibits discontinuities or bifurcations. 

An example of bifurcating behavior is the change in the circulation of the North Atlantic Ocean in long-term climate change projections. The thermohaline circulation (THC), or more formally, the zonally averaged meridional overturning circulation (MOC), refers to the circulation pattern of the North Atlantic Ocean in which warm surface water from the tropics travels northward, considerably warming mid and high latitudes in the Northem Hemisphere around the globe. This circulation is driven by deep water formation in the northern North Atlantic near Greenland, which is caused by the water becoming colder until it reaches a critical density that causes it to sink. As a possible consequence of climate change, it is hypothesized that warmer temperatures and increased freshwater runoff could prevent the water from reaching its critical threshold density, thus shutting off this circulation. 

The possibility of a collapse of the North Atlantic thermohaline circulation is one of the more severe potential impacts of climate change, and therefore is relevant to policy discussions [4]. A critical question, therefore, is: What is the probability of a THC collapse in the future?

One approach is to use simplified ocean models, which can reasonably be run for a large number of parametric assumptions [5]. However, for a more realistic representation of the ocean dynamics, one would ideally use a high-resolution three-dimensional ocean general circulation (GCM) model, coupled with a 3-D atmospheric GCM. A single simulation of several centuries with such models generally requires weeks to months on a supercomputer. Thus, even the small number of simulations (typically ~50 or more) required by methods such as Latin Hypercube Sampling [1] is prohibitive. Moreover, to inform policy, we need to know how the probability of a THC collapse will change with different policies, in addition to the reference case with no climate policy, requiring multiple sets of Monte Carlo simulations. 

To obtain the desired information from the more detailed models, some kind of reduced-form response surface model is needed that replicates the full 3-D dynamic behavior of the ocean, yet is simple enough to perform Monte Carlo on to obtain probability estimates. However, commonly used methods do not apply to a system with a bifurcation, and ocean circulation models are well-known to exhibit exactly this kind of behavior. 

In this article, we apply a commonly-used method for constructing optimal response-surface approximations for estimating the THC circulation from a 3-D ocean GCM. We will illustrate the challenges faced by this type of method, and demonstrate an alternative approach that is successful. The subsequent discussion frame directions for future research on more generalized approaches that can be applied to situations such as this one.
section_name: INTRODUCTION

section_content: Our coupled model of intermediate complexity consists of a three-dimensional ocean GCM [6] coupled to a zonally-averaged, statistical-dynamical atmospheric model [7], and a thermodynamic sea-ice model [8]. Further detail on the general coupled model can be found in [9, 10]. 

Our model's open passage through our idealized "Canadian Archipelago" plays an important role in the increased CO 2 simulations. Previous studies have speculated on the sensitivity of the ocean circulation and climate to freshwater discharge into the Arctic basin and subsequent flow into the Northern Atlantic [11] [12] [13] [14]. Our model employs a flexible river-routing scheme for anomalous runoff (as calculated in the atmospheric sub-component). In the southern hemisphere, for simplicity (and lacking a river network in this idealized topography) this runoff is distributed evenly over all ocean points. In the northern hemisphere, however, all anomalous runoff is diverted to the Arctic Ocean at 72-76°N between 96°and 260°in longitude. This diversion of anomalous runoff was necessary in order to achieve a complete collapse of the THC across a sizeable portion of our parameter phase space. Given this and other model idealizations, our model cannot be expected to give realistic information about when a collapse will occur. Rather, our goal is to study qualitatively how the collapse depends on the parameters. Such a study has previously only been carried out with two-dimensional models of the ocean basins [5]. 

For our climate change scenarios, the level of CO 2 is increased in the atmospheric model at a constant compounded rate for 100 years and then held constant at this resulting level. Thus, the rate of increase is proportional to the final change of forcing in the atmosphere. For the climate sensitivity parameter, different sensitivities are obtained by varying the strength of cloud feedback [7]. Varying the feedback allows the 2D atmospheric model to mimic the results of AGCMs with different sensitivities when coupled to a mixed layer ocean model, with a fixed ocean heat transport. Values of climate sensitivity shown throughout the article represent an equilibrium sensitivity of the atmospheric model coupled to a mixed layer ocean model for a doubling of CO 2 concentration. However, defined in this way, the climate sensitivity does not precisely match the climate sensitivity of the coupled climate model because of the interaction between the atmosphere and the dynamic ocean. 

We explore the uncertainty in the maximum overturning in the North Atlantic that is a consequence of uncertainty in two critical characteristics of climate system: the climate sensitivity and the rate of increase of CO 2 forcing. These uncertainties have previously been identified as primary determinants in ocean circulation changes [4, 5]. The assumed distribution for climate sensitivity, defined as the equilibrium warming resulting from a doubling of CO 2 concentrations, comes from [15], and is derived by updating expert priors with constraints from 20th-century observations. The probability density function (PDF) of the rate of CO 2 increase, driven primarily by anthropogenic emissions, is taken from [16], and is calculated from a Monte Carlo analysis of a macroeconomic model with uncertainty in economic growth rates and rates of energy efficiency improvement. Both PDFs are shown in Figure 1. The CO 2 forcing rate of increase is applied for the first 100 years of the simulation, and then CO 2 concentrations are held constant for the remaining 900 years of the of the simulation.
section_name: COUPLED CLIMATE MODEL DESCRIPTION

section_name: ALTERNATIVE METHODS FOR ESTIMATING PROBABILITIES

section_content: This section reviews the alternative methods for obtaining the uncertainty in an outcome from a deterministic computational model. Most simulation models are sufficiently complex that direct analytical solutions are not an option. The standard approach for uncertainty propagation is Monte Carlo simulation [17, 18], in which random samples are drawn from probability distributions of input parameters, the model is simulated for each random draw, and the frequency distribution of model outcomes provide the estimate of the probability distribution. The challenge to applying Monte Carlo comes when a model cannot be simulated thousands or tens of thousands of times. 

As discussed above, one approach is to use variance reduction methods for sampling from parameter distributions, so that fewer samples are needed for the estimated probability distribution of the outcome to converge. One popular and effective approach is stratified sampling, as in the Latin Hypercube Sampling (LHS) method [1, 19]. If the goal of the analysis is to estimate the probability of an extreme event, an alternative is to use Importance Sampling [20], which focuses on the low-probability region of interest. As mentioned previously, 3D ocean circulation models are likely to be too expensive for LHS, especially when separate Monte Carlo simulations must be performed for several different policy options. 

The other broad approach to estimating uncertainty from a computationallyintensive model is to construct a reduced-form model of the full model that produces a good approximation of the original model response with significantly less computation time. Reduced-form models can be further divided into two classes: theory-based or structural models and response surface approximations. Theory-based reduced-form models [21, 22] are simpler mathematical representations where the variables and equations still correspond to conceptual quantities and processes. This approach is primarily useful when transparency is critical for the reduced-form models behavior. The primary drawback is the extra time and effort required to develop a parsimonious closed-form model and the large number of runs of the original model to produce statistically acceptable parameter estimations. 

he other subclass of reduced form models is response surface approximations. In these methods, a mathematical representation of the full model's response surface is developed, focusing only on the uncertain parameters for the particular analysis and their relationship to the model outcome(s) of interest. There is a variety of methods for response surface approximation, ranging from simple linear models to more sophisticated techniques. The choices that distinguish between these methods are:

1. The choice of the basis function, the fundamental elements in the equation(s) to be fitted to the model responses, 2. The choice of which parameter values to evaluate the full model at and use to fit, and 3. The choice of solution method, given a set of data points from the model and a set of coefficients to solve for in the fitted equation(s). 

For example, standard linear approaches to response surface fitting [3] use firstor second-order polynomials of the uncertain parameters as a basis function, standard experimental design methods of choosing points for model evaluation, and minimize least-squared errors as the solution method to find the coefficients. 

An alternative method for response surface approximation is the Deterministic Equivalent Modeling Method (DEMM) [23, 24] -this method is also sometimes referred to as the Probabilistic Collocation method (PCM). It is equivalent to the Stochastic Response Surface Method (SRSM) developed by [2]. DEMM seeks to characterize the probabilistic response of the uncertain model output as an expansion in orthogonal polynomials. We describe DEMM in more detail below. 

There are several factors that determine which of the above methods is appropriate for any given situation, both the general class of approach (variance reduction vs. reduced-form model) and the particular choice (LHS vs. importance sampling). One important factor is the number of uncertain parameters under investigation. The number of simulations to obtain an accurate fit grows slowly for some methods (e. g., LHS) but expands rapidly for others (e. g., DEMM). Another critical factor is whether any prior information on the shape of the response within the range of uncertainty exists. Some methods are "black-box," no prior knowledge is required, while others (e. g., importance sampling) require some knowledge. DEMM is a good choice of method for estimating the uncertainty in the THC because: 1) it is black-box, requiring no prior knowledge of the shape; 2) the number of uncertain parameters is small (two); and 3) independent estimations of uncertainty are required for many different policy cases, which makes LHS infeasible.
section_name: Overview of Methods

section_content: Although any numerical computer model is itself deterministic, by positing uncertainty in a model parameter, the model's outputs become uncertain and thus can be thought of as a random variable. One useful representation for a random variable is an expansion of some family of orthogonal polynomials B N (x) with weighting coefficients a i :

where x is also a random variable of known distribution. Any family of orthogonal polynomials can be used, including Legendre, Laguerre, or Hermite. This expansion is sometimes referred to as a polynomial chaos expansion [25]. 

DEMM differs from the traditional approaches in all three steps that define a response-surface method. We first address the choice of the basis functions. Since a model output y is some function of its uncertain input parameter x, we can use information about the probability density of x to choose basis functions for the expansion. We can derive the set of orthogonal polynomials weighted by the density function of the parameter, according to the definition of orthogonal polynomials:

where dij

H i (x) and H j (x) are orthogonal polynomial functions of x of order I and j, P(x) is some weighting function, and C i is some constant (this constant is usually 1, and thus omitted, when the polynomials are normalized) in other words, the integral of the product of two orthogonal polynomials of different order is always 0. By using the probability density function of an input as the weighting function P(x), a set of orthogonal polynomials can be derived recursively. (The zero th order polynomial is always assumed to equal one. )

We next approach the method for estimating the weighting coefficients, a i. There is a class of methods designed for solving this problem known as the methods of weighted residuals (MWR) [26]. The residual at any realization x j of the random variable x, for some approximation í(x) of the function y(x) is simply the difference:

where R N ( ã, x j ) is the residual for an N-term expansion with weighting coefficients ã = {a 1, a 2,. . ., a N }. 

In general, MWR solves for N coefficients by solving the N relations:

Alternative schemes for MWR differ by the choice of the form of the weighting function, W j (x). Commonly used schemes include the least squares method, which chooses

, or Galerkin's method, which chooses W j (x) to be the derivatives of the approximation ¶ ¶ y a N j

. The difficulty with these schemes is that they require the explicit analytical form of the model in order to solve for the weighting coefficients. Because our goal is to approximate the uncertainty in a model output for any model, however complex, a method that allows the model to be treated as a "black-box" is preferable. This leads us to choose the collocation method, which uses the dirac delta function as the weighing function:

Since the integral of a function multiplied by a delta function is just the function evaluated at that point, solving (2) is equivalent to solving:

In other words, we simply solve for the set of a j such that the approximation is exactly equal to the model at N points, and thus only require the model solution at N points and not the explicit model equations. The final step in determining the polynomial chaos expansion to approximate the random variable is to choose the points x j at which we evaluate the "true" model y(x), in order to solve for the a i using equation (3). For this step, we borrow from the technique of Gaussian Quadrature, which uses the summation of orthogonal polynomials multiplied by weighting coefficients to approximate the solution of an integral. In Gaussian Quadrature, the optimal choice of abscissas at which to evaluate the function being integrated are the N roots of the Nth order orthogonal polynomial B N (x) [27]. Similarly in DEMM, to solve for the N coefficients in the expansion. 

we use the residual evaluated at the N roots of B N (x), the orthogonal polynomial one order higher than the highest order term. For multiple uncertain parameters, N roots are generated for each parameter to use as possible sample values. However, not all possible permutations of the N values for each parameter will necessarily be needed, depending on the number of terms in the expansion. Rather than combine sample values randomly, as in Latin Hypercube, we can use the probability density functions of the parameters to order the N possible values by likelihood. Then sample sets are formed by choosing permutations in decreasing order of joint probability, until the required number of sets has been formed. 

DEMM cannot find a sufficiently accurate approximation in every case. In particular, discontinuities in the response surface result in poor approximations. The approximation must be checked against model results at values of the uncertain inputs other than those used to solve for the coefficients. An optimal choice of points to check the approximation against the model is based on the roots of the next higher orthogonal polynomial than the one used to find points to solve at. The roots of the next higher order polynomial will always interleave the lower order roots [27], and so these will test the approximation at a maximal distance from the fit values while still spanning the highest probability regions. Moreover, if the expansion of order N results in an inaccurate fit, we already have the model results needed to solve the fit of order N + 1. Once the expansion for the probabilistic model response is solved and found to be reasonably accurate, the approximate probability density function of the response can be derived by applying Monte Carlo simulation to this expansion. 

DEMM and similar methods have been used successfully to explore the uncertainty in a variety of scientific, engineering, and economic modeling applications [23, 24, [28] [29] [30] [31] [32]. For many models, DEMM estimates multiple characteristics of the response distribution more efficiently than either modified sampling or traditional response surface approximation methods. DEMM's approach of representing the PDF of the uncertain response as an expansion of underlying PDFs, and of using probabilistic information in choosing the sample points for fitting the expansion, enable more efficient approximation of the overall response distribution relative to other methods.
section_name: The Deterministic Equivalent Modeling Method

section_name: RESULTS

section_content: The behavior of the maximum overturning for eight different parameter samples is shown in Figure 2. Note that for the first 100 years while CO 2 is increasing, the circulation slows in all cases, and does not collapse completely. But after several centuries the bifurcating behavior is apparent. For samples of either high climate sensitivity or rapid rate of CO 2 increase, ocean overturning continues to slow and shows no sign of rebounding within 1000 years. For samples with relatively low sensitivity and slow rate of CO 2 increase, the circulation recovers to close to present-day levels within a few centuries. Note that the transient behavior of the circulation in a simulation that does not recover (i. e., collapses) is continuous and smooth in the time dimension. The discontinuity is in the description of the circulation at one given point in time, for example in year 800, across all possible states of the world. The state of the circulation at some future time is the relevant outcome for policy studies.
section_name: Behavior of Ocean Model as Climate Sensitivity and CO 2 Forcing Changes

section_content: We first explore the application of DEMM to this problem. As described above, DEMM's use of orthogonal polynomials derived from the input PDFs is often superior to other response surface methods for non-linear surfaces, and has produced accurate estimates of probability distributions for a variety of applications including climate models. 

Figure 3 shows the sample points in parameter space used to fit and test, respectively, a 3rd order DEMM approximation. This requires 8 simulations of the coupled model used to solve for the coefficients (circles) and additional 10 simulations used to test the goodness of fit (crosses). Note that the sample points are designed to optimally span the joint density function of the input parameters. 

Before exploring response surfaces of ocean circulation strength, we first show the results for DEMM expansions of global mean surface air temperature (SAT) change. Third-order DEMM expansions for the parameter sets shown in Figure 3 result in approximations with sums of squared errors of less than 2% of the mean response value, accurately representing the response of the full climate model. Monte Carlo simulation is performed, drawing 10,000 random samples from the distributions for climate sensitivity and rate of CO 2 increase. The resulting PDFs of SAT change after 100 years and 1000 years are given in Figure 4. 

Unfortunately, unlike surface air temperature change, the DEMM expansions for maximum North Atlantic overturning have unacceptably large errors for all years beyond year 200 (Figure 5 ). This is not surprising, as the surfaces span the discontinuity between the region where the overturning recovers and the region where it does not (see Figure 3 ). A second approach is to still use DEMM, but to fit it piecewise on either side of the discontinuity. This requires first that we identify the threshold between the region in parameter space where circulation recovers and the region where it does not. A total of 62 simulations were performed and used to calculate the critical threshold for circulation recovery. We find that the threshold is best identified by s*r, the product of the sensitivity and the rate of CO 2 increase (Figure 6 ). When s*r < 1. 72, the circulation will recover, and when s*r > 1. 89, the circulation collapses and does not recover within 1000 years. 

Attempts to fit low-order DEMM approximations piecewise in each region of parameter space also fails to produce a reasonable representation of the ocean model's behavior. Figure 7 compares the best of the piecewise surfaces to the interpolated surface of the 62 GCM simulations. Monte Carlo simulations performed on DEMM approximations result in significant probability density for physically unrealistic values of maximum overturning below 0 and above 15 Sv. Further, piecewise fitting defeats the original purpose of selecting DEMM as a black-box method.
section_name: Response Surface Fits with Different Methods

section_content: To understand why any polynomial-based approximation will fail to yield a reasonable fit to the model, consider the shape of the model's response surface in Figure 7a. Note that the overturning strength, when fully recovered, levels out at around 10-12 Sv. Similarly, overturning strength, once fully collapsed, levels out at close to 0 Sv. Thus, the projection into either sensitivity or CO 2 rate parameter space, the maximum overturning function has the shape of a logistical S-curve. 

A low-order polynomial is unable to replicate this kind of S-curve shape, where function remains constant or approaches an asymptote above and below some critical values. As a demonstration, we apply DEMM to approximate the arctangent function, which exhibits this behavior. Treating arctan() as a black-box function, DEMM approximations are calculated, truncating terms at 3rd, 4th, 5th, and 7th order, respectively (Figure 8 ). Any low-order polynomial will have errors increasing exponentially in both directions beginning a short distance beyond the last model point used in the fit. A Monte Carlo with even low probability in these regions may yield large errors in the estimated PDF. Note that while a sufficiently large number of expansion terms in orthogonal polynomials could be found that would reasonably approximate this kind of function, it would require even more model simulations than one would need to directly simulate with Latin Hypercube Sampling, and thus would yield no advantage. 

The question becomes: is there an appropriate choice of basis function that WILL accurately replicate the model response across the parameter space? As described above, all response surface methods consist of a choice of basis function, a method of solving for coefficients, and a method of choosing points to evaluate the model for fitting. The problem here appears to be with the basis function choice. Having characterized the general shape of the response surface of the model, the ideal choice of basis function is one with the same logistical S-shaped curve. There are a number of functional forms with that shape from which to choose. One choice, from the example above, is the arctangent function. 

We use the 62 simulations of the ocean model to fit the function

where b is a shift parameter, a is an amplitude parameter, d is the inflection point parameter, s is the climate sensitivity, and r is the rate of CO 2 increase in % per year. Thus, we need to solve for three free parameters, b, a, and d, given a set of triplets (s, r, ovt). We solve for the parameters with ordinary least squares. The parameter values are given in Table 1. 

Fitting this equation produces a response surface that very closely resembles Figure 7a, and has extremely small errors of at most a few percent (Table 1 ). We then perform Monte Carlo simulation on this approximation, drawing 10,000 random samples from the distributions of climate sensitivity and forcing rate. The resulting PDF of overturning for year 1000 is shown in Figure 9. To estimate the probability of a collapse, we note that all parameter choices that recover have maximum circulations of 9 Sv or greater, while parameter choices that do not recover have maximum circulations of 8 Sv or less (Figure 7b ). By calculating the probability of a maximum overturning of 8 Sv or less, we estimate that the PROBABILITY OF THC COLLAPSE / 17 probability of a thermohaline circulation that collapses and does not recover within 1000 years is 13. 9%. 

This estimate is conditional on the assumed parameter distributions, but also importantly on the structural assumptions in the model. The true probability could be either higher or lower than this. More detailed studies are required with other coupled ocean-atmosphere GCMs for a range of assumptions to give better information on this likelihood. 

While the fit with 62 simulations achieves an acceptable level of accuracy, the goal is to develop a method with far fewer simulations if possible. We develop two more fits of equation 4 using the points chosen for a 3rd-order DEMM expansion. The first uses only the 8 parameter sets used by DEMM to fit, and the second uses all 18 parameter sets from Figure 3 which consist of the points used by DEMM to fit and the points used to check the fit. The 8 point fit has larger errors, but the 18 point fit is nearly as accurate as the 62 point fit. The estimate of the probability of THC collapse from the 18 point fit is 11. 6%, very close to the estimate from the 62 point fit. The results of a Monte Carlo on all three versions are shown in Figure 9.
section_name: A Successful Approximation Method

section_content: In this study, we have attempted to find a way of approximating the response of a coupled ocean-atmosphere general circulation model to changes in two critical uncertainties: climate sensitivity and the rate of CO 2 increase. In particular, our interest is in describing the relationship between these parameters and the likelihood of a collapse of the thermohaline circulation in the North Atlantic. Because this response is discontinuous with a bifurcation, it poses a particular challenge to developing an accurate reduced-form that is amenable to multiple rounds of Monte Carlo simulation. 

The solution to the methodological problem, while admittedly ad-hoc, points the way to new generalized techniques of response surface approximation. In the end, the obstacle to using existing methods was not so much the bifurcation, but the appropriate shape of the underlying basis functions. Although we leave the development of formal generalized methods to future work, needed improvements will be in the area of developing efficient methods for l) identifying the response surface shape characteristics; 2) choosing the appropriate basis functions for that shape, where the basis functions are chosen from a menu of options that include non-polynomial functions; and 3) identifying optimal points to sample the true model, given the choice of basis functions. 

This study also suggests a useful general approach for policy-focused studies of uncertainty in climate change. There is a hierarchy of complexity for climate models, ranging from simple box and 1-D models, to earth models of intermediate complexity (EMICs) which are often 2-D or 3-D with limited resolution, to full 3-D GCMs. One way to use this spectrum of available tools in studying the uncertainty in any climate change process is to study the process with an EMIC, develop an appropriate basis function for a response surface, and then conduct limited simulations with a full GCM to fit the response surface. This would be a hybrid approach between a theory-based and a response surface reduced-form model.
section_name: DISCUSSION

section_content: *This research was supported in part by the Methods and Models for Integrated Assessments Program of the National Science Foundation, Grant ATM-9909139 ; by the Office of Science (BER), U. S. Department of Energy, Grant Nos. DE-FG02-02ER63468 and DE-FG02-93ER61677 ; and by the MIT Joint Program on the Science and Policy of Global Change (JPSPGC). Financial support does not constitute an endorsement by NSF, DOE, or JPSPGC of the views expressed in this article.