Environmental Engineering Reference
In-Depth Information
The main program is given in the following:
function
par_estc
% transport parameter estimation with derivatives
global
xfit cfit T D c0 c1
% Example values for Marmara Sea Sediment Core
T = 6.3e11; % [s] 20.000 years
D = 1.0e-5; % [cm*cm/s]
c0 = 0; %
c1 = 619; % [mmol/l]
xmax = 4000; % [cm]
% specify fitting data
xfit = [0 20 40 60 100 120 140 160 215 255 275 300 375 450 525 600 750
1050 1200 1350 1650 1950 2250 2550 2700 3000 3450 3900];
cfit = [619 597 608 615 619 615 621 571 621 618 619 625 577 612 608 612
609 590 582 582 556 494 457 489 487 444 381 371];
v = fzero(@myfun,0.2e-8,options);
display (['Best fit for v = ' num2str(v)]);
x = [0:xmax/400:xmax];
h = 1./(2.*sqrt(D*T)); e = diag(eye(size(x,2)));
plot (xfit,cfit,'o',x,c0+0.5*c1*(erfc(h*(x-v*T*e'))+...
(exp((v/D)*x)).*erfc(h*(x+v*T*e'))),'-');
legend ('given','modelled');
xlabel ('depth [cm]'); ylabel ('chloride concentration [mmol/l]');
text(0.1*xmax,c1*0.65,['sedimentation velocity [cm/a]: '
num2str(v*3.15e7)]);
e = diag(eye(size(xfit,2)));
normc = norm(cfit-c0+0.5*c1*(erfc(h*(xfit-v*T*e'))+…
(exp((v/D)*xfit)).*erfc(h*(xfit+v*T*e'))));
text(0.1*xmax,c1*0.6,['norm of residuals: ' num2str(normc)]);
The corresponding M-file
'parestc.m'
is included in the accompanying
software.
The value for
T
represents 20,000 years. This is approximately the time when
a fresh water lake, located where the Sea of Marmara is found today, was flooded
from the rising Mediterranean. The value for
D
is a standard first guess for molecular
diffusivity, which is valid for many substances dissolved in water. It is assumed that
before the flooding the sediment pores were filled with water of low chloride
concentration. A more realistic low non-zero value could have been taken, but
this has no significant influence on the model results, as the inflow concentration for
saline water with 619 mmol/l is quite high. The maximum length for the simulation
corresponds to 40 m slightly exceeding the maximum depth of the measurement
locations.
What follows is already the
fzero
-command, which here is the optimization
routine. Starting value for the velocity is 2·10
-9
cm/s. All further commands concern
the design of the graphic in which measured data and the modelled curve can be
compared, and where the result of the estimation procedure is depicted in addition.
The output is shown in Fig.
10.7
.
