CC BY
35
Вычислительные технологии
Том 8, № 1, 2003
THE DIFFUSION MODEL FOR LONG-TERM COASTAL PROFILE EVOLUTION: NUMERICAL
VALIDATING*
A. V. Avdeev, E. V. Goryunov Institute of Computational Mathematics and Mathematical Geophysics
of SB RAS, Novosibirsk, Russia M. M. Lavrentiev (Jr.) Sobolev Institute of Mathematics of SB RAS, Novosibirsk, Russia
R. Spigler
Dipartimento di Matematica, Universita di "Roma Tre", Italy e-mail: avdeev@omzg.sscc.ru, mmlavr@nsu.ru, spigler@dmsa.unipd.it
В работе предложен численный подход к одновременному определению двух коэффициентов в обратной задаче для уравнений диффузного типа. Исследуемая модель была предложена в работе [1] для описания процессов долговременной эволюции берегового профиля. Итерационная процедура восстановления двух коэффициентов основана на минимизации функционалов невязки. Численные эксперименты с реальными данными показали эффективность работы предложенных алгоритмов.
Introduction
In the paper [1], a model of diffusion type was proposed to describe the evolution of the coastal profile. The diffusion coefficient in the governing equation (its physical dimension is length squared divided by time) corresponds to the time scale of shoreline change following a disturbance (wave action).
The aforementioned diffusion model can be described by the following basic equation:
d(SX) 32(5X) f 3(5X)\
Here 6X(z,t) represents the change of the cross-shore position (i.e., change in the depth at the distance z from the shore line) of the coastal profile and D(z) is the diffusion coefficient.
Following [1], we give some explanations for special cases of the term g(t, z, 6X, 0(5X)/dz). If g = S(z,t) (an external source function), it is possible to introduce the effects of random
*This research was supported in part by the Russian Foundation for Basic Research under grants 01-0564704, 00-07-90343, and 00-15-99092; the GNFM of the Italian INdAM; UNESCO under contracts UVO-ROSTE 875.629.9 and 875.704.0; and CASPUR-Rome.
© A. V. Avdeev, E.V. Goryunov, M. M. Lavrentiev (Jr.), R. Spigler, 2003.
forcing, along-shore transport gradients, and human interference such as nourishment and sand mining [1].
The linear choice g = B(z) d(SX)/dz or g = B(z)SX is also interesting in view of applications. In these models, the coefficient B(z) represents the speed of along-shore sand wave movement. We assume that g = B(z)SX in (1). So, we shall consider the following.
Inverse Problem. Given the function SX0(t) (the change of the cross-shore position at the point z = 0), find the coefficients D(z) and B(z) such that the solution SX(z,t) to the problem
d (SX ) d2(SX ) + _
= D (z) dz2 + B(z) SX, (2)
SX |t=0 = 0, (3)
d (SX ) = <po(t), SX |z=H = 0 (4)
dz z=0
satisfies the equation (surface measurements)
SX U = SXo(t). (5)
Note that the parameter H can be thought of as an estimate of the depth of closure, i.e., the location where the diffusive and transport phenomena virtually end.
1. Numerical reconstruction of two coefficients of the equation
1.1. Description of the algorithm
As the first step to solve the aforementioned Inverse Problem numerically, we apply the Fourier transform. Thus, the original dynamic problem is replaced by a Helmholtz equation with a complex-valued coefficient. The point is that recovering the space-dependent coefficients D(z) and B(z) in the frequency domain, we do not need to go back to the time domain. Here we have used the term frequency domain for retaining a formal analogy with inverse problems for the wave equation, in which case the frequency domain concerning Fourier images of solutions has a clear physical meaning (see, e.g., [2, 3]).
For numerical processing we apply to the problem (2)-(4) the formal Fourier transform under assumption that the coefficients are smooth: D(z),B(z) G C2(0,H), and u G [u^u2]. So, we consider the problem
d2V + B (z )- iu V = 0, (6)
dz2 D2 (z) ' W
dV
dz
F (u), V |z=H = 0, (7)
z=0
where V(z,u) = e-i^SXdt and F(u) = e-^Vo(t) dt.
The inverse problem we are interested in consists in reconstructing both functions D(z) and B(z) from the additional information
Vo(u) = V(0,u), ui < u < u2, (8)
where [u^u2] is the interval of available frequencies.
We solve the inverse problem (6)-(8) by minimizing the cost functional
iw 2 $[D,B]=/ \Vo(u) - K[D,B](u)\2 du + 3 sup |D - Dest| + 7 sup | B - Best (9)
J z z
where the operator K [D, B](u) maps the current "test" values of D(z) and B (z) into the trace of the solution of the boundary value problem (6), (7) at z = 0. Here 3 and 7 are some weighted regularization parameters, 0 < 3 < 1, 0 < 7 < 1, and Dest and Best are the estimated values of D and B, respectively, that can be obtained from physical measurements.
In the modification of the numerical algorithm which is used here, we operate with the gradient of the "uniform" cost functional (9), i.e., the functional with 3 = 0 and 7 = 0. In this case, after rather technical calculations (omitted here since, e.g., a similar approach can be found in [6]), we obtain the following formulas for the gradients of the cost functional with respect to D and B:
(VD$[D,B])(z) = — 4Re £ B^^F(u)[Vo(u) - F(u) G(z, 0; u)] x
x G(z,z; u) G(z, 0; u) du, (10)
(Vb$[D,B])(z) = 2Re / D-2(z)F(u) [Vo(u) - F(u) G(z, 0; u^ x
J w 1
x (G(z,z; u) (G(z, 0; u) du. (11)
Here G(z, Z; u) is the Green function of the problem (6), (7) and the bar denotes the complex conjugation.
Here we do not consider a rather nontrivial theoretical question of uniqueness as well as existence of the global minimum point of the cost functional (9). Instead, we refer the reader to the paper [4] in which similar questions are discussed for a more simple statement of the inverse problem.
To minimize the cost functional (9), the conjugate direction method was used (see [5]). 1.2. Numerical experiments for synthetic data
Computer codes were prepared in C++; and we used Mathematica 3.0 for verification and visualization. Only resources of a personal computer are needed. On each iteration, the function V(z,u) was computed by the so-called semi-analytical method described in [4].
We approximated the coefficients D(z) and B(z) by piecewise constant functions, i.e., 6 layers with equal width 200 meters. According to the results of auxiliary numerical tests, the regularization parameters 3 and 7 were chosen equal to 0.3 and 0.2, respectively. The results of reconstruction (obtained after 98 iterations) are shown in Fig. 1. In order to test the robustness of the algorithm, the
inversion data were artificially corrupted by adding a randomly distributed white noise. A normally distributed random noise with average fluctuations equal to 5% of the data amplitudes was added to the
inversion data.
The result of simultaneous identification of two coefficients is shown in Fig. 2. The initial approximations for D(z) and B(z) were the same as in the previous test, see Fig. 1.
To conclude this section we would like to stress that the proposed version of the numerical inversion algorithm has demonstrated quite reasonable performance and accuracy. Thus, a basis for real data processing has been established.
3.5
2.5
B
0 0.2 0.4 0.6 0.8 1 1.2
0 0.2 0.4 0.6 0.8 1 1.2
4
3
z
Fig. 1. Simultaneous reconstruction of both coefficients D(z) and B(z). True functions are shown by solid lines, the initial guess by dotted lines, and the result of reconstruction by dashed lines.
3.5
2.5
B
0 0.2 0.4 0.6 0.8 1 1.2
0 0.2 0.4 0.6 0.8 1 1.2
4
3
z
Fig. 2. Results of simultaneous reconstruction of two coefficients D(z) and B(z) from corrupted data.
2. Validating the diffusion model: numerical experiments for real data
Real data representing the long-term evolution of the cross-shore position of coastal profiles were collected over a period of 10 years from 1981 to 1991 at Duck, North Carolina, and are presented in Fig. 3. Real data SXr(z,t) (the subscript "r" stands for "real") for the cross-shore position consist of a 100 x 250 array. This means 100 observation points in the spatial variable z, each being an average over 80 meters, and 250 observation times, each being an average over 15 days. In addition, the source function f (t) was measured at the same observation times. This function represents the average height of waves over the observation periods and is also shown in Fig. 3.
It was supposed that the coastal profile evolution is described by the model equations (2)-(5), where, in equation (2), the function f (t) is added to the right-hand side.
First, the data were numerically recalculated in terms of the Fourier transform F in order to obtain the data Vr(z,u) = F[5Xr(z,t)] and t^(u) = F[f (t)]. More precisely, to solve the problem we use the Fourier-type representation
V(z,u)= i 5X(z, t) e-iWi dt, (12)
J 0
where T is the whole period of observation; with the inverse transform given by the formula
SX(z,t) = V(z,un) e-ni, (13)
n
where un = nn/T, n = 1, 2,... , N.
SXr
r 2 1 0 -1
z, km 0
2
1.75 1.5 1.25 1
0.75 0.5
50
100
150
200
250
Fig. 3. Real data for the cross-shore position and the source term. Then the cost functional was computed by the formula
C^N ¡-H
$
|V - Vsynth|2 dzdu,
(14)
1 J 0
where the data VSynth ("synthetic") were obtained as the solution to the direct problem
d2V
dz2
iu — B T r "0(u) V + w
D2
D2 '
dv
dz
F (u), V |z=H = h(u).
(15)
z=0
The boundary data were obtained as follows:
FM = , hM = Vr(H,w) (Az = 80 m).
0
First, the cost functional was studied in the case where both coefficients of the equation are constants: D(z) = const and B(z) = const. Then a contour plot of the cost functional on the plane (B, D) was calculated (Fig. 4). Each point in Fig. 4 corresponds to the value of the functional (14) for a pair of constants (B,D). The following step-sizes for numerical computations were chosen: AB = 10, AD = 0.02. The intervals for B and D were taken (0, 800) and (0, 0.5), respectively.
D F
400
B
800 B
D 0.4
0
0
Fig. 4. Contour plot and 3D visualization of the cost functional. The point shows the position of the global minimum.
We can see that the cost functional has a rather complicated structure even in the simplified case of the governing equation (2) with constant coefficients D and B.
Usual gradient methods provide very poor convergence in such cases and it is practically impossible to find the global minimum.
Therefore the global minimum of the functional (14) was found numerically by exhaustive search over a rough mesh. We denote the corresponding values of the parameters of the problem by Bmin and Dmin; Bmin = 400 and Dmin = 0.25. Next, the solution Vm\_n(z,u) to the direct problem (15), (16) (with constants Bmin and Dmin substituted for D and B, respectively) was computed.
Finally, the time-dependent profile, $Xmin1 (z, t) was determined after inverse Fourier transformation,
A 3D plot of the profile $Xmin1 (z,t) is shown in Fig. 5. For comparison with the measured data, we calculated the time integral of the relative error:
The average relative error is 15.8%. We can see that the rough assumption that the parameters of the problem are just constants, B = Bmin and D = Dmin, turns out to be in a rather good qualitative agreement with the measured data, see Fig. 3. It is important that in our case the relative error 5(z,t) = |(#Xr — $Xmin1)/$Xr | does not increase with the growth of time.
Then, using the procedure of minimization of the cost functional by the conjugate direction method we reconstructed the coefficients B(z) and D(z) under the assumption that they are piecewise constant functions (see Fig. 6). The range of the spatial variable was divided into 9 layers of equal width. The values D = Dmin and B = Bmin, which were obtained in the previous numerical test, are shown by dashed line.
The corresponding 3D plot of the profile $Xmin2(z,t) is shown in Fig. 7. Comparison with the measured data was made as above, see (17). The average relative error is 5.4%. We can
see (13).
(17)
z, km 0
Fig. 5. Results of numerical reconstruction of the cross-shore position for constant values of D and B.
800 -
600
400
200 -
Fig. 6. Results of reconstruction of D(z) and B(z) for 9 layers of equal width.
0
2
4
6
8
0
2
4
6
8
fiXmin2
2 1 0 -1
z, km 0
Fig. 7. Results of numerical reconstruction of the cross-shore position for piecewise constant D(z) and B (z) (Fig. 6).
see that piecewise constant coefficients provide a better agreement with the measured data, see Figs. 3, 5, and 7.
So, the results obtained can validate the diffusion model which was proposed by de Vriend and Capobianco in [1] for description of the long-term coastal profile evolution.
Acknowledgement
We are grateful to William Birkemeier and his coworkes at the US Army Field Research Facility in Duck, North Carolina, for supplying the high-quality data used in this research. We also thank Michele Capobianco (Tecnomare, Venice) for his interest and a number of useful suggestions.
References
[1] De Vriend H.J., Capobianco M., Chesher T. ET AL. Approaches to long-term modelling of coastal morphology: a review // Coastal Eng. 1993. Vol. 21. (1/3 [special issue]). P. 225-269.
[2] Aki K., Richards P.G. Quantitative Seismology: Theory and Methods. San Francisco: Freeman, 1980.
[3] Avdeev A.V., Lavrentiev M.M. (Jr.), Priimenko V.I. Inverse Problems and Some Applications. Novosibirsk: ICM&MG, 1999.
[4] Alekseev A.S., Avdeev A.V., Fatianov A.G., Cheverda V.A. Wave processes in vertically inhomogeneous media: a new strategy for velocity inversion // Inverse Prob. 1993. Vol. 9, No. 3. P. 367-390.
[5] KARMANOV V.G. Mathematical Programming. M.: Nauka, 1986 (in Russian).
[6] McGillivray P.R., OldenbourGH D.W. Methods for calculating Frechet derivatives and sensitivities for the nonlinear inverse probelms: a comparative study // Geophys. Prosp. 1990. Vol. 38, No. 5. P. 499-524.
Received for publication July 28, 2002
