-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathlec8.Rmd
270 lines (214 loc) · 12.2 KB
/
lec8.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
# Spatial temporal variation modelled by differential equations
### Deterministic spatial dynamic models
Deterministic models are based on the assumption that the processes
governing change are known. Given that knowledge, we need
* the state of a system (initial conditions), and
* the characteristics at boundaries of the system (boundary
conditions): what are the sources and sinks, when does what escape or enter the
modelled area.
for a (perfect) prediction of the changes over time, and in space.
Let us look at an example: air quality (fine particles, PM10).
### Model domain
For a model, we need a *model domain*, which is the spatial
area and temporal period over which we will define processes and
conditions. This could be e.g. Western Europe, 1980-2010, or NRW,
2000-2005, or the area not further than 100 m away from the crossing
Weseler Strasse-Bonhoeffer Strasse, Jan 8, 2008, 0:00-24:00. It should
be exactly quantifiable/traceable.
### Initial conditions
The initial conditions usually describe the variable of interest (the
concentration field) at the highest possible level of detail. In our
case this is the PM10 concentration field at the start of the modelling
period.
As this is a continuous field, we need some way to describe this and
usually the spatial domain is discretized into model usually
square, rectuangular or triangular model elements.
This discretization should match the level of detail
* at which we know initial conditions and
* at which we want to model features.
As an example: if we want to quantify the effect of individual
car fumes, spatial elements of 10 cm--1 m may work; if we want to
describe the effect of a complete streets something of 10m--100m
seems more appropriate. Smaller elements and time steps mean more
memory and CPU time requirements.
### Initial conditions-2
If we don't know initial conditions exactly, we may put the starting
point of the modelling domain further back in the past, and hope that
the effect of approximation will damp out as we model. (This assumes
we get the boundary conditions and processes right.)
### Boundary conditions
PM10 comes and goes. Sources are (i) emissions inside the model domain
(cars, households, industry), and (ii) material that enters the model
domain by movement of air bodies, but emitted elsewhere. We need these
source terms (points, lines or fields) in space, and over time.
Sinks are mostly air that moves out of the model domain, and wash out
(rain), dry deposition (your grandmother's white sheets turning black),
and ... inhalation. These terms are also needed, quantitatively.
### Processes
Particles move mostly for two or three reasons: by large-scale movement of
air (wind), by medium/small-scale movement of air (turbulence, dispersion)
and by themselves (diffusion; think Brownian motian of a single particle
in a gas).
As an example, take a look at the LOTOS-EUROS model
(http://www.lotos-euros.nl/) model documentation.
As you can read in the _model formulation and domain_, the model
uses external modelling results (interpolation or mechanistic modelling)
to get the atmospheric driving forces (height mixing layer, wind fields),
e.g. from FUB and ECMWF (http://www.ecmwf.int/).
Basically, the model _code_
* reads a lot of initial and boundary data,
* solves the differential equations and
* writes out everything that is of interest, such as the space-time concentration fields.
### Solving differential equations
The partial differential equation solved in LOTOS_EUROS is
the [continuity equation](http://en.wikipedia.org/wiki/Continuity_equation)
$$
\frac{\delta C}{\delta t} +
U\frac{\delta C}{\delta x} +
V\frac{\delta C}{\delta y} +
W\frac{\delta C}{\delta z} $$
$$
= \frac{\delta}{\delta x}(K_h \frac{\delta C}{\delta x}) +
\frac{\delta}{\delta y}(K_h \frac{\delta C}{\delta y}) +
\frac{\delta}{\delta z}(K_z \frac{\delta C}{\delta z}) +E+R+Q-D-W
$$
with (mostly quoted from the [reference guide](https://lotos-euros.tno.nl/media/10360/reference_guide_v2-0_r10898.pdf), page 6/7:
* $t,x,y,z$ time and the three space dimensions
* $C$ concentration of the pollutant (dynamic field)
* $U,V,W$ the large scale wind components in respectively west-east
direction, in south-north direction and in vertical direction
* $K_h,K_z$ the horizontal and vertical turbulent diffusion coefficients
* $E$ the entrainment or detrainment due to variations in layer
height.
* $R$ the amount of material produced or destroyed as a result of
chemistry.
* $Q$ is the contribution by emissions,
* $D$ and $W$ loss terms due to processes of dry and wet deposition respectively
To solve this equation, it needs to be discretized in space and
time. For this particular model, the spatial grid size is 0.5 degree
longitude $\times$ 0.25 degree latitude (meaning that grid cells
do not have constant area), and time step is 1h.
### Solving PDE's
The simples method to solve PDE's by discretization is, _finite difference_, uses a regular mesh size, $\Delta x$. The solution is computed at
location $j \Delta x$, with $j=1,...,N$.
In one dimension the first derivative uses one of the three approximations:
* backward:
$$\frac{\delta u}{\delta x}(j \Delta x) \approx \frac{u_j-u_{j-1}}{\Delta x} $$
* forward:
$$\frac{\delta u}{\delta x}(j \Delta x) \approx \frac{u_{j+1}-u_{j}}{\Delta x} $$
* centered:
$$\frac{\delta u}{\delta x}(j \Delta x) \approx \frac{u_{j+1}-u_{j-1}}{2\Delta x} $$
and for the second order derivative (centered):
$$\frac{\delta^2 u}{\delta x^2}(j \Delta x) \approx \frac{u_{j+1}-2u_j+u_{j-1}}{(\Delta x)^2} $$
### Diffusion equations, 1-D
Diffusion happens in space-time. Using a mesh in space-time, we can write
$u(j \Delta x, n \Delta t) \approx u^n_j$ with $n$ a superscript, not power.
We can approximate the PDE
$$\frac{\delta u}{\delta t} = \frac{\delta^2 u}{\delta x^2}, \ \mbox{with} \ u(x,0)=\phi(x) $$
by finite difference, in time:
$$\frac{\delta u}{\delta t}(j \Delta x, n \Delta t) \approx \frac{u_{j}^{n+1}-u_{j}^n}{\Delta t}$$
and space:
$$\frac{\delta u}{\delta x}(j \Delta x, n \Delta t) \approx \frac{u_{j+1}^{n}-u_j^n}{\Delta x}$$
Using forward difference for $t$ and centered for $x$, combining both, the
corresponding finite difference equation that approximates the PDE it is:
$$\frac{u_j^{n+1}-u_j^n}{\Delta t} = \frac{u^n_{j+1}-2u_j^n+u_{j-1}^n}{(\Delta x)^2}.$$
### Forward/backward, explicit/implicit
Solving
$$\frac{u_j^{n+1}-u_j^n}{\Delta t} = \frac{u^n_{j+1}-2u_j^n+u_{j-1}^n}{(\Delta x)^2}.$$
is trivial, as $n+1$ is only in the LHS. This means that for each $x$ we can solve
the equation explicitly, where we start is not important. They require, for stable solutions,
that $\frac{\Delta t}{(\Delta x)^2} \le \frac{1}{2}$. See examples.
If the equation were instead
$$\frac{u_j^{n+1}-u_j^n}{\Delta t} = \frac{u^{n+1}_{j+1}-2u_j^{n+1}+u_{j-1}^{n+1}}{(\Delta x)^2}.$$
then we have the unknown $u^{n+1}$ both left and right of the equal sign. This requires the solution
of a (sparse) set of coupled linear equations, and this solution is called _implicit_. It pays off: the
solutions are stable, and larger time steps can be chosen (provided of course, that change is close
to linear over the time step).
### Calibrating deterministic models
Models based on partial differential equations have parameters; think
of diffusion parameters, source and sink terms (boundary conditions),
and initial conditions. These need to be filled in, somehow.
Given that observation data on the model outcome are available, one way
to fill these in is to search for values such that the model predictions
best fit the data. We have seen methods for this; there is a long list
of further, possibly more advanced or efficient methods for finding
optimal parameter values, both in a deterministic (_optimum_) sense,
and in a stochastic (_distribution_) sense.
Also, choosing optimality criterium (least squares? least absolute
differences? combined criteria over multiple variables?)
### Difficulties in calibration
Problems that may occur with calibrating models are numerous.
One problem may be that the parameter we tune (optimize) is not constant
over space or time, but varies. This means that there instead of one
single value, there may be numerous. Their number may outnumber the
observations, and in that case there is little hope in finding realistic
values.
Another problem is that we may tune a parameter and get a better fit, but
that in reality we turned the wrong button, meaning we get a better
fit _for the wrong reason_. This may have disasterous effects when
using this *optimized* model in a prediction setting (such as future
forecasting, or scenario evaluation).
Automatic codes exist (e.g. PEST, or ''optim'' in R) that optimize
models, irrespective what the model is or does.
### More difficulties
Deterministic models use a temporal and spatial discretization. This
is a balance between CPU and memory costs, and the ability to fill the
discrete elements sensibly. Processes need to be *lumped*, meaning
that they cannot be taken into account because of the grid cell size
(think of convection above a forest, or a thunder storm, when grid cell
size is 50 km, and/or time step a day). Choosing a finer resolution,
the parameters, processes, boundary and initial values need to be filled
in with much more resolution (precision), and need disaggregation --
e.g. a country total emission may need to be assigned to 1 km $\times$
1 km grid cells.
### Dynamic parameter updating schemes
A probabilistic setting of a deterministic model is that of the
_Kalman filter_. This algorithm assumes measurements are a combination of
a true state and a measurement noise. Each component has its particular
error structure, expressed by a mean and covariance matrix.
For each new time step, the model predicts a new state, the observations
are compared to that new state, and the model errors are used to adjust
(improve) the model before predicting the next step.
Kalman filters are used a lot to optimize deterministic models, and come
nowadays in many flavours, e.g. depending on whether the model is linear
or not, and whether it is used forward in time, or in real-time.
_Particle filters_ do this in a stochastic setting.
### Simplified difference equation-type models
Often, the differential equations are simplified very crudely to the
state where only mass is preserved, but the solution no longer
approximates the true differential equation. Think of simple bucket-type
models in hydrology, where water bodies are buckets and soils are like sponges:
a soil grid cell drains always with an exponential decrease; a soil and
water body grid cells drain *instantly* when their maximum capacity
is exceeded, with the amount it is exceeded.
Despite the simplifications, these models can be more useful than
attempts to solve the full differential equation, because their data
demand can be more realistically met.
### Calibration vs. validation of models
_Validation_ of models involves assessing their validity, or value,
and is a principally different activity from _calibration_. Where
_calibration_ is an optimization, _validation_ is an assessment. The
validation may be quantitative ("the mean square prediction error
of the model is 0.75") or qualitative ("the model is fit for
its purpose").
_Validation_ may also involve model comparison: the comparison of
different models having different model _structure_, and can in
that sense be seen as an abstraction over calibration:
* calibration compares the models that have identical structure, but have different parameter values
* validation compares models that differ in structure
Model comparison, in th course, would e.g. be:
* comparison of AR(1), AR(2), ..., AR(n) to describe a particular time
series
* comparison of ordinary kriging, universal kriging with predictor
$X_1$, and universal kriging with predictor $X_1$ and $X_2$ for a
given interpolation problem.
In this case, the models compared are nested (i.e., AR(1) is
a special case of AR(2), etc), which may simplify the problem
to testing the significance of the extension (is the additional
parameter significantly different from zero?) If they are not nested,
comparison becomes harder, conceptually.
A famous model comparison on the other end of the extreme
is that of the _coupled model intercomparison project_, CMIP,
e.g. [CMIP](https://cmip.llnl.gov/); it involves large, international
groups of climate scientists who agree, prior to the experiment,
which models to compare.