Mark A Cunningham
email: mark.cunningham@nrc.gov

Mark A Cunningham Senior

Risk Advisor

Considerations of Uncertainty in Regulatory Decision Making

Keywords: risk assessment, nuclear safety, performance-based regulation, safety goals, probabilistic risk assessment

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Nuclear power plants in the United States historically have been licensed using Part 50 of Title 10 of the Code of Federal Regulations [1]. Implementation of Part 50 has been achieved, for the most part, using deterministic methods and acceptance criteria. These deterministic methods and acceptance criteria were established originally in the 1960's and 1970's and were intentionally conservative in recognition of uncertainties in both routine operations and potential accident conditions. Concepts used included:

-A set of "design basis" accidents (DBAs) that was intended to envelope conditions from a credible set of events, -A "single failure criterion," a qualitative approach to ensure that systems used to mitigate accidents were highly reliable, -A "defense in depth" philosophy that introduced barriers between radioactive material and workers and the public, and -Inclusion of safety margins, a traditional engineering approach for ensuring a robust design.

In 1975, NRC published its first probabilistic risk assessment (PRA) which examined two reactors designed using Part 50 from a different, more realistic, perspective [START_REF]Reactor Safety Study: An Assessment of Accident Risks in US Commercial Nuclear Power Plants[END_REF]. This study considered a broader set of possible accidents (relative to the DBAs) and estimated their occurrence frequencies, evaluated system reliability quantitatively, estimated the potential public health consequences, and measured, in effect, the effectiveness of the included defense in depth and safety margins. Public health risk was estimated, including an estimate of the uncertainties in this risk. The value of PRA, highlighted by the investigations of the 1979 accident at Three Mile Island, led to the performance of a number of PRAs on other plants and, in the late 1980's, an examination by all operating nuclear power plants for potential vulnerabilities using PRA techniques [START_REF]Individual Plant Examination for Severe Accident Vulnerabilities[END_REF].

Commission Policy on Use of Risk Assessment

By the mid-1990s, the nuclear industry had gained considerable experience with implementation of Part 50 and the results of PRAs. In 1995, NRC published a statement describing, among other things, the relationship between these two views of reactor safety [START_REF]USNRC: Use of Probabilistic Risk Assessment in Methods in Nuclear Regulatory Activities, Final Policy Statement[END_REF]. This policy statement directed the NRC staff to increase the use of PRA, setting a course for staff activities that has resulted in a significant expansion in its use.

The Commission's policy statement summarizes the value of risk assessment as follows:

The NRC has generally regulated the use of nuclear material based on deterministic approaches. Deterministic approaches to regulation consider a set of challenges to safety and determine how those challenges should be mitigated. A probabilistic approach to regulation enhances and extends this traditional, deterministic approach, by: (1) allowing consideration of a broader set of potential challenges to safety, (2) providing a logical means for prioritizing these challenges based on risk significance, and (3) allowing consideration of a broader set of resources to defend against these challenges.

With this perspective on the value, and relative roles, of traditional and risk methods, the Commission established the following as its policy:

Increase use of PRA technology in all regulatory matters to the extent supported by the state-of-the-art in PRA methods and data and in a way that complements the deterministic approach and supports the traditional defense-in-depth philosophy.

Use PRA, where practical within the bounds of the state-of-the-art, to reduce unnecessary conservatism in current regulatory requirements, regulatory guides, license commitments, and staff positions and to support proposals for additional regulatory requirements in accordance with 10 CFR 50.109 (Backfit Rule).

PRAs used in regulatory decisions should be as realistic as practicable and supporting data should be publicly available. Safety goals and subsidiary numerical objectives are to be used with appropriate consideration of uncertainties in making regulatory judgments on the need for new generic requirements.

NRC has made significant progress in implementing this policy in the past 15 years. In the late 1990's-early 2000's time frame, the NRC staff undertook a number of initiatives to better incorporate risk insights and performance considerations into its regulatory programs. In addition to regulatory changes, the NRC worked with other organizations (e.g., the American Society of Mechanical Engineers and the National Institute of Standards and Technology) to improve the technical infrastructure underlying risk assessments. These improvements included the development of consensus standards [START_REF]ASME/ANS RA-Sa-2009, Standard for Level 1/Large Early Release Frequency Probabilistic Risk Assessment for Nuclear Power Plant Applications, Addendum A to RA-S-2008[END_REF], the development of new methods [21], and performing research including developing better computational methods to validate these new assessment techniques [START_REF]NUREG-1824, Verification and Validation of Selected Fire Models for Nuclear Power Plant Applications[END_REF]. These initiatives resulted in fundamental changes to how the NRC conducts its licensing, inspection and rulemaking programs.

NRC's Commission has also directed the NRC staff to solicit input from industry and other stakeholders on performance-based initiatives, including areas that are not amenable to risk-informed approaches, to supplement the NRC's traditional deterministic system of licensing and oversight. It should be noted that deterministic 1 and prescriptive 2 regulatory requirements were based mostly on experience, testing programs and expert judgment, considering factors such as safety margins and the principle of defense-in-depth. These requirements are viewed as being successful in establishing and maintaining adequate safety margins for NRC-licensed activities. The NRC has recognized that deterministic and prescriptive approaches can limit the flexibility of both the regulated industries and the NRC to respond to lessons learned from operating experience and support the adoption of improved designs or processes.

The NRC has as one of its primary safety goal strategies the use of sound science and state-of-the-art methods to establish, where appropriate, risk-informed and performance-based regulations. The NRC issued a paper [START_REF]USNRC: White Paper on Risk-Informed and Performance-Based Regulation[END_REF] to define the 1 A deterministic approach to regulation establishes requirements for engineering margin and for quality assurance in design, manufacture, and construction. In addition, it assumes that adverse conditions can exist and establishes a specific set of design basis events and related acceptance criteria for specific systems, structures, and components based on historical information, engineering judgment, and desired safety margins. An example is a defined load on a structure (e.g., from wind, seismic events, or pipe rupture) and an engineering analysis to show that the structure maintains its integrity. 2 A prescriptive requirement specifies particular features, actions, or programmatic elements to be included in the design or process, as the means for achieving a desired objective. An example is a requirement for specific equipment (e.g., pumps, valves, heat exchangers) needed to accomplish a particular function (e.g., remove a defined heat load).

terminology and expectations for evaluating and implementing the initiatives related to risk-informed, performance-based approaches. That paper defines a performance-based approach as follows:

A performance-based regulatory approach is one that establishes performance and results as the primary basis for regulatory decision-making, and incorporates the following attributes: 1. measurable (or calculable) parameters (i.e., direct measurement of the physical parameter of interest or of related parameters that can be used to calculate the parameter of interest) exist to monitor system, including facility and licensee, performance, 2. objective criteria to assess performance are established based on risk insights, deterministic analyzes and/or performance history, 3. licensees have flexibility to determine how to meet the established performance criteria in ways that will encourage and reward improved outcomes; and 4. a framework exists in which the failure to meet a performance criterion, while undesirable, will not in and of itself constitute or result in an immediate safety concern. 3 Performance-based approaches can be pursued either independently or in combination with risk-informed approaches. After the paper's issuance, NRC continued to make progress on developing policies and guidance related to performance-based approaches and subsequently issued guidance documents [START_REF]USNRC: SECY-2000-0191, High-Level Guidelines for Performance-Based Activities[END_REF] [START_REF]USNRC: Guidance for Performance-Based Regulation[END_REF].

Perhaps the most significant programmatic adoption of risk-informed and performance-based considerations in the reactor area took place with implementation of the "reactor oversight process" (ROP) [START_REF]Recommendations for Reactor Oversight Process Improvements[END_REF]. The ROP, intended to focus agency reactor inspection resources in the most risk-significant areas of reactor operation, replaced the previous program with explicit consideration of risk and performance considerations. The normal "baseline" inspection program is focused on the more risk-important areas of plant operations. In addition, events or conditions at plants are assessed for significance using probabilistic risk models. The results of such assessments are used to direct additional oversight to plants with more significant findings.

A more recent reactor initiative that adopts a risk-informed and performancebased approach relates to fire protection, in which standards from the National Fire Protection Association (NFPA-805) were incorporated into NRC's regulation in 10CFR50.48(c) [START_REF][END_REF]. This regulation provides deterministic requirements that are very similar to those in NRC's traditional fire protection regulations, and also includes performance-based methods for evaluating plant configurations that provide a comparable and equivalent level of safety intended by the 3 Using the previous example (footnote 2), a performance-based approach might provide additional flexibility to a licensee on plant equipment and configurations used to accomplish a safety function (e.g., removing a heat load), but the performance criteria could not be the actual loss of a safety function that would result in the release of radioactive materials.

conservative deterministic requirements. The performance-based methods allow engineering analyzes to demonstrate that the changes in overall plant risk that result from these plant configurations is acceptably small and that fire protection defense-in-depth is maintained. 4 Defense-in-depth as applied to fire protection means that an appropriate balance is maintained between:

1. preventing fires from starting; 2. timely detection and extinguishing of fires that might occur; and 3. protection of systems, structures and components important to safety from a fire that is not promptly extinguished.

The adoption of NFPA 805, which is voluntary on the part of reactor licensees, provides a licensee with flexibility regarding how to implement its fire protection program while maintaining an acceptable level of fire safety.

In parallel, the NRC staff was incorporating risk insights into other regulatory areas. In the materials area, a staff document [START_REF]NUREG-1556, Consolidated Guidance about Materials Licenses[END_REF] was developed in the late 1990's to pull together into one place the various guidance documents written over the years for the wide variety of materials licensees. These documents allow license applicants to find the applicable regulations, guidance and acceptance criteria used in granting a materials license. Operational experience (performance) and risk insights guided the development of these documents. Over time the guidance has been revised to further incorporate risk insights, performance considerations and changing technology. A new revision to the series is under development to address security and other issues.

The materials inspection program was fundamentally revised in 2001both in terms of approach and frequency -in the Phase II Byproduct Material Review [13]. The inspection approach was modified to emphasize licensee knowledge and performance of NRC-licensed activities over document review. Inspectors now review a licensee's program against focus areas that reflect those attributes which are considered to be most risk-significant. If a licensee's performance against a given focus element during the inspection is considered to be acceptable, the inspector moves on to the next focus element. Performance concerns or questions lead an inspector to go deeper into that area. In addition, inspection frequencies were revised based on risk insights as well as licensee performance over time.

Developing a Strategic Vision

In early 2011, an NRC staff group was established [START_REF]Charter for Task Force for Assessment of Options for More Holistic Risk-Informed, Performance-Based Regulatory Approach[END_REF] to, in effect, reflect on the past 15 years of experience and to develop a "strategic vision and options for adopting a more comprehensive and holistic risk-informed, performance-based regulatory approach for reactors, materials, waste, fuel cycle, and transportation that would continue to ensure the safe and secure use of nuclear material." This group was established by Chairman Gregory Jaczko and is being headed by Commissioner George Apostolakis.

When this group was established, the agency was dealing with a number of challenging regulatory issues. In several instances, the technical understanding of the issue was relatively poor, with important uncertainties in key topics. Just after the group was established, the Fukushima nuclear power plants in Japan experienced a large earthquake and tsunami. The effects of these events, and their implications to the safety and regulation of US nuclear power plants, introduced additional challenges. Some of these challenges are discussed in more detail below.

Challenging Regulatory Topics

Performance of Emergency Core Cooling Systems. In 2004, NRC issued a communication [START_REF]USNRC: Generic Letter 2004-02, Potential Impact of Debris Blockage on Emergency Recirculation during Design Basis Accidents at Pressurized Water Reactors[END_REF] to all power reactor licensees requesting that each perform an evaluation to address: the identified potential susceptibility of pressurized-water reactor (PWR) recirculation sump screens to debris blockage during design basis accidents requiring recirculation operation of [emergency core cooling system] ECCS or [containment spray system] CSS and on the potential for additional adverse effects due to debris blockage of flowpaths necessary for ECCS and CSS recirculation and containment drainage.

Since that time, power reactor licensees have made modifications, as necessary, to address the 2004 issue. However, other related issues have also been identified, including the possible effects of chemical additives on the debris characteristics, potentially worsening the blockage potential, and the potential effects of debris entering the reactor core region and causing blockages there.

In this example, decision makers are provided technical information having uncertainties in several key technical areas, including the effect on debris accumulation, and possible cooling system blockage, of chemical interactions, as well the effect of debris potentially entering and blocking portions of the reactor core area. Experimental work is underway to provide additional data in both areas, but the applicability of the results to the different reactor types also complicates decision making.

NRC staff addressed the complex technical and regulatory issues in a paper [START_REF]Closure Options for Generic Safety Issue-191, Assessment of Debris Accumulation on Pressurized Water Reactor Sump Performance[END_REF] providing options to the NRC's Commission for decision. The paper recommended that the ongoing staff approach be continued, which included a determination "whether, given the conservatisms, nonconservatisms, and/or uncertainties in the various review areas, the licensee has demonstrated adequate strainer performance and therefore compliance with the regulations." The schedule for this considered the risk of different types of accidents, such that it would "address any outstanding issues associated with more likely and risksignificant smaller loss of coolant accidents (LOCAs) (14 inches and below) in the short term, but would allow more time to address issues associated with the low-likelihood larger break LOCAs (above 14 inches). In this way, the more risksignificant issues would be closed quickly, and licensees would have the flexibility to reduce the impact (cost and dose) of addressing the less risk-significant LOCAs through planning, testing, or refined analyzes." The paper also recommended a second approach to be used in combination that "would provide more flexibility to licensees for addressing larger LOCAs than is currently permitted" using existing guidance and that would "likely reduce the scope of modifications needed to address [the issue] for some plants and would be consistent with agency policy regarding risk-informed regulation." NRC's Commission subsequently approved the recommended approach [START_REF]Assessment of Debris Accumulation on Pressurized Water Reactor Sump Performance[END_REF] with "comments and clarifications." Earthquake Frequencies in the Central and Eastern United States. In 2010, NRC staff completed an assessment [START_REF]Safety/Risk Assessment Results for Generic Issue 199[END_REF] of new information related to potential earthquakes in the central and eastern United States (that part of the United States east of the Rocky Mountains). In some cases, the new information indicated that estimated frequencies of earthquakes increased relative to previous estimates. Not surprisingly, these estimates had considerable uncertainty.

Since that time, NRC staff have been working to determine what, if any, actions need to be taken by power reactor licensees. A progressive screening approach is being considered which would include comparisons with deterministic information used in the initial design and, if necessary, two alternative seismic risk assessment methods.

One Impact of the Fukushima Accident. In March 2011, the Fukushima nuclear power station in Japan experienced a very large earthquake and subsequent tsunami. The resulting damage is described in a number of documents, and is summarized as follows in an NRC report [START_REF]USNRC: Recommendations for Enhance Reactor Safety in the 21st Century[END_REF]:

As a result of the earthquake, all of the operating units appeared to experience a normal reactor trip within the capability of the safety design of the plants. The three operating units at Fukushima Dai-ichi automatically shut down, apparently inserting all control rods into the reactor. As a result of the earthquake, offsite power was lost to the entire facility. The emergency diesel generators started at all six units providing alternating current (ac) electrical power to critical systems at each unit, and the facility response to the seismic event appears to have been normal.

Approximately 40 minutes following the earthquake and shutdown of the operating units, the first large tsunami wave inundated the site followed by multiple additional waves. The estimated height of the tsunami exceeded the site design protection from tsunamis by approximately 8 meters (27 feet). The tsunami resulted in extensive damage to site facilities and a complete loss of ac electrical power at Units 1 through 5, a condition known as station blackout (SBO). Unit 6 retained the function of one of the diesel generators.

The operators were faced with a catastrophic, unprecedented emergency situation. They had to work in nearly total darkness with very limited instrumentation and control systems. The operators were able to successfully cross-tie the single operating Unit 6 air-cooled diesel generator to provide sufficient ac electrical power for Units 5 and 6 to place and maintain those units in a safe shutdown condition, eventually achieving and maintaining cold shutdown.

Despite the actions of the operators following the earthquake and tsunami, cooling was lost to the fuel in the Unit 1 reactor after several hours, the Unit 2 reactor after about 71 hours, and the Unit 3 reactor after about 36 hours, resulting in damage to the nuclear fuel shortly after the loss of cooling. Without ac power, the plants were likely relying on batteries and turbine-driven and diesel-driven pumps. The operators were likely implementing their severe accident management program to maintain core cooling functions well beyond the normal capacity of the station batteries. Without the response of offsite assistance, which appears to have been hampered by the devastation in the area, among other factors, each unit eventually lost the capability to further extend cooling of the reactor cores.

The current condition of the Unit 1, 2, and 3 reactors is relatively static, but those units have yet to achieve a stable, cold shutdown condition. Units 1, 2, 3, and 4 also experienced explosions further damaging the facilities and primary and secondary containment structures. The Unit 1, 2, and 3 explosions were caused by the buildup of hydrogen gas within primary containment produced during fuel damage in the reactor and subsequent movement of that hydrogen gas from the drywell into the secondary containment. The source of the explosive gases causing the Unit 4 explosion remains unclear. In addition, the operators were unable to monitor the condition of and restore normal cooling flow to the Unit 1, 2, 3, and 4 spent fuel pools.

From a decision making perspective, this accident raises issues with respect to the ability to predict the likelihood of very large earthquakes and to decide how large of an earthquake should be considered sufficiently likely that it must be considered in a nuclear power plant's design. In addition, this example introduces other important uncertainties. One key additional issue is the effectiveness of emergency response, relied upon in nuclear safety to help ensure that the potentially affected population near nuclear power plants would not be exposed to large amounts of radioactive material.

An NRC group provided recommendations to the NRC Commission [START_REF]USNRC: Recommendations for Enhance Reactor Safety in the 21st Century[END_REF] that included some that reflect on how to address issues with considerable uncer-tainty, and the relative role of deterministic and risk assessment methods. These recommendations are currently under evaluation.

Some Issues to Resolve

These example issues raise some important questions, including:

-What are the relative roles of traditional engineering approaches and risk assessment? -Are there better ways to collect and analyze new information? -Should performance based approaches be used to a greater extent to better reflect such new information? -How should decision makers include consideration of very unlikely events that could result in very large consequences?

As noted above, an NRC staff group is now considering a strategic vision for a more comprehensive and holistic risk-informed, performance-based regulatory approach. This group is considering questions such as these; it expects to complete its report in the spring, 2012.

Building upon the guidance in Regulatory Guide 1.174[START_REF]Regulatory Guide 1.174, An Approach for Using Probabilistic Risk Assessment in Risk-Informed Decisions on Plant-Specific Changes to the Licensing Basis[END_REF], "An Approach for Using Probabilistic Risk Assessment in Risk-Informed Decisions on Plant-Specific Changes to the Licensing Basis," Regulatory Guide 1.205[START_REF]Regulatory Guide 1.205, Risk-Informed, Performance-Based Fire Protection for Existing Light-Water Nuclear Power Plants[END_REF], "Risk-Informed, Performance-Based Fire Protection for Existing Light-Water Nuclear Power Plants," states: "Prior NRC review and approval is not required for individual changes that result in a risk increase less than 1 × 10 -7 per year for CDF (core damage frequency) and less than 1 × 10 -8 per year for LERF (large early release frequency -a measure of potential offsite health consequences). The proposed change must also be consistent with the defense-in-depth philosophy and must maintain sufficient safety margins. The change may be implemented following completion of the plant change evaluation."

DISCUSSION

Speaker: Mark Cunningham

Philip Starhill : Presumably the landscape of low probability events is much more diverse than that of higher probability events. How does that impact the utility of deterministic modeling of low probability events? Mark Cunningham : The estimation of very low probability events, such as very large earthquakes or floods, or the failure probability of very large pipes used in nuclear power plants, is a particularly challenging aspect of nuclear safety analyzes. Historically, nuclear safety organizations such as NRC have made conservative deterministic statements on what magnitudes of earthquakes or sizes of pipe breaks needed to be considered in plant designs. These statements included, of course, some consideration on the credibility of the specified earthquake or pipe break, which was tied to some perception at the time of their probability. Modern risk assessment methods supplement these statements by including more quantitative probability estimates. Expert interpretation of available information, and the translation of such interpretations in statements of probability, is used for low probability events. This approach, especially when a number of experts are used, has been found to be an acceptable method when sufficient data are not otherwise available. Van Snyder : A numeric limit is established for risk of death within a specified distance due to nuclear power plant operations. Has there been a quantification, and limit, of total system risk to the public at large, and a comparison of that total system risk to total system risk from alternatives, especially coal? Mark Cunningham : In 1986, NRC established "safety goals" for the operation of nuclear power plants. These goals were intended to ensure that public risk from nuclear power plant operations was a small fraction of other risks, including the risks from other forms of electric power generation. Risk assessments performed by NRC indicate that nuclear power plants meet the established safety goals. NRC has not made comparisons between, for example, risks from electric power generation from coal. Some such studies were performed, however, by staff at Oak Ridge National Laboratory a number of years ago.

William Kahan : If I remember correctly, burning coal releases Radon; and coal ash is very weakly radioactive, but there is a vast amount of it to be stored. When the totals of radioactivity liberated by coal mining are added up, do they bring the coal industry under the purview of the Nuclear Regulatory Commission? Should they? Mark Cunningham : The coal industry is not subject to regulatory review by NRC. Since its establishment in 1975 by Congress, the scope of NRC's statutory authority has been discussed by Congress to possibly include, for example, Department of Energy facilities using nuclear materials. To my knowledge, these discussions have not considered the coal industry.

Jeffrey Fong : My question is on stability of operation. After deregulation, electricity is being treated as a commodity purchased every six minutes. This creates instability and pressure on operators. Does NRC have new regulations to address this type of operational instability? Mark Cunningham : When deregulation was beginning to occur, NRC recognized that nuclear power plant operations could be impacted by deregulation of the electrical energy market, and reassessed its regulatory approach. To my knowledge, no new regulations were established specifically to address the possible effects of deregulation.

Jeffrey Fong : Do risk assessment standards include quantitative estimates of human errors during all phases of nuclear reactor operations? Mark Cunningham : Yes, risk assessments do include quantitative consideration of human errors.


Peter Challenor
email: p.challenor@noc.ac.uk

Using Emulators to Estimate Uncertainty in Complex Models

Keywords: emulator, Gaussian process, Bayes linear, sensitivity, uncertainty, calibration

The Managing Uncertainty in Complex Models project has been developing methods for estimating uncertainty in complex models using emulators. Emulators are statistical descriptions of our beliefs about the models (or simulators). They can also be thought of as interpolators of simulator outputs between previous runs. Because they are quick to run, emulators can be used to carry out calculations that would otherwise require large numbers of simulator runs, for example Monte Carlo uncertainty calculations. Both Gaussian and Bayes Linear emulators will be explained and examples given. One of the outputs of the MUCM project is the MUCM toolkit, an on-line recipe book for emulator based methods. Using the toolkit as our basis we will illustrate the breadth of applications that can be addressed by emulator methodology and detail some of the methodology. We will cover sensitivity and uncertainty analysis and describe in less detail other aspects such as how emulators can also be used to calibrate complex computer simulators and how they can be modified for use with stochastic simulators.

Introduction

The increase in computing power over the last few decades has led to an explosion in the use of complex computer codes in both science and engineering. Almost every area of science and engineering now uses complex numerical simulators to solve problems that could not have been tackled only a few years ago. Examples include engineering [START_REF] Bates | Achieving robust design from computer simulations[END_REF], climate science, [START_REF] Collins | The Community Climate System Model: CCSM3[END_REF], and oceanography [START_REF] Maltrud | An eddy resolving global 1/10 • ocean simulation[END_REF]. At the present time most of these simulators are deterministic but there is an increasing use of stochastic simulators as well [START_REF] Wilkinson | Stochastic Modelling for Systems Biology[END_REF]. These computer codes comprise many thousands (or even millions) of lines of code and take long times to compute even on the fastest supercomputers available today. Questions we would like to ask of the simulators include: for a given uncertainty in the inputs to the simulator what is the uncertainty in the outputs; which inputs have the most effect on the outputs; are all the inputs important; how can we relate the simulator to reality? In this paper we look at how emulators can be used to address these problems. For a more detailed examination of the background and theory of what is discussed here see [START_REF] Goldstein | Bayesian analysis for complex physical systems modelled by computer simulations: current status and future challenges[END_REF].

Uncertainty

Assuming for the moment that our simulators are deterministic rather than stochastic where does the uncertainty come from? We know that the predictions we make are not exact. (If you are tempted to think of your simulator as perfect, how much you would be prepared to bet on its result being the same as a physical experiment?) Some of the 'error' may be numerical, running on different computer architectures or in different precisions will give different answers. But for well written code these differences will be small. In general the uncertainty in our simulator outputs comes from two sources: from uncertainty in its inputs and from uncertainty in its structure. By inputs we mean all the external inputs to the simulator; these include the initial conditions, boundary conditions and parameters. Although these may have very different properties, for example initial conditions are often spatial fields while the parameters are normally collections of single numbers, we will generally treat them the same. Similarly we take a Bayesian approach and make no distinction between aleatoric uncertainty, arising from genuine randomness, and epistemic uncertainty which is a measure of our ignorance.

The input uncertainty can be thought of as the internal uncertainty within the simulator. The structural uncertainty is how our particular simulator relates to other simulators and more importantly to reality. This is discussed further below and in [START_REF] Goldstein | Bayesian analysis for complex physical systems modelled by computer simulations: current status and future challenges[END_REF].

Quantifying Uncertainty

Before we consider structural uncertainty and the relationship between simulators and reality lets first consider the input uncertainty. We take a Bayesian approach to the problem, but it can be reformulated in terms of frequentist statistics [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. The formulae generally stay the same but the justification is different.

We have some variables (y) that we are interested in and which we will call outputs. These are related to another set of variables (x), which we will call inputs. We can find the values of the outputs corresponding to the values of the inputs by running a complex computer program which we will denote by f . Mathematically we can write y = f (x)

We will assume for now that f is deterministic, so if we run f with the same set of inputs we will obtain the same set of outputs. We assume that any numerical error is small enough to be ignored. Some or all of our inputs might be uncertain. This might arise from genuine randomness, which we call aleatoric uncertainty; or the uncertainty could be a result of lack of knowledge, called epistemic uncertainty. The distinction between these two forms on uncertainty is not as clear as might appear. For example consider a coin toss. While the coin is in the air our uncertainty is aleatoric, but once it has landed, but is hidden from us, the uncertainty becomes epistemic. Because of this interchangeability we treat all sources of uncertainty the same and describe them with probability density functions (π(x)). If we want to know the uncertainty on y given π(x) we need to calculate the transformation of π(x) induced by f . If f is linear this is an easy problem and can be solved analytically, but the complex codes we are interested in are non-linear and an analytical solution does not exist. The naive solution is to use Monte Carlo methods: draw a sample from π(x), x i , propagate this through the program to produce y i = f (x i ). The resulting y i are then a sample from π(y) from which we can estimate π(y). Such methods are effective but large samples are needed, particularly for high dimensional x and y. This means that Monte Carlo methods are not viable for expensive computer codes. Our solution to this problem is to build a fast approximation to the full numerical simulator. This is known as an emulator. We not only want fast approximators, we want fast approximators that estimate their own error. It cannot be stressed too much how useful it is having knowledge of the uncertainty in the emulator estimate. When we come to validate our emulators it is invaluable as it gives us a measure against which we can gauge how far the estimator may be from the truth. Similarly it is very helpful in sequential design, where we can put the next point at the most uncertain current point and when we come to estimate the uncertainty of the outputs we need to include the uncertainty arising from the process of emulation.

Gaussian Processes and Emulators

Our requirement that we have an emulator that is fast and contains an estimate of its own uncertainty is satisfied by the Gaussian process, [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. Gaussian processes are very adaptable stochastic processes which can be used to fit nonlinear data. A Gaussian process (GP) is the infinite dimensional analogue of a Gaussian distribution. It is defined by a mean function and a covariance function. The mean function gives the expected value of the GP at any point. The covariance function then gives the covariance between any two points. The form of the covariance function dictates how 'smooth' a realisation of the GP is. For example if we use the 'exponential' form of the covariance function

cov(x 1 , x 2 ) = σ 2 e -|x1-x2| (1)
we get Brownian motion and any realisation does not possess any derivatives.

On the other hand if we use the squared exponential covariance

cov(x 1 , x 2 ) = σ 2 e -(x1-x2) 2 (2)
then all derivatives exist and we get a very smooth set of realisations. There are other forms of the covariance function, see [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] for further details. For uncertainty quantification a limited selection of covariance functions tend to be used. We rarely believe our simulators do not possess any derivatives so we tend not use the exponential form. The squared exponential is probably the most used but this can give rise to numerical problems (see below). Other options include the Matern and generalised exponential form

σ 2 exp(-||x 1 -x 2 || α ).
The smoothness of the GP is dictated by its covariance function but its large scale properties come from the mean function. In general a linear model is used as the mean function

µ(x) = h(x) T β (3)
where the h(x) are basis functions and the β's are coefficients. The h(x) can be any basis functions but are usually taken to be monomials {1, x, x 2 , . . .}. However more complex functions such as Fourier bases could be used if for example our output was on a circle. There is some discussion within the emulation community on how much effort should be put into building a good mean, or regression, model and how much the GP can be allowed to fit the large and small scale variation in the data. For example [START_REF] Lim | Design and analysis of computer experiments when the output is highly correlated over the input space[END_REF] make the case for not including a linear model term at all. However most practitioners do include regression terms even if they are only low order polynomials.

To analyse our simulator we use a GP as a prior and combine it with the simulator runs to produce a posterior emulator. Our prior has a mean function as given above

µ(x) = h(x) T β (4)
where the β are parameters that will need to have their own prior specified. We then specify a covariance function

cov(x 1 , x 2 ) = σ 2 c(x 1 , x 2 ; C) (5)
where σ 2 is a variance term and C is a matrix of 'smoothing' terms for the correlation function, c. For example if we have a squared exponential correlation function (sometimes called the Gaussian correlation function for obvious reasons) C might be a diagonal matrix such that

c(x 1 , x 2 ; C) = exp(-(x 1 -x 2 ) T C(x 1 -x 2 )) (6)
C does not have to be diagonal but it reduces the number of parameters to estimate and making this assumption does not appear to impact on our emulators.

Bayes Linear Methods

When we decided to use a Gaussian process as our emulator we made more assumptions than we needed to. The assumption that all the points on our stochastic process have a multivariate Normal distribution is not necessary. All we need assume is that second moments exist and we can then specify the mean and covariance functions without making any assumptions about the statistical distribution of the process.

Such methods that use Bayes theorem but do so in terms of first and second moments are known as Bayes Linear methods and are described in [START_REF] Goldstein | Bayes Linear Statistics, Theory and Methods[END_REF]. The clear advantage over a full Bayes solution is that the priors are also specified only in terms of first and second moments so we do not need to elicit full probability distributions from experts. The second advantage as we shall see is that the calculations are much simpler and faster. However nothing is free and without adding some distributional assumptions we cannot make realisations of the process or produce any form of probabilistic credibility or uncertainty limits. Because Normal distributions are defined by their first two moments there is often confusion between Bayes linear methods and making an assumption of Normality. A true Bayes Linear analysis refuses to make any distributional assumptions; although the results may look similar they are conceptually very different.

As with the GP emulator our basic form is

f (x) = j β j h j (x) + w(x) (7)
The regression terms are identical but the w(x) is not a Gaussian process but rather a general second order process, defined by its covariance function.

The equivalent to Bayes Theorem for Bayes Linear are the Bayes Linear update equations. If θ is a vector of our parameters (in our case this will be the β's, σ 2 and the length scales) and x are the results of our runs

E(θ|x) = E(θ) + Cov(θ, x)V (x) -1 (x - j β j h j (x)) (8)
for the adjusted expectation and by

V (θ|x) = V (θ) -Cov(θ, x)V (x) -1 Cov(x, θ) (9)
The full equations for the posterior moments of the Bayes Linear emulator are rather more complicated and are given in full in the Core BL Emulator Thread of the MUCM toolkit [START_REF] Mucm | MUCM Toolkit[END_REF] 6 MUCM and the Toolkit The toolkit is arranged in thirteen threads plus three large case studies. The threads include two 'core' threads on emulating with Gaussian processes and Bayes linear emulation. Further threads explore how to extend these two core procedures when there are multiple outputs, dynamic emulators, two level emulators (one more complex than the other), combining multiple, independent emulators and using and obtaining derivative information. The issues of model discrepancy, history matching, calibration and sensitivity analysis are also covered as well as the important problems of experimental design and screening. In the near future we expect at least one additional thread on emulators for stochastic simulators. Within what is a relatively short paper it is impossible to cover the breadth of material in the toolkit so we will concentrate on a few highlights.

Building an Emulator

In the MUCM Toolkit we have two threads that cover what is described as the core problem. The core problem is in many ways the simplest application for an emulator. We have a single output of interest, we do have multivariate inputs but no information on the derivative of the output. There are no real world data to compare the simulator to and we are not concerned with making inferences about reality; we are only interested in the simulator. This is not as uncommon as might be thought. Simulators, particularly relatively simple ones, are often used to study an idealised version of the system rather than to make inferences about the real world.

The steps in building an Gaussian process emulator are:

1. Set up the initial GP model. Decide on the form of the regression terms and the correlation function. 2. Decide what priors we are going to use for the GP parameters. 3. Design and run the initial experiment. This is sometimes called the training experiment. Design is an important issue. Computer time is usually limited and each simulator run is expensive so every run must be made to count.

Design has its own thread on the MUCM toolkit and is discussed further below. 4. Build the emulator. 5. Validate the emulator. This is a crucial step. Our first attempts at building good emulators are often failures and we need to establish, and possibly convince others, that we have built an emulator which we can have confidence in. There are two approaches. The first is to use a leave one out cross validation. Each point in the training set is left out in turn and an emulator is built using the remaining points. This reduced emulator is used to predict the left out point. The mean square error of the differences between these predictions and left out runs gives a global measure of the quality of the emulator, while a plot of the individual differences may reveal areas of input space where the emulator is performing badly. The advantage of leave one out is that we do not need to perform any more runs of the simulator.

The disadvantage is that we are not testing the full emulator but emulators built using reduced datasets. Our carefully crafted space filling design is continually compromised by having one point at a time omitted. The alternative is to run a second experiment with an independent set of simulator runs. These are then used to test the emulator. This can be done by simply comparing the difference between the predicted values and the truth (scaled by the predicted uncertainty) similar to the leave out method or there are more rigorous methods based on regression diagnostics that take into account the correlation between points [START_REF] Bastos | Diagnostics for Gaussian Process Emulators[END_REF]. Once we are satisfied that we have an emulator that validates we can combine the training and validation runs and create a single emulator.

The procedure for a Bayes linear emulator is similar but the specification of the priors is in terms of moments rather than full prior distributions.

Design

One of the first steps in any experiment is the design. This is as true of computer experiments as it is of field trials. Our methods have been developed with expensive simulators in mind, so our designs need to minimise the number of runs required. However we also we wish to cover the complete input space. This seems like an impossible task. If we were to use a traditional factorial design with only 2 levels for each input the number of runs required would be 2 p where p is the number of inputs. Two levels of each input gives us very poor coverage and 2 p becomes unaffordable for even a small number of inputs.

We therefore have to look at other designs. There are two main families of design: those based on the Latin Hypercube and those based on Quasi-Monte Carlo sequences.

The traditional design is the Latin Hypercube [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]. In a Latin Hypercube design we first decide how many evaluations of the simulator we can afford in total, let this be n. We then divide the range for each input variable into n equal sections. The simulator is evaluated once, and only once, in each of these settings. This means that we have good marginal coverage of each of the variables. The Latin Hypercube design is now produced by permuting the numbers 1, ..., n for each variable separately. Note the randomisation is to produce a design rather than to randomise for external factors as it is used in field experiments. All Latin Hypercubes have good marginal properties but are not necessarily the space filling. There is no algorithm for an optimal space-filling Latin Hypercube. In general an additional space-filling criteria, maximising the minimum distance between points [START_REF] Johnson | Minimax and maximin distance designs[END_REF] or minimising the sum of inverse distances for example, is imposed on the design. Alternatively we can use Latin Hypercubes based on orthogonal arrays [START_REF] Owen | Orthogonal arrays for computer experiments, integration and visualization[END_REF].

An alternative to the Latin Hypercube is to use a low discrepancy sequence to define the design points. Such sequences were originally devised to efficiently compute multidimensional integrals [START_REF] Niederreiter | Random number generation and quasi-Monte Carlo methods[END_REF]. These sequences are space filling but for small numbers of simulator runs there can be problems with certain projections having 'holes' in them. Examples include the Halton sequences [START_REF] Halton | On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals[END_REF] and Niederreiter nets. For the design of computer experiments the most commonly used low discrepancy sequence is the Sobol sequence [22,?].

Current work in the design of computer experiments explored in the MUCM toolkit [START_REF] Mucm | MUCM Toolkit[END_REF] includes sequential design where we use an initial space filling design to learn about the system and then use this information to guide us on where future runs should be carried out.

An Example

As a simple example consider a energy balance model of the Earth's climate. The Earth is reduced to a line of grid boxes all ocean with a single box below to represent the deep ocean. If it becomes cold enough in a box sea ice forms this has a different albedo to water so a different amount of radiation is reflected from sea ice. As the surface water becomes colder it becomes denser and can sink at a variable location in the North. Deep water upwells back to the surface at a fixed location in the South. Heat is transferred between the surface grid boxes by both advection and diffusion. The system is driven by incoming short wave radiation from the Sun. As in the real world this is greater at the equator than the poles. The total amount of incoming radiation is set by the Solar constant. As this is a simple example we only vary one input, the Solar constant, and we only use a single output, the mean surface temperature. Note as is often the case the output we emulate is a function of the simulator not a state variable.

A space-filling Latin Hypercube in 1-d is simply a set of evenly spaced points. Using 6 points we get the emulator shown in panel (a) in Fig. 2. The dashed line shows the expected value of the emulated and the shaded region gives 95% uncertainty bands. Because we do not have a nugget term in the emulator the uncertainty collapses to zero where we have simulator runs. This makes perfect sense as here we know the value of the simulator output. Note that because this simulator is very cheap to run we can run it across all of input space and plot the true simulator value across the range and this is shown by the solid line in the figure.

We now do a further three runs of the simulator. Fig. 2 shows the standardised residuals for each of these new points. This is the distance between the expected value from the emulator and the true value divided by emulator standard deviation. If the emulator were correct we expect these numbers to have Normal distributions (for simplicity we are ignoring correlation here; for a full solution see [START_REF] Bastos | Diagnostics for Gaussian Process Emulators[END_REF]). One of the points is outside the limits shown in the figure, indicating that we do not have a good emulator. We therefore combine our existing points, both the six original training points and the three new validation points, into a single dataset and build a new emulator. This is shown in Fig. 2 panel (c). We now need to validate this new emulator. A further three runs are carried out and this time all three points are within the bounds for the standard used residuals. The emulator in panel (c) therefore validates. However we now have three additional runs which we can use to build an even better emulator. This, our final emulator, is shown in panel (d) of Fig. 2.

Sensitivity Analysis

One important application of emulators is in sensitivity analysis. Although computer simulators may have large numbers of inputs, often the outputs are dependent on only a few. Formally we can look at sensitivity analysis as answering the question: if the inputs x are changed by a small amount δx what is the effect on the simulator output f (x). One way of looking at this problem is to vary each input in turn, run the simulator and see what effect the change has on the output. This is known as one at a time sensitivity analysis. If we could guarantee that all the inputs were completely independent of each other it might not be a bad idea, but we are dealing with large, complex, non-linear simulators and it would be foolhardy to make such an assumption. Any approach to measuring the sensitivity of a simulator must acknowledge that there will interactions between the inputs and should at least estimate not only the effect of single inputs (main effects) but also at least the first order interactions.

Traditionally this has been tackled by looking at the derivatives of f (.) w.r.to x. This gives the local sensitivity since ∂f (x)/∂x depends on x and may change radically as we move around input space. An alternative is to use variance based sensitivity analyses, see for example [START_REF] Saltelli | Sensitivity Analysis: Gauging the Worth of Scientific Models[END_REF]. [START_REF] Oakley | Probabilistic sensitivity analysis of complex models: a Bayesian approach[END_REF] extend these methods so that they can be used with emulators.

History Matching and Calibration

After sensitivity analysis, probably the most important application of emulators is in comparing simulators with data and calibrating the simulators. We make a distinction between simulator comparison, often referred to as history matching, where we simply rule out areas of input space that are incompatible with the data and calibration where we estimate the 'best' input values from the data. Both rely on the concept of 'model discrepancy'. The simulator is attempting, in most cases, to simulate some property of the real world. But if we measure that property do we expect the simulator to give an exact fit to the data? The answer is almost certainly not. In building simulators we make assumptions, we parameterise processes and there are real world processes we do not include, or even no about. All these omissions and approximations mean that we do not expect our simulator to explain the data perfectly. We describe this difference between the real world and the simulator as the model discrepancy. Model discrepancy will change as the input values change and often will have to be elicited from experts rather than being estimated from the data. Statistically it is possible to think of model discrepancy in two ways: as a bias between the simulator/emulator and the data, or as an increase in variance around the simulator/emulator.

First consider the problem of history matching. Some data are collected and we wish to know which values if any of the inputs could have produced them. Translating this into the negative language of statisticians, who only ever reject hypotheses and never accept them, we are looking for input values that are implausible given the data. This is done via an implausibility measure

I mp = all x (x -x emul ) 2 σ 2 data + σ 2 emul + σ 2 discrep
where x is a data point; x emul is the emulator value corresponding to that data point; σ 2 data is the variance of the data; σ 2 emul is the variance of the emulator; and σ 2 discrep is the model discrepancy expressed as a variance. Values of implausibility greater than 3 (or 5 depending on how conservative we want to be) are ruled to be implausible. Additional simulator runs are then done within the input space not yet deemed implausible (sometimes to as NROY, Not Ruled Out Yet, space). By refining the emulator and adding additional runs in waves within NROY space the volume that hasn't yet been deemed implausible falls rapidly. Table 1, taken from [START_REF] Vernon | Galaxy Formation: a Bayesian Uncertainty Analysis[END_REF], shows the volume of NROY space at each wave of the experiment as well as the number of runs in each wave.

If we want to go beyond history matching to model calibration we need to reformulate the problem. We follow [START_REF] Kennedy | Bayesian calibration of computer models[END_REF]. First we split the inputs into two: control inputs and calibration inputs. The calibration inputs are those inputs we are trying to calibrate. Control inputs on the other hand are inputs that control the simulation but which can't or won't be calibrated. For example in an environmental simulator we may have inputs that give the spatial position of the outputs, these would be control inputs. We now say that the simulator f (x con , x cal ) is the sum of reality y(x con ) and discrepancy d(x con ). Note that in this formulation both reality and the discrepancy are functions only of the control inputs not the calibration inputs. We have some measurements of reality, z(x con ) which are given by z(x con ) = y(x con ) + [START_REF] Kennedy | Bayesian calibration of computer models[END_REF] We can now build an equation that links the data (z) and the simulator.

z(x con ) = f (x con , x cal ) + d(x con ) +

We now use Gaussian processes not only to build an emulator for f (x con , x cal ) but also for d(x con ). The posterior of x cal is the calibrated distribution inputs. For full details see the toolkit [START_REF] Mucm | MUCM Toolkit[END_REF] or [START_REF] Kennedy | Bayesian calibration of computer models[END_REF].

12 Beyond the Core Problem

The core problem as discussed above gives us a basis on which to expand to more complex problems. Examples dealt with in the toolkit include:

1. Multiple outputs where we are interested in more than a single output from the simulator. 2. Dynamic emulators, where we are interested in an output that is itself changing over time. 3. Multiple level emulators, where we have a number of simulators ranging from a coarse, fast simulator through a hierarchy to a slow, but more accurate, simulator. 4. Derivatives. Often we have derivative information available from the simulator in the form of an adjoint model. We can use this information to improve our emulator. However there is a trade off between the extra expense in calculating the derivatives and producing more runs of the simulator without derivatives. In other cases we may not have derivative information from our simulator but we are interested in the form of the derivatives. For instance we may be interested in local sensitivity. Even with the use of automatic differentiating compilers [START_REF] Hascoët | Computing adjoints by automatic differentiation with Tapenade[END_REF] the production of an adjoint for any reasonably sized simulator is a major undertaking. An alternative is to produce an emulator for the derivative. This is relatively easy as the derivative of a Gaussian process is another Gaussian process. Validating the derivative emulator is more difficult if we do not have an adjoint to compare with.

Conclusions

I hope that I have managed to show in this short introduction that both Gaussian process and Bayes Linear emulators are powerful tools in the quantification of uncertainty. The MUCM Toolkit [START_REF] Mucm | MUCM Toolkit[END_REF] is a good reference for these methods in sufficient detail to allow code to be developed. The toolkit is not static and is updated approximately quarterly. In the near future we expect it to be extended to cover stochastic simulators

Fig. 1 .

 1 Fig. 1. The top page for the MUCM toolkit

Fig. 2 .

 2 Fig.2. Panel (a) shows the expected value of the emulator (dashed) and the true simulator based on six runs (o). The shaded region is a 95% uncertainty band around the expectation. Panel (b) shows the values for 3 additional validation runs with the emulator uncertainty. Note that one of the points is outside the uncertainty limits. Panel (c) shows the emulator built using all nine points (the original training set plus the three validation points). This emulator validates with an additional three points (not shown). The final emulator based on 12 points (6+3+3) is shown in panel(d)

  The Managing Uncertainty in Complex Models (MUCM) consortium consists of five UK research institutions (University of Sheffield, University of Durham, University of Aston, London School of Economics and the National Oceanography Centre). As part of its research activities MUCM has set up an on-line toolkit. This consists of more than 300 pages describing most aspects of emulation. The toolkit is not a software package. The best analogy is a recipe book, one of those good recipe books that encourage you to experiment. There are worked examples so you can check that your code works and there are links to existing software packages. Part of an example page is shown in Fig.1. The URL for the toolkit is http://www.mucm.ac.uk/toolkit.

Table 1 .

 1 Table showing the reduction in NotRuled Out Yet space by application of an implausibility measure to the Galform simulator of galaxy formation. For details see[START_REF] Vernon | Galaxy Formation: a Bayesian Uncertainty Analysis[END_REF]

	Wave Runs % Space
	1	993 14.9%
	2	1414	5.9 %
	3	1620	1.6 %
	4	2011 0.26 %

Acknowledgements This paper wouldn't have been possible without the help and support of the MUCM team, in particular Ioannis Andrianakis (NOC), Ian Vernon (U. Durham) and Hugo Maruri-Aguilar (QMUL). John Shepherd (U. Southampton) kindly provided the simple climate model used in section 9.

DISCUSSION

Speaker: Peter Challenor

Maurice Cox : This seems a very useful toolkit. I understand you are on your second grant (MUCM2). I would be interested in hearing plans for your maintenance of the toolkit beyond MUCM2.

Peter Challenor : That is a very pertinent question that is at the top of our agenda. The MUCM2 project runs until October 2012. Maintaining the toolkit frozen from that point would not be expensive and even further developments could be done for little cost and we are looking for funding to put the toolkit on a long term sustainable basis. We would like to encourage non-MUCM participants to contribute to the toolkit and make it a true community resource.

John Reid : You have said nothing about parallel programming. It strikes me that you scope for "embarrassingly parallel" (ideal) execution.

Peter Challenor : Because we rely on ensembles of simulator runs it is quite correct that have an embarrassingly parallel problem. However some of the simulators we are working with are so large that we do not have the computer resources to make use of this. Parallel computing brings up some interesting questions in the design of experiments. I mentioned sequential designs above. Traditionally sequential designs would involve additional single simulator runs, but if we have access to parallel computing it is much more efficient to use 'batch sequential' designs where we run n additional simulations at a time.

John Rice : Have you used your simulator to explore the possible sources of the unusual climate changes observed in historical data? For example, the "little ice age" that occurred a few centuries ago or the major climate change that occurred about 12000 years ago and which appears to have occurred several times earlier with a periodicity of (as I recall) about 125000 years? Historically speaking, we appear to be near the end of a long warm period if this cycle is persistent.

Peter Challenor : The simple Earth radiation balance simulator we use in the example is not suitable for reproducing the history of the Earth's climate. It is much too simple and is lacking too many process. For a good description of the Earth's climate over the last 20000 years see [START_REF] Shakun | A global perspective on Last Glacial Maximum to Holocene climate change[END_REF] Will Welch : The MUCM Toolkit provides recipes for developing code. To test that code it would be useful to have test cases complete with data and results. What plans do you have to provide such test cases?

Peter Challenor : This is a very good point and it is our intention to supply such test cases with most, if not all, the pages. We currently have worked examples for some pages; the example in section 9 is taken from the coreGP pages. We currently have nine test cases plus the three large case studies. Over the next year I hope we will see many more being produced.


Alexandru Cioaca

Adrian Sandu
email: sandu@vt.edu

Eric De Sturler

Emil Constantinescu

Efficient Computation of Observation Impact in 4D-Var Data Assimilation

Keywords: Four dimensional variational data assimilation, observation impact, nonlinear optimization, preconditioning, adjoint model, iterative linear solvers

Data assimilation combines information from an imperfect model, sparse and noisy observations, and error statistics, to produce a best estimate of the state of a physical system. Different observational data points have different contributions to reducing the uncertainty with which the state is estimated. Quantifying the observation impact is important for analyzing the effectiveness of the assimilation system, for data pruning, and for designing future sensor networks. This paper is concerned with quantifying observation impact in the context of four dimensional variational data assimilation. The main computational challenge is posed by the solution of linear systems, where the system matrix is the Hessian of the variational cost function. This work discusses iterative strategies to efficiently solve this system and compute observation impacts.

Introduction

This paper discusses a framework for computing the impact of observations in four dimensional variational (4D-Var) data assimilation. The purpose of this framework is to reveal the degree of usefulness of the observational data in the process of correcting the a priori knowledge of the state of the system. From the analysis of observation impact, certain decisions can be inferred about the strategy to collect observations. The main computational challenge associated with estimating observation impact is the solution of a linear system where the matrix is the Hessian of the 4D-Var cost function. This Hessian matrix is usually very large and accessible only through matrix-vector products. A high accuracy solution is typically sought. In applications such as real-time sensor deployment, this solution needs to be obtained quickly. This paper discusses the particular characteristics of the linear system involved in observation impact calculations. We study the use of several iterative linear solvers that use only Hessian-vector products, and propose several inexpensive preconditioners to improve convergence. Numerical tests using shallow water equations illustrate the improvements in efficiency resulting from this methodology.

The paper is organized as follows. Sect. 2 introduces variational data assimilation. Sect. 3 zooms in on the topic of observation impact in the context of data assimilation. Sect. 4 presents the shallow water equations test problem and its implementation. Sect. 5 gives details of the data assimilation experiment. A discussion of the main computational issues and the proposed solutions can be found in Sect. 6.

Data Assimilation

Data assimilation is the process by which measurements are used to constrain model predictions [START_REF] Daley | Atmospheric data analysis[END_REF][START_REF] Kalnay | Atmospheric modeling, data assimilation and predictability[END_REF]. The information from measurements can be used to obtain initial and boundary conditions that approximate better the real state of the model at a particular moment. Variational data assimilation allows the optimal combination of three sources of information: a priori (background) estimate of the state of the system, knowledge of the interrelationships among the variables of the model, and observations of some of the state variables. The optimal initial state x a (also called the analysis) is obtained by minimizing the cost function

J (x 0 ) = 1 2 (x 0 -x b ) T • B -1 • (x 0 -x b ) + 1 2 N k=1 (H k (x k ) -y k ) T • R -1 k • (H k (x k ) -y k ) (1)
x a = Arg min J

The first term of the sum quantifies the departure from the background state x b at the initial time t 0 . The second term measures the distance to the observations y k , which are taken at N times t k inside the assimilation window. When assimilating observations only at the initial time t 0 , the cost function is known as 3D-Var. Otherwise, when the observations are distributed in time, the cost function is known as 4D-Var. In this study we will focus on 4D-Var. The blockdiagonal background error covariance matrix B is built to take into account the spatial correlations of the variables, as well as the periodic boundary conditions. H k is the observation operator defined at assimilation time t k . It maps the discrete model state x k ≈ x(t k ) to the observation space. R k is the observations error covariance matrix. The weighting matrices B and R need be synthesized in order to have a fully-defined problem and their quality will influence the accuracy of the analysis. The efficient numerical minimization of (1) requires the gradient of the cost function or second-order derivative information when available. 4D-Var usually relies on adjoint sensitivity analysis to provide information about the first and second-order derivatives of the objective function. For instance, the gradient of the cost function can be given by the first-order adjoint model, while the secondorder information can be computed under the form of a Hessian-vector product from the second-order adjoint. Variational data assimilation is an example of PDE-constrained optimization as it is a problem of nonlinear optimization where the the minimization of the cost function is constrained by a numerical model associated with a set of PDEs. The minimizer of the 4D-Var cost function can only be computed iteratively using gradient-based methods, because an analytical solution is almost impossible to derive. When using 4D-Var to correct initial and boundary conditions in a real-time operational setting, we have a limited number of model runs that can be used before the time window expires. This usually implies the minimization process is halted after a number of iterations, so the global minimum might not be attained. Although the most significant decrease in the cost function usually happens during the first iterations, it is important to take into account that the computed solution might not satisfy the optimality conditions. For more details, see the studies in [START_REF] Alekseev | Comparison of Advanced Large-scale Minimization Algorithms for the Solution of Inverse Ill-posed Problems[END_REF][START_REF] Cioaca | Second Order Adjoints for Solving PDEconstrained Optimization Problems[END_REF].

Observation Impact in 4D-Var Data Assimilation

Methods for computing observation impact have been initially developed for the framework of 3D-Var data assimilation in [START_REF] Baker | Diagnostics for evaluating the impact of satellite observations[END_REF][START_REF] Baker | Observation and Background Adjoint Sensitivity in the Adaptive Observation-targeting Problem[END_REF]. The extension of this theory to address 4D-Var data assimilation was only recently accomplished and the details can be found in [START_REF] Daescu | On the Sensitivity Equations of Four-dimensional Variational (4D-Var) Data Assimilation[END_REF][START_REF] Daescu | Adjoint Sensitivity of the Model Forecast to Data Assimilation System Error Covariance Parameters[END_REF].

If the 4D-Var problem is defined and solved accurately, then the analysis x a should provide a better forecast than the background x b . Since one does not usually have access to the real state of the system (otherwise we would not have to perform data assimilation), the accuracy of the analysis is verified through a forecast score against a solution of higher accuracy. Let M be a nonlinear forecast model that discretizes a system of time-dependent partial differential equations, while we will denote its tangent linear and first-order adjoint models with M and M T . A forecast score e(x a ) is defined on M as a short-range forecast error measure:

e(x a ) = (x a f -x v f ) T C(x a f -x v f ) (2)
where x a f = M t0→t f (x a ) is the model forecast at verification time t f , x v f is the verifying analysis at t f and C is a matrix that defines the metric in the state space.

Following the derivations in [START_REF] Daescu | On the Sensitivity Equations of Four-dimensional Variational (4D-Var) Data Assimilation[END_REF], the equation of forecast sensitivity to observations is:

∇ y e(x a ) = R -1 HA∇ x e(x a ) (3)
where

A = [∇ 2 xx J (x a )] -1 (4)
denotes the inverse of the Hessian matrix of the cost at x a and it is assumed to be a positive definite matrix. The same paper [START_REF] Daescu | On the Sensitivity Equations of Four-dimensional Variational (4D-Var) Data Assimilation[END_REF] provides an approximation to the impact of observation, which is defined as the change the observation brings to the value of the 4D-Var cost function. If we consider the innovation vector as δy = y -h(x b ), then a first-order formula would be δe = e(x a ) -e(x b ) = (δy) T R -1 HA∇ x e(x a )

(5)

The elements of the innovation vector are individual observations that get multiplied with their corresponding sensitivity. For those observations that contributed to the decrease of the cost function, the impact (individual value of δe) will be a negative value. The matrix vector product in the matrix A can be computed as a linear system of the form

A -1 • µ 0 = ∇ x e(x a ) (6)
The value of µ 0 is associated with the sensitivity at the initial time. By initializing the tangent linear model with this vector of sensitivities and integrating it from this time to the final time of the assimilation window, we can obtain the sensitivities to observations corresponding to the rest of the assimilation windows.

The problem of computing observation impact has to be solved in real-time since these results are used in applications such as sensor network deployment. Numerical models that simulate the atmosphere are very expensive from a computational standpoint and their adjoint models even more so. Most of the computation bulk is in computing the solution of the linear system A -1 •µ 0 = ∇ x e(x a ). This linear system has to be solved with a high degree of accuracy because its solution µ 0 is a vector representing a "supersensitivity". Any errors that it contains can be propagated and amplified by the numerical models used for forecast or sensitivity analysis. At the same time, there is also a constraint on the number of model runs one can execute before the results are due for operational use. The Hessian matrix of a system associated with a numerical model underlined by PDEs is usually sufficiently complicated that an iterative solver generally will not converge to a satisfactory solution in a few iterations unless it is tweaked accordingly. Accuracy of the solution and fast convergence of the solver are two of the main challenges when computing observation impact in data assimilation, along with means of computing the second-order information.

The Shallow Water Equations Test Case

We consider the two dimensional Saint-Venant PDE system that approximates fluid flow inside a shallow basin (also known as "Shallow Water Equations" (SWE)):

∂ ∂t h + ∂ ∂x (uh) + ∂ ∂y (vh) = 0 ∂ ∂t (uh) + ∂ ∂x (u 2 h + 1 2 gh 2 ) + ∂ ∂y (uvh) = 0 ∂ ∂t (vh) + ∂ ∂x (uvh) + ∂ ∂y (v 2 h + 1 2 gh 2 ) = 0 .
The spatial domain is square shaped (Ω = [-3, 3] 2 ), and the integration window is t 0 = 0 ≤ t ≤ t F = 0.1. Here h(t, x, y) denotes the fluid layer thickness, and u(t, x, y), v(t, x, y) are the components of the velocity field. g is the standard value of the gravitational acceleration. The boundary conditions are periodic in both directions. For ease of presentation, we arrange the n discretized state variables in a column vector

x =   ĥ ûh vh   ∈ R n . (7)
In order to solve these equations, we implemented a numerical model that makes use of a finite volume-type scheme for space discretization and a fourthorder Runge-Kutta scheme for timestepping. This method was introduced by Liska and Wendroff in [START_REF] Liska | Composite Schemes for Conservation Laws[END_REF]. In the following, forward solves will be denoted by FWD.

The tangent linear model (TLM), first-order (FOA) and second-order (SOA) adjoints were generated through automatic differentiation of the forward model using TAMC [START_REF] Giering | Recipes for Adjoint Code Construction[END_REF][START_REF] Giering | Tangent Linear and Adjoint Model Compiler[END_REF]. TAMC is a source-to-source translator that generates Fortran code for computing the derivatives of a numerical model. For more details about automatic differentiation see [START_REF] Griewank | On Automatic Differentiation[END_REF]. The FOA and SOA are integrated backwards in time and they require the values computed at each timestep for the FWD and TLM, respectively. There are two ways to generate these numerical values. The first one is to run the FWD and TLM in advance and to store their evolution in time; this approach has the disadvantage of taking up a lot of storage space. The second one is to recompute the FWD and TLM at each timestep where the FOA and SOA need their values. Theoretically this is necessary at each timestep, which means the latter approach saves storage space at the expense of increasing demand for computation.

A question that arises naturally is how large is the overhead of the code generated through automatic differentiation. Table 1 illustrates the CPU times of TLM, FOA and SOA models, normalized with respect to the CPU time of a single forward model run. It is seen that a full SOA integration is about 3.5 times more expensive than a single first-order adjoint run, while this FOA takes 3.7 times longer than the forward run. The time for SOA can be considered to be a large value as it takes longer time to compute the SOA than to approximate it through finite differences (two FOA runs ). However, for the purpose of this test they are computationally feasible. Please note that these values apply only for our particular implementation (method of lines) and the overhead can greatly vary with the complexity of the numerical method or the automatic differentiation tool used.

We mentioned that each time the tangent linear model or first-order adjoint are run, the forward model is run beforehand for storing its evolution. Similarly, in order to compute the second-order adjoint one must first run the forward, first-order adjoint and tangent linear models. As such, the actual timing values are the ones displayed in the right hand side of the table.

Data Assimilation Scenario

Prior to computing observation impact, we have to perform the 4D-Var data assimilation experiment with the SWE model presented in the previous section. Following are the parameters of the Data Assimilation System (DAS) we are using:

-The 2-D grid of the SWE model is divided into 40 grid points on each direction (1600 in total). The number of time steps for the model runs is fixed at N = 1000. -The reference solution is generated component-wise: the height component h will be a Gaussian pulse of amplitude A = 30 while the velocity fields u and v are assigned a value of 2 at each grid point. -The background solution x b is generated by applying a correlated perturbation on the reference profile for h, u and v. -The background error covariance B was generated for a standard deviation of 5% with a nondiagonal structure and correlation distance of five grid points. This will help the 4D-Var method to spread the information from observations in each cell based on information passed from its neighbors. -The model was run with the reference as initial condition in order to generate the synthetic observations y k . The observation frequency is set to once every 20 time steps. In order to simulate the effect of noise over observations, we apply a normal random perturbation to the perfect synthetic observations. -The observation error covariance R is a diagonal matrix, based on the assumption that the observational errors are uncorrelated. The standard deviation of these errors was set to 1% of the largest absolute value of the observations for each variable.

The observation operator H is configured to select observations for each variable at each point of the grid, hence we say the observations are dense in space. In an operational setting, the operator H is enforced by the availability of observations provided by meteorological stations and other types of sensors.

In order to minimize the 4D-Var cost function, the L-BFGS-B iterative solver will be employed [START_REF] Zhu | L-BFGS-B: A Limited Memory FOR-TRAN Code for Solving Bound Constrained Optimization Problems[END_REF]. The optimization parameters of interest were assigned the following values:

-The stopping criterion is set to

∇J (x 0 [k] ) < 10 -6 • max(1, x 0 [k]
). -Wolfe conditions parameters are set to c 1 = 10 -4 and c 2 = 0.9.

-A maximum number of 500 iterations is allowed for each optimization.

6 Computing Observation Impact

Second-Order Information

The matrix of coefficients is the Hessian of the 4D-Var cost function, evaluated at the analysis. For large-scale models like the atmosphere, this matrix cannot be manipulated in its full explicit form due to high computational demands, storage restrictions and I/O bottlenecks. The dimensions of such Hessian matrices, associated with models of many variables, could reach orders of magnitude of 10 6 or more. In practice, one usually tries to evaluate directly the action of the Hessian matrix on a vector, an approach known as "matrix-free".

Evaluating the Hessian or its action on a vector is itself another computational issue. The methodology of adjoint models provides ways for the computation of first-order and second-order derivative information. Second-order adjoint models are considered to be the best approach to computing Hessian-vector products. Unfortunately, they are not that popular because of the tedious process of deriving them. At the same time, the adjoint models usually have higher computational demands than the forward model, unless they are built around a clever strategy to reuse computation from the forward model [START_REF] Cioaca | Second Order Adjoints for Solving PDEconstrained Optimization Problems[END_REF]. When one does not have access to the second-order adjoint but only to the first-order adjoint, Hessian-vector products can be computed by using finite differences. This approach might be more economical in terms of computation than using secondorder adjoints, but is also less reliable.

These two arguments indicated to us that there is only a limited class of methods we can use in order to solve, precondition or transform the system. More specifically, we will use those methods that need only matrix-vector products. By elementary analysis, the Hessian of any reasonable system is symmetric. However, this symmetry can be lost due to numerical inaccuracies. At the same time, we cannot guarantee the Hessian is positive definite when the 4D-Var firstorder optimality condition is not met and the resulting analysis is far from the true minimizer. While using iterative methods on our linear systems seems to be feasible, it appears to be significantly more difficult to obtain preconditioners, even simple ones based on the diagonal. Also, we can exclude factorizations or direct solve methods since these usually require the full matrix or at least a significant portion of it.

In Figs. 1234we present some graphical representations of the explicit form of the Hessian, built through the second-order adjoint methodology from matrixvector products with the e i unity vectors. For the scenario presented in Sect. 5 the size of the Hessian is 4800 rows by 4800 columns. The first 1600 rows correspond to the variable h, the next 1600 rows to u and the last 1600 to v. The Hessian thus generated proved to be symmetric. Regarding the minimization of the 4D-Var cost function, we let L-BFGS run for as many iterations as needed to obtain a solution close to the optimum (around 500 iterations). As a result, the Hessian was indeed positive definite. Although not shown here, we can confirm that early termination of L-BFGS after 50 or 100 iterations produced an analysis whose corresponding 4D-Var Hessian was no longer positive definite. As said before, the mere purpose of generating the full Hessian is to analyze it offline. In an operational setting it would not be possible, due to constraints imposed by the computational resources.

The first thing noticed about the Hessian is its somewhat regular structure that looks near block-Toeplitz. Although the diagonals are not constant, their values vary inside a well-defined interval. At the same time, we could notice that these diagonals occur at fixed distance from the main diagonal, spaced at every 40 or 1600 rows/columns. These patterns can be seen in Figs. 3 and4, which are graphical representations of the upper-left corner of the matrix. The patterns repeat themselves across the whole matrix. Since the 4D-Var Hessian is the inverse of the covariance matrix of the analysis [START_REF] Gejadze | On Error Covariances in Variational Data Assimilation[END_REF], we can attribute the patterns in the Hessian to the covariance of each grid point with the grid points in its discretization scheme stencil. Thus, each point of the grid is correlated to the points to the East and the West (corresponding to adjacent rows/columns in the Hessian) but also to the points found to the North and the South (corresponding to rows/columns situated at a distance of 40) and with the other variables at the same points (corresponding to a distance of 1600 rows/columns). This can be verified in Fig. 3 where the diagonal breaks down after 40 cells and is displaced, due to periodic boundary conditions.

The eigenvalues of the Hessian are displayed in Fig. 5, sorted in ascending order. Based on their values, the condition number for the matrix is of 7.2095e+3. The magnitude of the condition number does not indicate the matrix is extremely ill-conditioned, but the eigenvalues are not clustered together and many of them are close to zero.

Solving the Linear System

Next, we will present our results and comments on solving the linear system with Krylov-based iterative methods when no preconditioner is used to accelerate the convergence. For our experiments we used MATLAB and its built-in linear solvers. We chose three algorithms for our comparison of their efficiency at solving the system. These three methods are "generalized minimum residual" (GMRES), "conjugate gradients" (CG) and "quasi-minimum residual" (QMR). There are many other variations of these methods but these three are each representative of one class of such methods. While GMRES and QMR operate on general matrices, CG needs the matrix to be SPD which might not be true all the time for the 4D-Var Hessian.

As we can see from Fig. 6 and Table 2, all three solvers converge to a solution whose norm of the residual is close to zero. GMRES and QMR converge to the same solution and they take the same iterates. This means they are equivalent when applied to this problem. CG converges to a more accurate solution from the perspective of its error norm, but the residual norm is larger. This can cause confusion in practice because we would not have access to the reference solution. Also, the norm of the residual for CG fluctuates, instead of decreasing monotonically. Based on these results, we decided to use GMRES for subsequent tests. All solvers start with a vector of zeroes as initial solution.

Preconditioning

The convergence speed of Krylov-based linear solvers depends on the spectrum of the system matrix. These algorithms perform better when the eigenvalues are clustered together, but as seen in Fig. 5 the eigenvalues of our Hessian matrix are scattered across various orders of magnitude. Also, many of these eigenvalues are close to zero which makes them more difficult to be estimated by the solvers. Such a matrix is called "ill-conditioned" and this explains why neither algorithm presented in the comparison of the previous section converged in 100 iterations to a solution that is highly accurate.

In order to alleviate slow convergence, a linear system can be preconditioned by multiplying it with a matrix M (called preconditioner) that approximates the inverse of the system matrix or some of its features, but is easier to invert or factorize. One way to use this matrix is to multiply the linear system with it to the left, which leads to an equivalent linear system. This system should be easier to solve if the preconditioner approximates the inverse of the system matrix, since their multiplication will yield a new system matrix that is better conditioned. In the ideal setting when the preconditioner is exactly the inverse of the original system matrix, the new system matrix becomes the identity matrix and the linear solver will converge in one iteration. However, finding the exact inverse of a matrix is usually a more difficult problem then solving the linear system.

The problem of finding a preconditioner for our system is the fact that we do not have access to the full system matrix. This means we have to exclude the possibility of using certain preconditioning techniques such as incomplete factorizations, wavelet-based, or variations of Schur complement. Moreover, we cannot just use matrix-vector products to construct basic preconditioners such as the main diagonal or a diagonal band, unless we construct the full matrix. Knowing or anticipating the structure of the matrix can be of great help when trying to devise a preconditioner. For instance, we know our Hessian matrix is the inverse of the covariance matrix of the 4D-Var analysis. If the data assimilation process is accurate enough, one can observe a near block-diagonal structure corresponding to each variable. On a 3D grid these blocks correspond to the vertical levels of each variable. This was noted by Zupanski in [START_REF] Zupanski | A Preconditioning Algorithm for Large-scale Minimization Problems[END_REF] who hinted at the possibility to approximate the diagonal in a block-wise fashion, where the diagonal values for a certain block are equal to a constant value. We will try to approximate the diagonal by using Hessian-vector products to "probe" the matrix. A straightforward way of accomplishing this for our three-variable model is to run three Hessian-vector products with unity vectors that extract one column (row) of the Hessian at one time and then use the value corresponding to the diagonal element for all diagonal elements in that block. For example, consider three unity vectors for our 4800 Hessian that have the value 1 at position 1, 1601 and 3201 respectively, and zeros everywhere else. The Hessian-vector product realized with these three vectors will extract the columns number 1, 1601 and 3201 which correspond to the three different variables in our Hessian structure. Building an approximation of the diagonal means using the value found at coordinates 1,1 for the entire first block (so up to coordinates 1600,1600), the value found at coordinates 1601,1601 for the entire second block and so forth. This approximation can be refined by probing for more elements from the same block. If there are many blocks that have to be probed and the computational burden increases significantly, one can employ coloring techniques to probe for more than one element with the same matrix-vector product. We shall refer to this technique as "block-diagonal preconditioner".

Another characteristic of covariance matrices is their diagonally dominant structures and under certain favorable cases, strongly diagonally dominant. In this former case, an approximation of the diagonal can be computed with just one matrix-vector product, where the seed vector is composed only of ones. This is equivalent to computing the sum of each row, which is a good approximation of the diagonal if the diagonal element supersedes the sum of all other elements from the same row. This is the second preconditioner we tried in our experiments.

The Hessian matrix can also be approximated from data collected throughout the minimization process of 4D-Var. Iterative solvers such as L-BFGS work by generating approximations of the Hessian evaluated at the analysis in a certain subspace and minimize the cost function upon it. This approximation is very efficient as it preserves the properties of the matrix (positive definiteness and symmetry) and because it is designed for limited storage. We built a preconditioner by reusing the approximation of the Hessian generated over the last 10 iterations of L-BFGS. Our tests showed this to be as accurate as if using more iterations, a result that confirms the theory. The norm of the error against the reference solution and that of the residual are shown in Table 3. We also present the results obtained with the exact diagonal, although this will not be available in practice, as stated before. The preconditioner obtained from the sum of elements on each row did not improve the convergence at all. This is due to the fact that in our case, the Hessian was not strongly diagonally dominant so this preconditioner was far from approximating the diagonal. The L-BFGS preconditioner brought a slight improvement in the solution and has the advantage that no extra computation is required. The block-diagonal preconditioner improved significantly the accuracy of the solution and behaved as good as the exact diagonal.

Conclusions

This paper presents an adjoint-based framework to compute observation impact in 4D-Var data assimilation. The observation impact calculations need to be performed accurately and rapidly in real-time sensor network deployment problems. The main computational task is solving a linear system whose matrix is the Hessian of the 4D-Var cost function, evaluated at the analysis state. This matrix is typically very large, as each row corresponds to one model state, and it is only accessible through matrix-vector products.

The main contributions of this work are to outline the characteristics of the linear system, to investigate iterative linear solvers, and to propose three efficient preconditioners. Two of the preconditioning methods approximate the diagonal of the Hessian from matrix-vector products, while the third method uses data generated during the 4D-Var minimization to provide a quasi-Newton approximation of the Hessian. Our study shows that these inexpensive preconditioners accelerate convergence. Future work will include the development of improved preconditioners for the linear system, and will study the use of additional iterative methods such as multigrid.

Fig. 1 .Fig. 2 .

 12 Fig. 1. The values found on the first row of the Hessian matrix

Fig. 3 .Fig. 4 .

 34 Fig. 3. Density plot of elements in the (100 × 100 minor) of the Hessian matrix

Fig. 5 .Fig. 6 .

 56 Fig. 5. The computed eigenvalues of the Hessian matrix (in ascending order)

Table 1 .

 1 Ratio between the CPU time of adjoint models and the forward model

	FWD	1
	TLM	2.5 FWD + TLM	3.5
	FOA	3.7 FWD + FOA	4.7
	SOA 12.8 FWD + TLM + FOA + SOA 20.0

Table 2 .

 2 Norm of residual and error for solutions computed with different iterative methods (100 iterations)

	Solver	GMRES CG	QMR
	Norm of residual 0.222	0.893 0.222
	Norm of error	0.357	0.267 0.357

Table 3 .

 3 Norm of residual and error for solutions computed with different preconditioners (100 iterations)

	Preconditioner	None Diagonal Block-diagonal Row sum L-BFGS
	Norm of residual 0.222 0.009	0.009	62.346	0.222
	Norm of error	0.357 0.005	0.005	5.301	0.339

Acknowledgements. This work was supported by National Science Foundation through the awards NSF DMS-0915047, NSF CCF-0635194, NSF CCF-0916493 and NSF OCI-0904397.

DISCUSSION

Speaker: Adrian Sandu

Van Snyder : There are two possibilities to assimilate remote sensing data into chemical transport models (CTMs). One is to incorporate the measured quantity, e.g., radiances. The other is to incorporate geophysical quantities that result from separate analyses. The latter is a nonlinear parameter estimation problem, which depends upon having a good starting point. Further, the problem is ill-posed, which requires a priori information or Tikhonov for stabilization. If the CTM initializes and stabilizes the remote sensing problem, whose results are then assimilated into the CTM, is there a positive, negative, or neutral effect on the solution quality?

Adrian Sandu : This is a very good question. In principle the direct assimilation of measured quantities is to be preferred, since it does not introduce additional errors/biases in the data. In order to do this, one needs good models of the instrument, e.g., that produce radiances from the concentration fields computed by the CTM.

In practice, there are very many instruments that produce data, and for which good models are not available. For this reason assimilation of geophysical quantities, derived from the raw data through an off line estimation process, is often employed. To the best of my knowledge no comprehensive tests have been carried out to date in order to quantify the impact of assimilating geophysical quantities in lieu of measured quantities.