An Approach to Uncertainty of a Long Range Air Pollutant Transport Model

Working Paper

A N A P P B O A C E T O U N C E B T ~ OF A LONG W G E AIEt POUUTILNT TEANSPOBT YODEL

Joseph Alcamo J e n y Bartnickf

December 1985 UP-85-88

International Institute for Applied Systems Analysis

A-2361 Laxenburg, Austria

AN APPROACH TO UNCERTAINTY OF A LONG RANGE AIR POLLUTANT TRANSPORT MODEL

Joseph Alcamo J e r z y Bartnicki

December 1985 UP-85-88

Worktnq h p e r s are interim r e p o r t s on work of t h e I n t e r n a t i o n a l Institute f o r Applied Systems Analysis a n d h a v e r e c e i v e d only lim- ited review. Views o r opinions e x p r e s s e d h e r e i n d o n o t neces- s a r i l y r e p r e s e n t t h o s e of t h e Institute or of i t s National Member Organizations.

INTERNATIONAL INSTITUTE FOR APPLIED SYSTEMS ANALYSIS 2361 Laxenburg, Austria

PREFACE

One of t h e goals of IIASA's Acid Rain p r o j e c t is to create a model t h a t could b e used in negotiations about control s t r a t e g i e s f o r acid deposition between European countries. To t h a t end i t is necessary t h a t t h e model builders p r e s e n t t h e model u s e r s a c l e a r p i c t u r e of t h e credibility of t h e model. One way to maximize credibility would b e to create a v e r y complex model with as many as possible (mostly non-linear) relationships. O u r stra- tegy h a s been a n o t h e r one: c o n s t r u c t a simple model and evaluate i t s uncer- tainties. Thus uncertainty analysis forms a n important p a r t of t h e Acid Rain p r o j e c t ' s r e s e a r c h agenda. This p a p e r d e s c r i b e s a g e n e r a l framework f o r o u r uncertainty analysis. Moreover t h e a u t h o r s have applied t h e framework to t h e atmospheric submodel of our RAINS (Regional Acidification Informa- tion and Simulation) model. I am convinced t h a t this p a p e r not only i s a sub- s t a n t i a l contribution to evaluation of t h e credibility of RAINS, b u t t h a t i t is also of importance f o r t h e f u r t h e r development of t h e long r a n g e t r a n s p o r t model which i s i n c o r p o r a t e d in RAINS and h a s been built by t h e Norwegian Institute of Meteorology u n d e r t h e Co-operative Programme f o r Monitoring a n d Evaluation of t h e Long-Range Transmission of Air Pollutants in E u r o p e

(EMEP).

This p a p e r i s t h e product of a collaboration with t h e Institute f o r Meteorology and Water Management in Warsaw (Poland) u n d e r a study con- tract "Analysis of Uncertainty in Modeling Atmospheric Processes".

Leen Hordijk Pro Ject L e a d e r

ACKNOWLEDGEMENTS

Wolfgang Schiipp developed t h e computer program f o r t h e Monte Carlo simulation described herein. Maximilian Posch conducted t h e analysis of EMEP source-receptor matrices. The a u t h o r s also wish to acknowledge t h e contributions of Roe1 van A a l s t , Ruwin Berkowicz. Gode Gravenhorst, Syl- vain Joffre, Annikki MgkelB, Gijran Nordlund, S e r g e i Pitovranov and Joop den Tonkelaar to t h e ideas p r e s e n t e d in t h i s paper. The a u t h o r s are grate- ful f o r the g r e a t amount of assistance provided by Anton Eliassen and Jor- gen Saltbones through t h e Meteorological Synthesizing Center-West. Leen Hordijk and Maximilian Posch gratefully reviewed t h e e n t i r e manuscript of t h i s paper. Work described in t h i s p a p e r w a s supported, in p a r t , by a colla- borative agreement between IIASA and the Institute of Meteorology and Water Management in W a r s a w , with Dr. J e n y Pruchnicki as Polish coordina- t o r of t h e agreement. Vicky Hsiung typed t h e manuscript of t h i s p a p e r .

ABSTRACT

This p a p e r presents a preliminary framework f o r analyzing uncertainty of a long range a i r pollutant transport model. This framework w a s used t o assess EMEP m o d e l uncertainty. The uncertainty problem is defined in a decision-making context and a distinction is made between uncertainty analysis, sensitivity analysis, and m o d e l calibration/verification. A taxon- omy is introduced to organize uncertainty sources. The taxonomy includes:

model s t r u c t u r e , parameters, forcing functions, i n i t i a l s t a t e and model operation. These categories are f u r t h e r subdivided into d i a g n o s t i c and fbrecasting components. To limit t h e number of uncertainties f o r quantita- tive evaluation, some uncertainties are "screened". Methods are introduced to evaluate uncertainties. These include (1) Monte Carlo simulation of com- posite parameter, forcing function and initial state uncertainties, and (2) statistical analysis of EMEP source-receptor matrices. Preliminary results of applying this methodology to t h e EMEP model are presented.

-

^vii

-

TABLE OF CONTENTS

P R E F A C E

ACKNOWLEDGEMENTS ABSTRACT

1. INTRODUCTION

1.1 U n c e r t a i n t y a n d M o d e l C r e d i b i l i t y

1.2 S e n s i t i v i t y A n a l y s i s a n d Calibration/Verification 2. PROPOSED FRAMEWORK

3. PROBLEM E'ORMULATION

3.1 T i m e a n d Space Scales 3.2 E M E P M o d e l D e s c r i p t i o n

3.2.1 D e t e r m i n i n g A i r T r a j e c t o r i e s

3.2.2 M o d e l A r e a , E m i s s i o n s a n d M e t e o r o l o g i c a l D a t a

3.2.3 SO2 a n d SO: C o n c e n t r a t i o n

i i i v v i i

1 3 5 9 11 11 1 5 16 1 7

3.2.4 Deposition of Sulfur 4. INVENTORY OF UNCERTAINTY

4.1 Taxonomy

4.2 Application t o EMEP Model

5. SCREENING AND RANKING OF UNCERTAINTY

6. METHODS TO EVALUATE DIAGNOSTIC UNCERTAINTY 6.1 Model Calibration/Verif ication

6.1.1 I n t e r p r e t a t i o n of Model Calibration/Verification 6.1.2 Data Observations and Uncertainty Estimates 6.2 Model S t r u c t u r e

6.2.1 The Linearity Question 6.2.2 Non-Linear Coefficients

6.3 Monte Carlo Analysis of Composite Uncertainty ( P a r a m e t e r Forcing Function Initial S t a t e )

1 J

6.3.1 Overview

6.3.2 Frequency Distribution of Forcing Functions 6.3.3 Frequency Distribution of P a r a m e t e r s 6.3.4 Algorithm f o r Composite Uncertainty 6.3.5 An Example

6.4 Future Work

7. METHODS TO EVALUATE EVRECASTING UNCERTAINTY 7.1 Model S t r u c t u r e : The Linearity Question 7.2 Forcing Functions: Interannual Meteorological

Variability

7.2.1 "Climate Change" Approach 7.2.2 "Past Variability" Approach 7.3 F u t u r e Work

8 . APPLICATION OF UNCERTAINTY TO DECISION MAKING 8.1 P a r a l l e l Models

8.2 U n c e r t a i n t y R a n g e s 9. CONCLUSIONS

REFERENCES

AN APPROACH TO UNCERTAINTY OF A LONG RANGE AKFl POLLUTANT TRANSPORT MODEL

Joseph Alcamo and Jerzy Bartnicki

1. INTRODUCTION

Along with t h e recognition of regional and interregional a i r quality problems, came t h e need f o r new tools to analyze t h e s e problems. Among t h e s e new

tools

are atmospheric long range t r a n s p o r t models which help to

establish t h e relationship between pollutant emissions and t h e i r deposition hundreds o r thousands of kilometers away. The importance given to t h e s e models by t h e scientific community i s c l e a r f r o m r e c e n t national and inter- national publications (see, e-g., OECD (1979), U.S. National Research Coun- cil (1983), MOI (1982)).

A key issue in using t h e s e and o t h e r a.k pollution models for decision making (any mathematical model, f o r t h a t matter) i s t h e credibility of t h e model's results. An essential a s p e c t of this credibility i s how w e l l t h e model u s e r understands t h e model's uncertainty. This p a p e r p r e s e n t s a framework

to comprehensively treat t h e uncertainty of long r a n g e t r a n s p o r t of a i r pol- lutants models (sometimes r e f e r r e d to as LRTAP models) and applies t h i s framework to t h e analysis of uncertainty of t h e so-called EMEP model*

(Eliassen and Saltbones, 1983). From a l a r g e r perspective, w e believe t h a t t h e framework p r e s e n t e d h e r e i n c a n b e generally applied to o t h e r t y p e s of environmental models. Throughout t h e p a p e r w e (1) discuss key issues con- c e r n e d with uncertainty analysis, (2) p r e s e n t numerical examples of dif- f e r e n t aspects of t h i s analysis based on preliminary r e s u l t s from t h e IIASA Acid Rain P r o j e c t , (3) denote f u t u r e work t h a t will b e conducted within t h e frame of t h e IIASA P r o j e c t . Since t h i s r e s e a r c h is only in i t s e a r l y s t a g e s , w e intend t h i s to be a discussion paper.

In t h i s p a p e r w e are specifically i n t e r e s t e d in determining the uncer- t a i n t y of w i n g model results i n a decision-making contezt. Our goals f o r t h e uncertainty analysis include:

1. To quantify, where possible, t h e combined uncertainties of many d i f f e r e n t uncertainty s o u r c e s , i.e. determine t h e uncertainty r a n g e of model calculations.

2. To determine under what conditions t h e model performs best.

3. To make more explicit t h e assumptions behind model p a r a m e t e r s , forcing functions, etc.

4. To identify t h e s o u r c e s and r e l a t i v e importance of uncertainties as a guide to model use and s e t t i n g r e s e a r c h priorities.

he

EhIEP model is described i n S e c t i o n 3.2 of t h i s paper.

The analysis r e p o r t e d in t h i s p a p e r builds on previous work on model uncertainty in t h e fields of decision analysis (cf. Howard and Matheson, 1983) and econometrics (cf. Griliches and Intrilligator, 1983), as _{w e l l} as investigations in water quality modeling (cf. Fedra, 1983 a n d Beck and V a n S t r a t e n , 1983) and ecological modeling (cf. G a r d n e r et a l , 1982). Compared t o t h e s e fields, much less quantitative analysis h a s been conducted on atmospheric model uncertainty. A notable exception i s t h e work done at Carnegie-Mellon University ( Morgan et al. 1984). Also, a r e p o r t from an American Meteorological Society Workshop outlines some key issues in t h e quantitative assessment of atmospheric models (Fox, 1984). Unfortunately a review of t h e aforementioned work i s outside of t h i s p a p e r ' s scope.

1.1- Uncertainty and Model Credibility

Model credibility i s based on s e v e r a l ill-defined c r i t e r i a . One c r i t e r i o n is t h e s c i e n t ~ c b&s of t h e model equations, i.e. t h e soundness of t h e physical/chemicaVbiological c o n c e p t s behind t h e model. Another i s v e t i j t - c a t i o n and v a l i d a t i o n , generally meaning t h e comparison of model r e s u l t s with observations and t h e examination of model behavior to see if i t i s real- istic. Still a n o t h e r way to e n h a n c e t h e credibility of model r e s u l t s i s to p e r - form sensitivity analysis. Collectively, t h e s e a p p r o a c h e s make model u s e r s more confident in using a model y e t they d o not specifically a d d r e s s t h e question of t h e u n c e r t a i n t y of model results. In t h i s s e n s e model u n c e r - t a i n t y i s t h e d e p a r t u r e of model calculations from c u r r e n t or f u t u r e "true values". Mathematically, o u r meaning of uncertainty c a n b e e x p r e s s e d as t h e following.

Let us assume t h a t a n environmental model can b e expressed

*

as:

Y

= G(X)

where

Y

=

(y , y n ) is a n output vector (model results)

X

=

(zl,

...,

z,) in an input vector (input model variables)

6

is an o p e r a t o r (usually a differential).

Since t h e input vector usually contains variables which a r e dependent on space and t i m e , t h e output vector is also a function of space and time. In addition, output variables depend on s o m e constants in time and space, i.e.

parameters.

If w e assume now, t h a t "true values" of t h e output variables are represented by vector Y , t h e model uncertainty can be defined as:

C = Y - r

where:

I t should be mentioned here. t h a t it is extremely difficult to compute t h e complete uncertainty vector because, among o t h e r reasons, "true values" are illusive. There are ways however to circumvent this problem.

Repeated comparisons of observed versus m o d e l computations (model

.Even though t h e model d e f i n i t i o n used i n t h i s paper is n o t thd t h e most general possible, it is s t i l l general enough f o r most o f t h e environmental models.

calibration/verification) yields insight t o r , though in sections 1.2 and 6.1 w e discuss drawbacks to this approach. Another strategy, which is dis- cussed in Section 6.3, is to assess t h e uncertainty of t h e

X

vector in (1.1), and compute a new Y. This provides a n indirect estimate of c. Other stra- tegies are reviewed in t h e text.

The equation (1.2) used h e r e to define uncertainty is related as w e l l to model calibration/verificatioh. However, t h e important difference is, t h a t in t h e case of model calibration/verification, t h e components of t h e vector

T

have to be measurable, while in case of uncertainty, this i s not necessary. In this sense o u r definition of uncertainty i s more general.

1.2. Sensitivity Analpis and CalibrationNerification

Though "sensitivity analysis" and "model calibration/verifi~tion" are relevant to a model's uncertainty, both approaches have limitations.

S e n s i t i v i t y a n a l y s i s in t h e conventional sense is difficult

*

to perform f o r two or more variables and tends to emphasize extreme events. It is dif- f e r e n t from model uncertainty because sensitivity analysis is interested in t h e incremental changes of model results caused by an incremental change in input variables. In f a c t , t h e objective of m o s t s e n s i t i v i t y analyses Is, of course, to determine t h e relative importance of one independent variable compared to another; not how much model calculations d e p a r t from reality.

In this sense sensitivity analysis is an essential p a r t of model development.

Mathematically w e can e x p r e s s sensitivity analysis as a procedure f o r com-

laather than add y e t another definition of sensitivity analysis we quote a published do- finition: "Sensitivity analysis involves

...

maklng a series of runs with a model and noting the magnitude of the changes in results as assumptions, parameters and initial conditions are changed in an orderly fashion." (McLeod, 1982. p. 96).

puting matrix S:

The elements of t h e sensitivity matrix S

=

[si,] a r e given by t h e relation:

Model calibration/veriftcatioh, i.e. comparison of model output with observations h a s t h e following limitations in assessing model uncertainty:

(i) Observations a r e o n e n u n r e l i a b l e . Eliassen and Saltbones (1983) p r e s e n t one example of analytical e r r o r s in sulfate d a t a used to check EMEP model calculations.

(ii) Model o u t p u t is n o t n e c e s s a r i l y "observable" in n a t u r e , espe- tidy iJ the model describes an aggregated s y s t e m . For models with l a r g e temporal/spatial resolution such as t h e EMEP model, it is difficult t o r e l y on comparisons of model output with observa- tions. S t r i c t l y speaking, since t h e EMEP model computes SO2 gas and SO; in r a i n o v e r 150 km long orthogonal coordinates and a 1 km v e r t i c a l mixing layer, model output should be checked with observations averaged o v e r t h e same spatial scale. This class of e r r o r i s termed a g g r e g a t i o n e r r o r and has been dealt with in s o m e detail in t h e ecological modeling l i t e r a t u r e (Gardner, et ^d, 1982). A r e l a t e d problem o c c u r s when an important model output i s virtually unobservable as in t h e case of total sulfur deposition.

(iii) C e r t a i n cause-eflect r e l a t i o n s m a y not be r e a d i l y observable. An example of this is the relationship between sulfur emissions from a p a r t i c u l a r country and i t s deposition at a particular location in

Europe. Though wind sector malysis may help to quantify this relation for s h o r t periods of time, i t is difficult to d o so o v e r a longer time scale, say one y e a r . Nevertheless this time scale and relationship i s computed by t h e EMEP model and is of p a r t i c u l a r importance in decision making.

(iv) Agreement @ model o u t p u t with d a t a does not settle t h e ques- t i o n of model u n c e r t a i n t y w h e n t h e model i s u s e d jbr forecast- i n g purposes. For example, model agreement with observations does not address t h e impact of interannual meteorological varia- bility on t h e uncertainty of model forecasts.

(v) Sometimes model parameters c a n be " a r t ~ c i a l l y t u n e d " s u c h t h a t model o u t p u t closely agrees with data. Under t h e s e cir- cumstances i t may a p p e a r t h a t t h e model has little or no uncer- tainty, although t h e uncertainty has simply been t r a n s f e r r e d to t h e uncertainty in choosing t h e correct parameters for forecast- ing purposes.

(vi) It is often d m c u l t to assemble d a t a for a comprehensive r a n g e of environmental conditions. Even though w e test t h e model against d a t a f r o m s e v e r a l time periods, w e still may have l o w con- fidence t h a t w e have covered a representative r a n g e of environ- mental conditions.

Despite t h e preceding caveats, model calibration/verification remains t h e only s u r e "benchmark" of a model's relationship to reality. A s such,

model c a l i b r a t i o n / v e ~ c a t i o n together with sensitivity a n a l y s i s is

n e c e s s a r y a n d us*l t h o u g h i n s y r j t d e n t in evaluating environmental

model uncertainty. In the following sections we propose a comprehensive framework to assess model uncertainty which incorporates elements of both model calibration/verif ication and sensitivity analysis.

2. PROPOSED FRAYLEWORK

A comprehensive approach to analyze long range t r a n s p o r t m o d e l uncertainty should include t h e following:

(i) R o b l e m F o r m u l a t i o n

-

Despite t h e trivial nature of this step it is surprising how often investigators discuss uncertainty of a m o d e l without specifying t h e time and space scales of interest. In Section 7.2 of this p a p e r w e present a n example of t h e dependence of m o d e l uncertainty on t h e temporal-spatial dimensions of t h e prob- l e m . Before proceeding with an uncertainty analysis it is there- f o r e vital to carefully formulate the problem of interest.

(ii) I n v e n t o r y of U n c e r t a i n t y

-

In this s t e p w e assemble and classify the sources of uncertainty for f u r t h e r analysis. Our goal _is to list as comprehensively as possible every major source of uncer- tainty. To do this w e propose a taxonomy of m o d e l uncertainty in Section 4.1 of this paper.

(iii) Screening a n d R a n k i n g

~

U n c e r t a i n t y

-

Virtually every m o d e l used to describe a real system w i l l have a very large number of uncertainties. To l i m i t t h e sources of uncertainty f o r quantitative evaluation w e try in this step to identify t h e most important sources. This is accomplished through conventional sensitivity analysis or qualitative judgement and need not have time-space scales identical to those in step number one.

(iv) E v a l u a t i o n qf U n c e r t a i n t y

-

The sources of uncertainty which remain after s t e p (iii) can be evaluated by a number of different quantitative techniques. Sections 6 and 7 describes some

a p p r o a c h e s being t a k e n in t h e IIASA Acid Rain P r o j e c t ' s analysis of EMEP model uncertainty.

(v) Application to Decision Making

-

^{Once an} estimate of u n c e r t a i n t y i s d e r i v e d in s t e p (iv), we s t i l l must i n t e r p r e t t h i s estimate in a way useful to decision making. F o r example, we could e x p r e s s t h e u n c e r t a i n t y of EMEP calculations of s u l f u r deposition as s p a t i a l v a r i a t i o n s of deposition isolines, or as deposition r a n g e s a r o u n d individual isolines. Alternatively we could apply an "average"

u n c e r t a i n t y estimate to e a c h EMEP g r i d element. These a n d o t h e r a l t e r n a t i v e s are a d d r e s s e d in Section 8 of t h i s p a p e r .

3. PROBLEM FORKULATION

3.1. Time and S p a c e Scales

The d e g r e e t o which uncertainty can vary depending on spatial- temporal scales is illustrated in Figure 7.1 which summarizes a n analysis of uncertainty in computed sulfur deposition due to interannual variation of precipitation and wind patterns. Since we have specified above that w e are interested in "determining t h e uncertainty of using model results in a decision-making context", w e must now clarify t h e t i m e and space scales relevant to decision-making. First. w e assume t h a t w e are interested in a specific source-receptor relationship f o r sulfur emissions. sulfur deposition and a i r concentration. Next w e assume t h a t t h e country-scale is t h e a p p r o p r i a t e spatial-scale f o r sulfur emission sources because (1) most countries in Europe r e p o r t t h e i r sulfur emissions as country totals, although a f e w r e p o r t additional spatial information, (2) most proposed international control policies (for example, t h e 'Thirty Percent Club") refer to country-scale sulfur emissions. The EMEP grid element is a n a p p r o p r i a t e spatial scale f o r r e c e p t o r sulfur deposition since a coarser resolution would b e unsuitable f o r analyzing known spatial variations of environmental impact (such as f o r e s t damage) which occurs within coun- tries. In addition, since a model f o r analyzing international control policies in Europe should cover all of Europe, a spatial scale much smaller than 150 km may increase t h e number of computational s t e p s to an unacceptable level. Moreover, t h e spatial resolution of meteorological data in Europe is

4 2 also approximately 10 km

.

The time scale of t h e source-receptor relationship should take into account that confidence of any a i r pollution m o d e l increases with t h e averaging period of results

* .

In addition, t h e time s t e p should be compatible with t h e long t i m e period and broad spatial coverage needed f o r policy analysis. With t h e s e considerations in mind, a n annual time s t e p i s taken to b e an appropriate scale. This t i m e s t e p i s also appropriate for assessing forest damage since m o s t field studies r e c o r d annual pollutant deposition or air concentration.

W e may summarize t h e discussion to this point by specifying t h e source-receptor time resolution as one year, countty-scale as t h e spatial resolution f o r sulfur emissions, and EMEP grid element as t h e spatial reso- lution f o r sulfur deposition and a i r concentration. The relationship of interest, therefore, between deposition and sulfur emissions can be expressed as:

d i j = s t aij

where

dv ⁼

^total sulfur deposition at grid element j due to country i (g

s

m -2 y r -I)

st

=

total sulfur emissions from country i ( t S y r -I) atj

=

element of source-receptor matrix

W e define our uncertainty r of deposition as

'%

( 3 . l a )

where d i j i s t h e "true" deposition.

- ~ e an example, one EMEP revlew s t a t e s that the model "continued t o demonstrate i t s ef- fectiveness in modelling air concentrations and depositions when averaged over seasons or years" (WMO, 1903).

W e are also interested in t h e uncertainty of t h e total deposition at grid element j (where b j is background deposition):

f o r j=l.. m and

E ~ ,

=

d j

-

^{d j (3.2b)}

The same form of equations (3.1a) through (3.2b) can be applied to t h e o t h e r EMEP state variables (e.g. SO2 air concentration). These o t h e r state vari- ables will be introduced in t h e next section.

It follows from t h e above t h a t w e are interested in t h e uncertainty of computed annual sulfur deposition at various locations In Europe, where these locations are defined by EMEP grid elements. This can be expressed e i t h e r as a n uncertainty range around a linear source-receptor relation- ship (Figure 3.1) o r a frequency distribution of computed sulfur deposition (Figure 3.2). In summary. equations (3.1) and (3.2) define our uncertainty problem. Figures 3.1 and 3.2 illustrate this problem graphically.

SLlLFUR EMISSIONS

COUNTRY l1 I l1

1

^{UNCERTA 1 RANGE}

LINEAR RELATIONSHIP EMISSIONS

-

DEPOSITION

SULFUR DEPOSITION AT LOCATION " J ~

Figure 3.1. Source receptor uncertainty.

FREQUENCY OF OCCURRENCE

SULFUR DEPOSITION AT LOCATION ' I J ~

Figure 3.2. Frequency distribution of computed deposition.

3.2.

W

Model Dercription

A general form of t h e EMEP model is:

where B c

-

B t

=

change in concentration with time.

U ,V

=

orthogonal wind velocities.

B c B c

- - ⁼

orthogonal concentration gradients.

B z ' 6 y

4 ⁼

change in concentration due to chemical reactions and sink processes.

4i

=

pollutant emissions.

The EMEP model uses a Lagrangian approach to solve equation (3.3).

Concentrations of SO2 and SO; are computed along a moving fx-ame of reference. The computation procedure consists of two steps. Trajectories are f i r s t calculated, and then mass-conservation equations are solved on these traJectories to compute t h e concentrations at the r e c e p t o r point. An additional procedure is used f o r computing dry and w e t deposition of sulfur.

The theoretical formulation of t h e EMEP model is described by Eliassen and Saltbones (1975) and Eliassen (1978). This model is similar to t h e one used in t h e OECD program (OECD, 1979). The main difference is that t h e EMEP model is based on trajectories followed for 96 hours instead of 48 hours, and grid size of 150 km instead of 127 km.

3.2.1. Determining Air Trajectories

A trajectory can be considered as the path of an air parcel followed by t h e wind. In the EMEP model two-dimensional trajectories are calculated which neglect vertical motion of the air. The wind field from t h e 850 hPa level is assumed to be t h e transport wind within the mixing layer.

Petterssen's method (Petterssen, 1956) was chosen f o r numerical computa- tions of the trajectories. If

z

is the position of t h e trajectory at t i m e f , t h e next position E

+

&-is calculated using the wind field (Z.f) as follows. Let

&-, be the f i r s t estimate f o r t h e position increment :

&-,

=

C(Z.f)

Af

The i-th estimate &-< f o r & is:

N e w estimates f o r & are computed until:

where r is a small positive number equal 0.003 in t h e EMEP model. If t h e condition (3.6) is satisfied f o r i-th estimate then:

ls

=

Gf (3.7)

The t i m e s t e p

At

is 2 hours, which means that each trajectory is represented by a set of 49 discrete points, including t h e r e c e p t o r point.

This procedure is sufficiently fast and in m o s t cases condition (3.6) is quickly satisfied.

3.2-2. Yodel Aren, M o n a and Meteorological Data

The coverage of t h e EMEP model i s shown on Figure 3.3. I t c o v e r s a l l Europe, a l a r g e p a r t of t h e Atlantic Ocean and a small p a r t of N o r t h e r n Africa. The numerical g r i d system has 39 points in t h e x-direction and 37 in t h e y-direction. As w a s mentioned e a r l i e r , t h e grid size i s 150 km. The grid elements are identified by t h e coordinates ( i s j ) . The relation between g- graphical latitude p, and longitude X and a point ( i , j ) i s given by t h e equa- tions:

p = w - - 360 Arctan r

7r R 7r (3.8)

-(I _d

+

^sin-) ₃

where

r

=

$(i

-

³¹²

+ -

^3712'

( t h e coordinate of t h e Northern Pole is (3,37))

R

=

^6370km

-

r a d i u s of t h e E a r t h d

=

1 5 0 km

-

^{grid size}

A l l meteorological and emissions d a t a are given in t h e g r i d sysiem denoted by equations (3.8) and (3.9). The meteorological d a t a a r e : wind field at 850 h P a level

-

e v e r y 6 h o u r s (with l i n e a r interpolation in-between), and precipitation f o r t h e last 6 hours. In t h e routine computations emission data were taken from an inventory p r e p a r e d by Dovland and Saltbones (1979).

Seasonal variation of emission is introduced into t h e model calculations. I t h a s a s h a p e of sinusoidal function with amplitude 302 and maximum in t h e beginning of January.

Figure 3.3:Area covered by the EMEP model calculations. Tralectories are followed from arrival points within heavy line.

Concentrations of SO2 and SO: and dry and w e t deposition of sulfur are computed f o r t h e e n t i r e grid system.

3.2.3. S& and SO: Concentrations

Emissions are computed by l i n e a r interpolation in each of 49 points from t h e f o u r n e a r e s t g r i d points, and o c c u r r e n c e of precipitation is checked. Having t h i s information, equations f o r SO2 and SOT c a n b e solved. Denoting SO2 concentration by cl and SO: concentration by c 2 (both measured in sulfur units), we can write t h e s e equations in t h e follow- ing form:

The o p e r a t o r

-

D i s t h e total time derivative, Q is sulfur emission p e r d t

unit area and time. Values for all o t h e r symbols and p a r a m e t e r s in equations (3.10) and (3.11) are given in Table 3.1.

Table 3.1. Parameter values in t h e EMEP long-range t r a n s p o r t model (from Eliassen and Saltbones, 1983).

Notation Explanation Parameter value Parameter unit

d Deposition velocity for SO2 8x10 -3 m s - I

vds Deposition velocity for SO: 2 x i 0 " m s Y 1

h Mixing height 1000 m

kt Transformation rate of SO2 to SO: 2 x 1 0 ~ s -I kw W e t deposition r a t e of SO2,

used only in grid elements and six-hour periods when i t rains a Additional dry deposition in

t h e same grid square where emission occurs

Part of sulfur emission assumed to be emitted directly as sulfate Overall decay rate f o r SO:

Proportionality coefficient in equation (3.13)

In Finland and Norway In o t h e r countries

Background concentration in equation (3.13)

In Finland and Norway 0.27 m 9 s l - l

In o t h e r countries 0.40 m 9 s l - l

Equations 3.10 and 3.11 a r e ordinary linear equations solved f o r a particular trajectory. In regions where precipitation o c c u r s e i t h e r con- stantly o r not at all, t h e r e is also t h e analytical solution f o r these equations presented by Eliassen (1978).

3.2.4. Deposition of Sulfur

D r y deposition of sulfur is computed by applying deposition velocities t o SO2 and SO; concentrations:

d d

= GI.

Zld

+

^{C 2} Zldr) T where:

d d =dry deposition of sulfur during time T T

=

period of t h e transport (T= usually 1 y e a r

in t h e EMEP model).

and o t h e r variables are as previously defined.

In t h e routine model w e t deposition is not calculated directly from t h e mass-conservation equations 3.10 and 3.11, because of t h e constant k, rate. It is estimated by a n indirect method instead, in which t h e mean con- centration of sulfur in precipitation E3 is estimated from t h e computed mean concentration of sulfate during t h e rain 62 using a linear empirical rela- tionship:

where C12 and E3 are averaged over time T. The precipitation-weighted mean

c2

is calculated from

where pi is t h e amount of precipitation observed on day i , c2,( is t h e corresponding calculated daily mean a i r concentration of SO4= and P is t h e total amount of precipitation during time T in a specific grid element. Days without precipitation d o not contribute to

Cz.

The empirical proportionality coefficient a in (3.13) corresponds t o a scavenging r a t i o f o r anthropogenic sulfate. The constant b accounts for background concentration in t h e rain. The values of a and b are a l s o given in Table 3.1.

The value of t h e w e t deposition in t h e model d, i s computed as:

d,

= 6, . ^P

and total deposition of sulfur (1, is:

d t

=

dd

+

^{ti,,, (3.15)}

Units of d d ,d, , and d t a r e in g

-

^{m -2.} In o r d e r t o compute t h e m a s s dep*

sited in a grid element, t h e values of t h e deposition must be multiplied by t h e area of t h e grid element.

4.1. Taxonomy

A f t e r formulating o u r uncertainty problem in Section 3 of this paper.

w e now wish to assemble and classify t h e sources of uncertainty. To assist in this classification w e propose t h e following taxonomy of m o d e l uncertainty:

(1) Model Structure (2) Parameters (3) Forcing Functions (4) Initial State (5) M o d e l Opel-ation.

Uncertainty due to Model SYructure results f r o m imperfect or inaccu- l-ate representation of reality by a model. In this sense m o d e l s t r u c t u r e i s taken as t h e collection of model variables and parameters together with t h e i r relationships.

P b r a m e t e r s are defined as those variables which are constant in e i t h e r t i m e o r space, are usually estimated or confirmed as p a r t of t h e m o d e l calibration, and are meant to approximate a more complicated pro- cess.

Forcing j b n d i o n in this taxonomy is a m o d e l variable which inherently changes in time and space, serves as input for m o d e l calibration, and i s assumed to be w e l l known ( o r at least b e t t e r defined) compared to parame- ters.*

*Forcing f i n c t i a corresponds to the concept of input distudhann i n systsms science terminology and crogcnmrs variable i n econometric terminology.

Initial State uncertainty results from t h e e r r o r in assigning boundary and initial conditions.

Finally, uncertainty due to Model m e r a t i o n r e f e r s to e r r o r s in t h e solution techniques of model equations o r in processing model input and out- put. e.g. numerical e r r o r arising from approximation of differential equa- tions and interpolation of model input and output. The sum of forcing func- tion and initial state errors can also b e termed i n p u t uncertainty.

Each of t h e above c a t e g o r i e s can b e f u r t h e r sub-divided into t w o addi- tional classes: diagnostic uncertainty and forecasting uncertainty. Diag- nostic uncertainty p e r t a i n s to model use in simulating p a s t and c u r r e n t conditions. Forecasting uncertainty arises when t h e m o d e l i s used to esti- mate f u t u r e conditions.** Each source of uncertainty (according to t h e model taxonomy, p r e s e n t e d above) h a s both a diagnostic =d forecasting component.

Before proceeding with t h e application of t h e above taxonomy to t h e EMEP model, w e note t h a t t h i s taxonomy i s hierarchically organized as illus- trated in Figure 4.1. This figure notes t h a t uncertajnties due to parameters, forcing functions, initial state and model operation depend on model s t r u c - t u r e . A s a n example. let us assume t h a t w e are uncertain of t h e e x a c t value of t h e dry deposition velocity ud in t h e EMEP model, but can estimate i t s interval as [vd

1.

W e t h e n estimate t h e uncertainty of computed sulfur depo- sition by using, f o r example, a Monte Carlo technique described in Section 6.3. This computed uncertainty will depend on t h e form and content of t h e

==Other inveetigators use different terms to make the same dietinction. For example, Beck (1983) uses Zdmt*cation and R e d i c t i o n .

m o d e l equations, i.e. t h e model structure. Thus it is unlikely t h a t the uncer- tainty of Model 'A' w i l l be exactly t h e same as t h e uncertainty of Model 'B' even if both m o d e l s use duplicate environmental conditions and parameters values, etc., as input. This idea is illustrated in Figure 4.2. In o t h e r words, quantitative estimates of uncertainty due to parameters, etc., pertain only to a particular model.

Figure 4.3 also notes that parameter, etc., uncertainty depends on . ^t

environmental conditions. This i s also obvious if w e consider t h a t the uncer- tainty of vd will have a small influence on computed sulfur deposition if con- ditions are very wet, i.e. if w e t deposition i s t h e predominant sulfur removal mechanism. For d r i e r conditions t h e r e v e r s e will be true. This implies t h a t conclusions about model uncertainty must always include information about t h e environmental conditions under which these uncertainty estimates were made. This leads to the concept of a "frequency distribution of uncertainty"

and "expected value of uncertainty", illustrated in Figure 4.4.

Of course if t h e d e p a r t u r e of model calculations from observations is relatively constant for many different environmental conditions, then w e may suspect t h a t r w i l l also not vary very much f o r different environmental conditions.

4.2. Application to EMEP Model

The diagnostic and forecasting uncertainties due to model s t m c t u r e w i l l be t h e same if t h e system doesn't change, 1.8. if t h e m o d e l contains t h e dominant variables and interrelationships of the real system for both future and c u r r e n t conditions. However, if for example t h e alr concentnitions of co-pollutants such as Og or NOz increase such that they affect t h e transfor-

MODEL STRUCTURE

ENVIRONMENTAL C O N D I T I O N S

PARAMETERS FORCING I N I T I A L MODEL FUNCTIONS STATE OPERAT I O N

Figure 4.1. Hierarchy of model uncertainty

DEPOSITION AT LOCATION

n ~ n '

Figure 4.2. Uncertainty (e) depends on model structure.

r\

I \

RIDEL m

Fa?

FREQUENCY OF

OCCURRENCE MODEL

m

Fa?

E N V 1 m w C T A L COTDI

/

Figure 4.3. Uncertainty (c) depends on environmental conditions.

FREQUENCY OF

OCCURRRJCE

Figure 4.4. R e q w n c y distribution of u n c s r h t n t y

(r)

f o r a mnge of sn- vironmenhl conditions.

mation of SO2 or XI:, t h e n a new model structure may be r e q u i r e d with new v a r i a b l e s and p a r a m e t e r s . This implies t h a t forecasting and diagnostic uncertainty d u e to model s t r u c t u r e will b e different.

F o r t h e EMEP model, uncertainty due to model structure (diagnostic c a s e ) includes; but i s hot l i m i t e d to:

(1) Simplification of air chemistry

-

including, t h e question of linearity.

(2) Assumption of single v e r t i c a l Layer.

(3) Simplification of d r y S deposition p r o c e s s . (4) Simplification of wet S deposition p r o c e s s .

(5) Assumed immediate complete mixing of emissions into mixing layer.

( 6 ) Omission of horizontal diffusion.

( 7 ) Omission of vertical advection and related phenomena (e.g. frontal movements and deep convection).

(8) Omission of exchange between boundary layer and free atmosphere.

(9) Omission of shallow convection.

(10) Omission of orographic effects.

Model parameter uncertainty should also be t h e same f o r both diagnos- t i c and forecasting cases unless the system changes. As a n example of how t h e "system could change", let us assume t h a t t h e SO2 wet deposition

rate,

k,, i s currently oxidant-limited. In this case w e should expect t h e uncer- tainty of k,

to

change in the future if t h e background level of oxidant increases, i.e. such t h a t kk, is no longer oxidant-limited. For the EMEP model, uncertainty w i l l arise f r o m the parameters listed in Table 2.1.

In comparison

to

parameter uncertainty, t h e m is a clear difference between diagnostic and forecasting uncertainty f o r jbrcing jbncfions. In t h e diagnostic case uncertainty arises f r o m o u r i n t e r p r e t a t i o n of t h e actual forcing functions, l.e. either data is incomplete o r w e must transform i t to make i t compatible with t h e model. W e can illustrate this point by con- sidering t h e use of precipitation d a t a as a forcing function of t h e EKEP model. The density of precipitation stations from which these data are derived is very crude compared

to

EMEP's spatial coverage. Consequently, these data must be interpolated before they can s e r v e as _input to the EMEP m o d e l . This "interpretation" is a n example of diagnostic error due

to

t h e m o d e l ' s forcing functions. Other e r r o r of this nature arises from estimation of S emission and wind velooity fields.

On t h e o t h e r hand, t h e forecasting aspect of forcing function uncer- tainty r e f e r s to o u r inability to accurately f o r e c a s t f u t u r e forcing func- tions. In o t h e r words, w e can only estimate t h e magnitude of f u t u r e precipi- tation, S emissions and wind velocity. Interannual meteorologic variability and f u t u r e climate change are p a r t of t h i s c a t e g o r y of uncertainty. The s o u r c e s of forcing function uncertainty in t h e

EMEP

model are summarized in Table 4.1.

T h e r e is also a difference between diagnostic and forecasting uncer- tainties related to i n i t i d state u n c e r t a i n t y . As with t h e forcing functions, uncertainty arises in t h e diagnostic case because we are unable to accu- r a t e l y translate actual boundary and initial conditions into our model. I.%, we cannot input t h e e x a c t SO2 and S3: boundary and initial concentrations into our model. The sources of poesible initial state uncertainty are sum- marized in Table 4.2.

F o r t h e forecasting case, uncertainty arises because w e are unable

to

exactly estimate t h e initial states at t h e beginning of t h e forecasting period.

Uncertainty due to t h e final category of uncertainty in t h e

EMEP

model, model operation includes: ⁽¹⁾ input-output processing, (2) t r a j e c - tory computations, and (3) solution of

EMEP

equations. One type of input- output processing e r r o r is t h e interpolation of input emissions' data. A s an

Table 4.1. Forcing functions in t h e EMEP long-range transport model.

Symbol Q

Explanation Unit

Sulfur emissions k t yr Components of t h e ^m s 1 t m s p o r t wind

vector

Precipitation

Table 4.2. Initial S t a t e Uncertainties in EMEP Model.

Horizontal boundary conditions Vertical boundary conditions Initial conditions

example, emissions' input ^data are illustrated in Figure 4.4. During tsaJec- tory calculations, however, emissions' data are interpolated as in E'igure 4.5 which causes s o m e error in t h e input to equations (3.10) and (3.11)*.

Uncertainty in t h e trajectory calculations arises from t h e so-called h t t e r s s m method described in Section 3.2.

Uncertainty due

to

solution of EMEP equations refers

to

t h e technique for solving equations (3.10) and (3.11).

Finally, w e summarize t h e above uncertainties of t h e EMEP model in Table 4.3.

.1n t h i s example we call t h e transformation of actual sulfur e m i s s i o n s t o t h e input data i n Figure 4.4 jbrcing finction uncertainty, and the interpolation of t h e m input data by t h e model from Figure 4.4 t o Figure 4.5 a s model operation uncertainty. W e t e r m it operation unartointy b.craw it roferr Lo an intrrnd opermtlon of a .p.dnc IPDd.L

1 2

3 4 5

6 7 8

9 GRID POINT

NlMBER

Figure 4.5. Form of emission data which are input ta the EMEP model.

-

SULFm

EMISSIONS -

,

I I

1 2 3

4 5

6 7

3 9

GRID POINT MMBER

Figure 4.6. Interpolation of emission data from Figure 4.5 for model computations.

Table 4.3. A summary of some EMEP model uncertainties.

DIAGNOSTIC UNCERTAINTY

1. MODEL STRUCTURE A s listed in section 4.2 2. PARAMETERS Estimation e r r o r s of

parameters in Table ^3.1

3. FORCING FUNCTIONS Estimate of current magnitude and spatial distribution of sulfur emissions

'Smoothing" e r r o r s and measurement e r r o r s of

of precipitation 'Smoothing" e r r o r s and measurement e r r o r s of wind velocities 3. INITIAL STATE Estimation and

approximation e r r o r s :

-

boundary conditions

-

initial conditions 4. MODEL OPERATION Input/output

processing

FORECASTING UNCERTAINTY

Changes in co-pollutant concentrations

Changes in co-pollutant concentrations

Forecasted sulfur emissions

Interannual meteorological variability (precipitation and wind patterns)

Long term cllmate change

Future boundaryand and initial conditions

Forecasting uncertainties same as Diagnostic

uncertain ties Trajectory calculations

for processing wind velocity data

5. SCREENING AND RANKING OF UNCERTAINTY

The goal of t h i s screening e x e r c i s e is to limit t h e number of uncertain- t i e s which must b e quantitatively evaluated in t h e n e x t s t e p of t h e uncer- tainty analysis. To d o s o w e d r a w on sensitivity analyses conducted by EMEP and o t h e r investigators (cf. Eliassen and Saltbones, 1983 and Anon, 1983).

as w e l l as additional calculations, reviewed in a s e p a r a t e p a p e r (Bartnicki a n d Alcamo, forthcoming).

A s pointed o u t in Section 1.2, t h e r e are difficulties in translating r e s u l t s of sensitivity analysis to conclusions about uncertainty since sensi- tivity analysis and uncertainty analysis (as defined in t h i s p a p e r ) have dif- f e r e n t goals. W e t h e r e f o r e t a k e a pragmatic a p p r o a c h , and r a t h e r than eliminate any uncertainties from f u r t h e r consideration, we assign them

to

categories of f i r s t and second p r i o r i t y . W e will b e conservative and assign

to

t h e second p r i o r i t y only those uncertainties where t h e r e is strong evi- dence t h a t they are less important than o t h e r uncertainties. Remaining uncertainties are considered t o have f i r s t priority. A s will b e seen, most uncertainties are placed in t h e f i r s t priority category, though f u r t h e r sen- sitivity analyses may permit u s in t h e f u t u r e

to

i n c r e a s e t h e number of

"second priority" uncertainties.

Model S t r u c t u r e

-

The following uncertainties described in Section 4.2 are assigned

to

a lower priority: (1) Assumption of immediate complete mix- ing of emissions into mixing layer, (2) Omission of horizontal diffusion, and (3) Omission of shallow convection,

(1) T h e r e are physical r e a s o n s why t h e u s s u m p t i o n of immediate complete m i t i n g of e m i s s i o n s i n f o the m i z i n g l a y e r would not add l a r g e

uncertainty t o t h e EMEP calculations. During t h e day, especially with con- vective conditions, pollutants a r e mixed relatively quickly a f t e r emission.

The c h a r a c t e r i s t i c time in which a p a r c e l of pollutants r i s e s

to

t h e top of mixing l a y e r i s less than one h o u r (Lamb, 1984). which is less than t h e com- putational time step. Therefore, t h e assumption of complete initial mixing of pollutants c a n b e justified f o r daytime t r a n s p o r t . During t h e night, although stable conditions usually inhibit v e r t i c a l mixing, lateral mixing still o c c u r s because of t h e different heights of emissions. Even if t h e above arguments are not s t r o n g enough to put this uncertainty into a lower prior- ity, in p r a c t i c e t h i s phenomenon is parameterized by coefficient a in equa- tion (3.10), t h e local deposition coefficient. Consequently we t a k e this model s t r u c t u r e uncertainty into account by investigating t h e parameter uncer- tainty of a.

(2) The o m i s s i o n of h o r i z o n t a l d m s i o n i s considered less important than o t h e r uncertainties because of t h e smaller scale e f f e c t s of this diffu- sion compared

to

t h e scales t r e a t e d by t h e EMEP model. Considering t h e l a r g e initial size assumed f o r a p a r c e l of a i r pollutant in t h e EMEP model (150 x 150

x

1 km) w e d o not expect horizontal diffusion to affect t h e mixing of pollutants during t h e lifetime of a typical 96 h o u r t r a j e c t o r y . In support of t h i s conclusion, Prahm and Christensen (1977), using a n Eulerian one- l a y e r model similar to t h e EMEP model, found small changes in computed SO2 and SO; air concentrations (around 3X) when they compared models with and without horizontal diffusion.

(3) Shallow convection intensifies pollutant mixing within t h e mixing l a y e r and also chemical transformation of pollutants. The f i r s t effect would influence t h e value of a in equation (3.10) and t h e second effect, kt in t h e

same equation. W e can t a k e this uncertainty into account, t h e r e f o r e , by investigating t h e parameter uncertainty of a and kt.

The remaining model s t r u c t u r e uncertainties listed in section 4.2 a r e placed in t h e f i r s t p r i o r i t y category.

Model Arrameters a n d Forcing ALnctions

-

All of t h e s e uncertainties are c u r r e n t l y considered very important.

I n i t i a l State

-

Of t h e uncertainties of this type listed in Table 4.3, w e may consider t h e uncertainty due

to

t h e vertical b o u n d a r y condition

to

b e contained in t h e t h e uncertainty of parameter b in equation (3.13). We can t h e r e f o r e t r a n s f e r t h i s t y p e of initicrl s t a t e uncertainty to parameter uncertainty. Consequently "vertical boundary condition" as a s e p a r a t e uncertainty h a s been placed in a lower priority category.

Model m e r a t i o n

-

Bartnicki e t al. (forthcoming) p r e s e n t evidence t h a t uncertainty due t o t h e trcy'ectory ccrlcdation method is relatively unim- portant. They examined analytical v e r s u s numerical t r a j e c t o r i e s f o r a n artificial rotational wind and found t h a t a f t e r 96 hours of travel, t r a j e c t o r y positions differed by less than 15 km.

The s o u r c e s of uncertainty c u r r e n t l y assigned 'Second Priority" a r e presented in Table 5.1. By default all o t h e r uncertainties have a higher priority.

-

³⁷

-

Table 5.1. "Second Priority" Model Uncertainties

MODEL STRUCTURE

Assumption of immediate complete mixing of emissions into mixing layer.

Omission of horizontal diffusion.

Omission of shallow convection.

2. MODEL PARAMETERS

-

^none

3. MODEL FORCING FUNCTIONS

-

^none

4. INITIAL STATE

Vertical boundary condition.

5. MODEL OPERATION

Trajectory calculation method.

6. METHODS TO

EVALUATE

DIAGNOSTIC UNCERTAINTY

In this section w e concentrate on t h e evaluation of diagnostic uncer- tainties assigned "first priority" a f t e r screening and ranking of Section 5.

Unfortunately, a discussion of all important diagnostic uncertainties is out- side t h e scope of this paper, though w e t r y and highlight some of t h e more important sources.

6.1. Model C a l i b r a t i o n N e r i f i c a t i o n

In Section 1.2 w e reviewed t h e drawbacks to model calibration/verification in assessing model uncertainty. W e also pointed out t h a t model calibration/verification is necessary and useful though insuffi- cient f o r this task. I t i s necessary, as noted e a r l i e r , because without data comparisons w e have no standard with which to compare model output with reality. I t i s u s m l because t h e d e p a r t u r e of model output from observa- tions i s a measure of t h e magnitude of model diagnostic uncertainty. The goodness of this measure naturally depends upon the amount of data and r a n g e of environmental conditions that t h e model i s tested against. Since diagnostic uncertainty varies with environmental conditions as noted in Sec- tion 4.1, t h e more environmental conditions, i.e. data sets, t h e model i s tested against, t h e closer w e come to t h e "expected value" of diagnostic uncertainty. Using t h e EMEP model as an example, a single comparison of annual average model output with observations provides only one value of model uncertainty under specific environmental conditions. This com- parison does, however, give us an idea of t h e possible maximum diagnostic uncertainty.

6.1.1. Interpretation of Hodel Calibration/Verification

There is no straightforward way to translate t h e results of model cali- bration and verification into estimates of E , model uncertainty. In Figure 6.1 we illustrate t h r e e possible ways in which model calibration and verification can b e combined in modeling practice. In Case I, model parameters are adjusted s o that model output a g r e e s with Data S e t A ("calibration"). Using these calibrated parameters and new forcing functions, model output is com- pared to Data S e t B ("verification"). W e denote t h e d e p a r t u r e of model out- put from observations as clo and Elbl respectively. In Case 11, the model parameters are again adjusted s o that model output a g r e e s with Data S e t A.

This exercise is repeated with new forcing functions s o t h a t model output also a g r e e s with Dah S e t B. W e have t h e r e f o r e "calibrated" t h e model separately to two independent data sets and obtain two independent parame- ter sets. In Case 111, w e are interested in finding t h e single parameter set which fits both Data Sets A and B simultaneously. In o t h e r words, this pr*

cess yields a single parameter set for t h e two data sets.

In general, given an identical model and a n identical calibration pr*

cedure f o r all cases, w e expect

E l a

=

' 2 0 ' 2 b

<

&3a E3b

<

E l b (6.1) Each epsilon is an estimate of diagnostic u n c e r t a i n t y since we assume o u r forcing functions and initial states a r e input to t h e model. Individually they are not necessarily good estimates of average diagnostic uncertainty.

Even though E~~ and E~~ from Case I1 are smaller than E~~ and E3b f r o m Case 111, uncertainty w a s "conserved", i.e. "apparent" diagnostic uncer- tainty has decreased but we have increased forecasting uncertainty

Figure 6.1. Three w a y s of combining model calibration and verification.

L b

C A S E I

DATA

SET

A

t "CALIBRATION" --

^{-P SET}

^1A

MODEL

t "VERIFICATION"

/

DATA

SET

B

C A S E 1 1

DATA

SET

A

t

^N

CALIBRATION"

PARPMETER SET

ZA

MODEL

t "CALIBRATION" -

PAMMETER SET

2B

DATA

SET

B

C A S E I 1 1

DATA

SET

A

t "CALIBRATION

^It

MODEL PARAMETER SET

3

DATA t

SET

B ' t ~ ~ ~ ~ ~ O N " /

3

1

because we d o not know which of t h e two parameter s e t s in Case 11 to use f o r forecasting.

In conclusion, caution must b e used in interpreting comparisons of model output with observations.

6.1.2. Data Observations and Uncertainty Estimates

L a t e r in this section we describe how to obtain quantitative uncertainty estimates of t h e EMEP model. The simplest check of t h e s e uncertainty esti- mates would be, of course simply to compare t h e observations with t h e com- puted frequency distribution. In t h e hypothetical example in Figure 6.2a w e compare observed annual SOz air concentration with t h e frequency distri- bution of computed SO2. W e expect, f o r example, t h a t 90% of t h e time an observation such as this would fall within t h e frequency distribution's 90%

confidence interval. Though a single comparison of this nature proves lit- tle, a comparison of observations versus computed frequency distributions at five stations would serve as a check on o u r procedure f o r analyzing diag- nostic model uncertainty. (This of course also depends on t h e accuracy of t h e data.) In fact, t h e probability t h a t all five observations a r e outside t h e 90% confidence intervals of t h e frequency distributions is ( 0 . 1 ) ~

=

0.001 %.

The probability of two o r more being outside t h e 90% confidence interval i s 1%. This s e r v e s as a way

to

check o u r procedure with data.

6.2. Model Structure

There i s only one satisfactory way t o evaluate t h i s uncertainty and t h a t is of course to compare different model s t r u c t u r e s . If alternative models are available then these comparisons can b e performed with identical input

t

^OBSERVED

9

^{A I R}

CONCENTRATION

Figure 6 . 2 ~ . Hypothetical comparison of observation with computed fre- quency distribution of SO2 air concentration.