- Email: [email protected]
Application of a three-dimensional hydrodynamic model for San Quintin Bay, B.C., Mexico. Validation and calibration using OpenDA
Accepted Manuscript Application of a three-dimensional hydrodynamic model for San Quintin Bay, B.C., Mexico. Validation and calibration using OpenDA Mariangel Garcia, Isabel Ramirez, Martin Verlaan, Jose Castillo PII: DOI: Reference:
S0377-0427(14)00236-2 http://dx.doi.org/10.1016/j.cam.2014.05.003 CAM 9671
To appear in:
Journal of Computational and Applied Mathematics
Received date: 1 October 2013 Revised date: 28 April 2014 Please cite this article as: M. Garcia, I. Ramirez, M. Verlaan, J. Castillo, Application of a three-dimensional hydrodynamic model for San Quintin Bay, B.C., Mexico. Validation and calibration using OpenDA, Journal of Computational and Applied Mathematics (2014), http://dx.doi.org/10.1016/j.cam.2014.05.003 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
mariangelapril2014.tex Click here to view linked References
Journal of Computational and Applied Mathematics 00 (2014) 1–9
Comp. and Applied Math.
4th International Congress on Computational Engineering and Sciences
Application of a Three-Dimensional Hydrodynamic Model for San Quintin Bay, B.C., Mexico. Validation and Calibration Using OpenDA. Mariangel Garciaa,1, Isabel Ramirezb,2, Martin Verlaanc,3, Jose Castilloa,4 a Computational b Centro
Sciences Research Center, San Diego State University, San Diego, C.A, U.S de Investigacion y Educacion Superior de Ensenada, B.C, Mexico c Deltares, Delft, The Netherlands
Abstract A 3D hydrodynamic model (Delf3D) was applied for San Quintin Bay (SQB), B.C., Mexico. Calibration and validation were conducted using measured bathymetry, wind, water surface elevation and velocities fields. The calibration period took place during the winter of 2010. Model predictions were evaluated graphically and statistically against field observations in order to quantify the accuracy of the predictions and evaluate the success of the model calibration. A 2D Calibration of tide constituents at the boundaries, seabed roughness and depths was done to evaluate the improvement in the velocities for the 3D Model. Comparisons between the model predictions and the field observations of the water surface elevation at interior stations indicated that the model had been successfully calibrated. It also indicated that the model predictions were highly correlated with the observed water surface elevations and velocities. This connection between the observed and simulated values was detected through graphical comparisons and root-mean-square errors and correlation coefficient. The objective of this study was to demonstrate that OpenDA could be used effectively to rapidly calibrate a hydrodynamic model for smaller sub-meso scales like the size of SQ Bay. Keywords: Calibration, OpenDA, Data Assimilation, Delft3D, San Quitin Bay, Hydrodynamic.
1. Introduction
improving the starting position of the forcing for a forecast, so that the estimates are different in each cycle. OpenDA software was designed to quickly and efficiently facilitate the use of a tool which enables the reduction of uncertainty of a given model. The generic OpenDA data assimilation environment [1] provides both a range of filtering routines, as well as uncertainty and sensitivity analysis routines. It also supports assimilation of recent observations to improve model forecasts, and allows the user to carry out sensitivity analysis and simultaneous parameter optimization of model parameters in a given model application. It is an updated, open-source version which includes the extended functionalities of the existing DATools system, which has been successfully used in data assimilation of current and salinity profiles [2], for flood forecasting purposes [3] as well as for tidal sensitivity analysis[4]. However OpenDA has been applied before to calibration of coastal seas, but not yet to the submeso scales. A sensitivity analysis of the tidal representation in Singapore Regional Waters was done using 2D WAQUA by [5] but the boundary conditions were not calibrated. In this paper we are using Delft3D for the hydrodynamics model focus in a local scale problem where the open-boundary plays a very important role as well as the dynamics at the entrance.
Numerical ocean models have proven to be essential for the creation of forecast estimates, as well as weather forecasts for coastal areas. However, model errors are currently inevitable, due to uncertainties in both the initial and boundary conditions, as well as in the oversimplified or entirely neglected physical processes and mathematical approximations used in the systems equations. Furthermore, despite numerical weather prediction applications, available coastal ocean observation data is relatively scarce, and often does not capture complex coastal dynamics. Measurement error is another source of variability that is often not well understood, making it difficult to characterize them. Two methods of reducing the uncertainty within the models are presented here. The first method are calibration, with a goal of tuning a set of parameters that are fixed in time, for example: depth, amplitude, phase of tidal constituents and roughness. The second method is data assimilation, with a goal of 1 [email protected] 2 [email protected] 3 [email protected] 4 [email protected]
1
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
The region of study selected was San Quintin Bay (SQB), located South of Baja California., Mexico. SQB forms an interesting ecosystem, which lends itself for the development of aquaculture. In particular, excellent results have been obtained in the growth and survival of the Japanese oyster Crassotrea gigas [6]. One of the most important biological processes for cultivating this organism is its reproductive cycle, which is principally regulated by the amount and quality of its food, as well as by water temperature [7]. The hydrodynamic must be known with accuracy in order to ensure efficient ecological monitoring of the bay. The hydrodynamics of SQB has been studied by [9], [10], and [11]. One of the first results found was by [9] who applied a one-dimensional numerical model to reproduce the hydrodynamics of the lagoon. Observations of tides, currents, temperatures, salinity and winds were used in that research, however the numerical model that was used proved to be of little predictive value for this location. Similar values calculated and observed values of water surface elevations and associated velocities were obtained by applying a two-dimensional hydrodynamic model, which determined that the dynamics of SQB are mainly controlled by the tide [12]. A principal component analysis was used by [11] to show that the effect of the tide and wind on the current could be separate in the data, as well as in the model results. A 3D simulation using ELCOM was done by forcing the boundary using predicted tides with satisfactory results [11], however for this area automated calibration has not been applied before as well the prediction capabilities of the model. This paper presents the details regarding how OpenDA is used to calibrate the parameters at any depth or friction, as well as the components of the tides in SQB. One of the objectives of this work is to show when the numerical model is unable to capture the physics of its region of study. In such cases, the correlation is low, and OpenDA plays an important role in the models calibration. The performance of the model is also quantified in terms of root-mean square errors (RMSE), and the correlation coefficients (r) follow the recommendations of the ASCE Task Committee [8].
2
Figure 1: Three-second resolution bathymetry interpolation of San Quintin Bay, B.C, Mexico figure 2 and the data available at each station is described as follows: • The first station which is named ID: North is located north of the bay, and comes from a predicted tide from MARV0.9 [15]. Since there are few stations to measure, this data is used as a observed data to calibrate Delft3D model for SQ. • The second station which is named ID: Met, is located in the north-west section of the bay, in the N 30.43358 W 115.98988 position. On this site, a weather station was available providing the atmospheric data such as air temperature, wind velocity components, atmospheric pressure, humidity, and rainfall gathered during the 2010 campaign. From this station wind fields were used to force the 3D Model. • The last station located at the bay entrance is named ID: Mouth and the data is provided by the Mexican Navy Secretary. The data contains hydrostatic pressure used for calibration and currents to validate the 3D run. From hydrostatic pressure, it obtained water levels at this location.
2. Study Region and Monitoring Data The bathymetry far from the coast were taken from ETOPO2, while the bathymetry near the shore were taken during campaign of June 2005 [13]. The topographic component and coastal lines were recorded from the digital fusion of level curves, at scale 1:50.000, as generated by the National Geophysical Data Center (NGDC). To optimize our data the Kriging, an interpolation method was applied[14] . The final topographic result is the generation of a barometric digital model with a special resolution of 3 seconds in the latitude and longitude in the horizontal, and 15 sigma-levels in the vertical. The results of the interpolation can be seen in figure 2. To validate the results, time series of water level and currents were used. The curren data were provided by the Navy Secretary from Mexico (SEMAR). This data was taken from the wet season of 2010 starting on August 27th to September 20th with 30 minutes intervals. Their locations are shown in
Before we could do any calibration, the data was treated to remove the noise and outliers due to equipment error. In order to smooth the velocities fields series at the ID:Mouth location, a moving average of 5 steps was applied. For hydrostatic pressure the seasonality in the data was removed to obtain water level. 3. Model Set up 3.1. Hydrodynamics Model The numerical model used in the present study is Delft3D, which has a long history of successful application to oceanic, coastal, and estuarine waters. This model is a multi-dimensional 2
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
3
Figure 2: A zoom in of the study region, including location points for the observed data stations: North, Mouth, and Met. San Quintin Bay, Mexico. (2D or 3D) hydrodynamic (and transport) simulation program which calculates non-steady flow that result from tidal and meteorological forcing on a boundary fitted grid. In 3D simulations, the vertical grid is defined following the sigma coordinates approach [16]. The numerical model solves a coupled system of differential, prognostic equations which describe the conservation of mass, momentum, heat, and salinity. Delft3DFLOW solves the Navier Stokes equations for an incompressible fluid under shallow water, as well as the Boussinesq assumptions. In the vertical momentum equation, the vertical accelerations are neglected, which then leads to the hydro-static pressure equation. In 3D models, the vertical velocities are computed from the continuity equation. This set of partial differential equations, in combination with an appropriate set of both initial and boundary conditions, is solved on a finite difference grid. (A detailed description of the model equation and its numerical algorithm can be found in [16]). Delft Dashboard was used to set up the numerical domain, initial and boundary conditions, and time discretization. Delft Dashboard is a standalone, Matlab-based graphical user interface [16]. The interface is coupled with the Delft3D modeling suite (Deltares), which allows users to carry out hydrodynamic, wave, morphodynamic, and water quality computations.
position from which the tide is observed. The general formula for the astronomical tide [16] is:
H (t) = A0 +
k X i=1
Ai Fi cos(ωi t + (V0 + u)i − Gi )
(1)
in which: H(t) is the water level at time t, A0 is the mean water level over a specific period number of relevant constituents, i is the index of a constituent, Ai the local tidal amplitude of a constituent, Fi is the nodal amplitude factor, ωi represent the angular velocity, (V0 + u)i are the astronomical arguments and Gi is the improved kappa number (= local phase lag). Fi and (V0 + u)i are time-dependent factors which, together with ωi , can easily be calculated and are generally tabulated in the various tidal year books. V0 is the phase correction factor which relates the local time frame of the observations to an internationally agreed upon celestial time frame. V0 is the frequency dependent and Fi and ui are the slowly varying amplitude and phase corrections (also frequency dependent). For this study, the astronomical tide, H(t), due to the eight primary harmonic constituents (A and G) of the semidiurnals (M2, S 2, N2, K2) and the diurnal (K1, O1, P1, Q1) was obtained from a global model of ocean tides. This global inverse tide model (TPX0 7.2) was developed by Oregon State University, and is used to initiate Delft3D [17]. On the other hand, a classical tidal harmonic analysis including error estimates using [18] was done to the pressure time series at the bay entry in order to find the more predominant components of the tide in the region (results are shown in Figure 3). An estimation of the tidal amplitude and phase with
3.2. Model Forcing Function Boundary conditions can specified in terms of astronomical components. The observed tidal motion can be described in terms of a series of simple harmonic constituent motions, each with its own characteristic frequency (angular velocity). The amplitudes A and phases G of the constituents vary with the 3
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
4
starting point, in such a way that they produce the corresponding results within state variables (the procedure in empirical models). When working with an empirical model, the values of the parameters must be calibrated, beginning with a sample of the input values, and output of the model and a purpose function (which value should be minimized). In the present study, OpenDA′ s semi-automated parameter estimation method, called Does Not Use Derivatives (DUD), is utilized. DUD is a derivative-free algorithm for nonlinear least squares [19]. It evaluates and optimizes uncertain model parameters by minimizing a generalized form of a least squares or goodness of fit (GoF) criterion which is formulated in the time domain. DUD is an optimization algorithm which does not use any derivative of the function being evaluated. It can be seen as a Gauss-Newton method, in the sense that it transforms the nonlinear least square problem into a well-known linear square problem. The difference is that, instead of approximating the nonlinear function by its tangent function, DUD uses an affine function for the linearization. A general idea of the DUD implementation is given by the following: Figure 3: Spectrum analysis for water level at the bay entry. Phase (Top) and Amplitude (Bottom) estimation of the NeapSpring significant tidal constituents.
• Start running first guess and modify each parameter;
95 % interval confidence shows that the predominant components are K1, 01, M2 and S 2 . Table 1 shows the parameter values phase and amplitude estimated from observation, as well as those used as an open boundary condition for the SQB model. This information is used to compare the improvement of the calibration within these parameters.
• If this is an improvement, update linearization with new point;
• Linearize the model around these values and solve linear problem;
• Or, do a line-search (only until there is improvement). In tidal modeling, the water level is the key model variable [4] and the built-in GoF formulation is now configured to read as:
3.3. Optimization and Calibration
Parameter M2 S2 K2 K1 O1 M4 MS 4
MAT Tide Observer Amplitude Phase 0.538 134.436 0.3147 131.844 0.0836 172.0200 0.2297 214.343 0.2025 190.2559 0.0087 292.2344 0.0046 327.7652
max max X max X 2 1X obs sim (tn ) /σ2Hobs (2) (tn ) − Hr,s,n Hr,s,n 2 r=1 s=1 n=1
R
Once a model has been chosen, it must be applied to a specific problem. First, the values of the parameters required for the model must be obtained. These can be measured directly within the field (the procedure in a physical model), or can be obtained via the use of optimization techniques. In the latter case, we begin with the unknown values of the variables at the
Go f =
S
N
in which H(t) is the water level measured at time t, sim refers to results obtained from model simulations over the simulations period [0, T ], obs are observed values, Nmax is the number of time steps in the time series, S max is the number of stations in region r, Rmax indicates the blocks for which observations are included, while σ2Hobs denotes the uncertainties assigned to the observations (here: tidal prediction values). σ2Hobs is set at 0.05 m. After calibrating the model, the parameter values should have a physical meaning, however if this is not the case, it could be that the model has predictive influence for the set of data utilized in the calibration.
TPXO 7.2 Amplitude Phase 0.485 129.599 0.215 127.896 0.062 121.610 0.308 203.941 0.196 189.309 0.003 272.662 0.003 301.000
4. Numerical Results
Table 1: Comparison of significant tidal constituents from observation at the Bay Entry and model tidal TPXO 7.2 at the boundary.
The model grid was Cartesian, with 15 layers of sigma coordinates, and a 3-second resolution in the horizontal ( 90 m). Since the main force in the bay is coming from the tide one major boundary-forcing function is applied at the west side of the 4
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
domain, and the global inverse tide model (TPX0 7.2) is used to initiate Delft3D. Since the entrance to the bay is protected by a sand bar of 3 km long, the open boundary was selected to be far from the entrance to allow the tide to approach the lagoon in a natural way. To begin the computations, it is necessary to specify the initial conditions for the water level elevation, velocity and wind force. The initial water surface was assumed to be horizontal, and the velocity components were set to zero throughout the model domain. To run the model proceeded in the following manner:
5
parameter the standard deviation chosen was 0.002. For depths, the standard deviation chosen was 0.05 m, and the new value calculated for the parameter will be added in the entire domain. For tide constituents, the standard deviation chosen for the amplitude was 0.10 m, and the phase was 10 degrees. It is important to note from table 2 that even small changes in the tidal constituents at the boundary can improve the model considerably. The maximum variation in amplitude were 5 degrees for O1 and 0.29 m in phase for M2 For roughness and depths, improvement was small; maximum of 1.5 % in the two stations (Table 3). Since the decrease of uncertainty in the model resulting from the calibration of the friction and depth parameters was not as considerable as expected and the semidiurnal tide is the main forcing of the hydrodynamics on the bay in a two-week period, it is necessary to extend the simulation to a longer period, we then proceeded to calibrate the semidiurnal and diurnal tidal constituents to determine if there was any significant improvement. In fact, there was a noteworthy reduction in the MSRE (see table 3), for station ID: Mouth (up to 12 %) and station ID: North showed an improvement of up to 101 %. Friction was also added to this experiment by changing the roughness coefficient, but the variability did not decrease compared with the last experiment. In both cases, 7 and 14 days simulations, it was observed that the calibration was better adjusted in the north station than in the entrance station of the bay. This may be the result of a bias between observation and the model. A final test was implemented to remove the bias, the bias is computed explicitly and removed from the cost function, result in table 3 shown an improvement of 26.83 % at the mouth station and a 169.48 % at the north station. It is important to note from table 2 that even small changes in the tidal constituents at the boundary can improve the model considerably. The maximum variation in amplitude were 5 degrees for and 0.29 m in phase for After the model was satisfactorily calibrated by removing the bias, a final forecast of 10 days. After the model was satisfactorily calibrated by removing the bias, a final forecast of 10 days (Figure 4) shows how well the model captures the dynamics for Water Level. At North Station the model did very well, but for south station it looks as if there is still some bias between the model and observation. It could be the manner in which the Water Level was taken from hydro-static pressure at ID:Mouth station, more than the calibration itself. To validate the model, the correlation coefficient was calculated for each station. The original model, after calibration and forecast against observations, (Table 4) shows a
1. Run 2D mode: Water Level • 7 days initial guess run, check friction and depths for calibration. • Extend to 14 days calibration for tide constituents at the boundaries. • Chose the best calibrated model.
• 10 days forecast with the best calibrated model. • Validate calibration.
2. Run the 3D mode: Velocities fields • 7 days Initial guess adding wind.
• 14 days calibrated model using the parameters obtained in the 2D run. • 10 days forecast with the calibrated model.
• Compare Initial run with calibrated model and forecast. 4.1. Water Level For water level calibration, the stratification can be neglected at present and the model run in 2D mode. The parameters to consider are the friction, bathymetry, and tidal boundary conditions. The temporary calibration was done as follows: to save computational time, 7 days of simulation was run first, and the friction and depth parameters were calibrated using the data from north of the bay and at the bay entrance. The 7 days of simulation lasted 30 minutes, and it was very simple to run this first set of tests. Also, to minimize uncertainty in the optimization process, the calibration was begun one day after the model was started. The calibration using OpenDA is based on an iterative process of optimization until the parameters of defined tolerance are completed. In our study, the criteria were as follows: • stop criterion 1, imain > maxit = 10 • stop criterion 2, diff < abstol = 1.0 • stop criterion 3, relDiff < reltol = 0.01
Parameter M2 S2 K1 O1
• stop criterion 4, linearized cost relative error < 0.01 To apply the DUD algorithm, first it is necessary to predetermine an observation standard deviation to the GoF function. For calibration of roughness OpenDA changes the Manning coefficient values in the regions specified and parameters are modified with a percentage that is constant within a region. For this
TPXO 7.2 Amplitude Phase 0.485 129.599 0.215 127.896 0.308 203.941 0.196 189.309
After Calibration Amplitude Phase 0.456 124.792 0.200 124.907 0.316 203.348 0.195 183.758
Table 2: Parameters changes before and after calibration at the boundaries. 5
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
% ObsID
Mouth.waterlevel
North.waterlevel
Simulation 2010 Campaign (2D Mode) Parameter Before Calibration (14 days) Depth (7 days ) Roughness-u/roughness-v (7 Days) M2, S2, K1, 01 (14 Days) Roughness, M2, S2, K1, 01 (14 Days) Roughness, M2, S2, K1, 01, Bias Removed (14 Days) Forecast (10 Days) Before Calibration (14 days) Depth (7 Days) Roughness-u/roughness-v (7 Days) M2, S2, K1, 01 (14 Days) Roughness, M2, S2, K1, 01 (14 Days) Roughness, M2, S2, K1,01, Bias Removed, (14 Days) Forecast (10 Days)
STD 0.1681 0.0385 0.02962 0.1446 0.1446 0.1449 0.0567 0.1394 0.0465 0.0437 0.0557 0.0554 0.0539 0.1227
6
RMSE 0.1837 0.15630 0.1549 0.16327 0.16329 0.1449 0.1931 0.1451 0.088350 0.08487 0.0720 0.0720 0.0538 0.1451
IMP % n/a 0.308153 1.168351 12.54576 12.53115 26.8300 n/a n/a 1.7113 5.8739 101.4792 101.7253 169.4833 n/a
Table 3: Sensitivity results of calibrating parameters depth, roughness-u/roughness-v, tidal constituents M2, S2, K1, 01. Observation Vs Model Output. Water Level SQ Bay 2010 at NORTH Station Observatios Vs Model Output
Observed − North Station Model Calibrated Model Forecast
1.5
WL (m) North
1 0.5 0 −0.5 −1 −1.5 27−Aug
29−Aug
31−Aug
02−Sep
04−Sep
07−Sep 09−Sep time−−−−>
11−Sep
13−Sep
Water Level SQ Bay 2010 at SOUTH Station Observatios Vs Model Output
15−Sep
18−Sep
20−Sep
Observed − South Station Model Calibrated Model Forecast
1.5
WL (m) Bay Entry
1 0.5 0 −0.5 −1 −1.5 27−Aug
29−Aug
31−Aug
02−Sep
04−Sep
07−Sep 09−Sep time−−−−>
11−Sep
13−Sep
15−Sep
18−Sep
20−Sep
Figure 4: Comparison of instantaneous water levels at the North station (Top Figure) and Bay Entry (Bottom Figure). 14 days simulation from the calibrated model starting on Aug 27, 2010 , and 10 days forecast after calibration is done (black line) very high correlation, increasing after calibrations and for prediction as well at the North station. The model gets better until 0.99 correlation, almost linear for forecast.
previous research [11] wind was found to be a predominant component in the dynamics of the bay. Wind from ID: Met station was set up in the model. The observed currents data was taken from station ID: Mouth measured at different depths. A vertical depth average will be used to validate the model against observation. Since seven days of simulation in 3D mode takes 8 hours of real-time in a single machine, we validate the first run only with 7 days for validation before any calibration is done, then run the model previously calibrated until September 20th. Figure 5 shows a comparison of instantaneous depth aver-
4.2. Velocity Field In order to create the 3D model for our study, 7 days of simulation were run with boundary conditions determined prior to calibration. As the major dynamic within the bay is a reflected depth between 1 and 15 meters, the 15 layers of sigma coordinates were set up every 1 m (for the first 10 meters). From 6
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
Model Set Up 2D Mode 3D Mode
Before Calibration r 0.9559 0.9681 0.8064 0.8185
Station Mouth.waterlevel North.waterlevel V x Depth average Mouth.velocity Vy Depth average Mouth.velocity
Calibrated Model r 0.9617 0.9759 0.8777 0.7610
7
Forecast r 0.9595 0.9925 0.8862 0.86441
Table 4: Correlation Coefficient r for Model Output vs Observation. Before calibration, calibrated model and forecast for water level and depth average velocities fields age velocity V x Eastward and Vy Northward with depth average currents observations, and also the forecast performed. To validate the model, a linear fit was performed for the components of velocity in the eastern and northern directions, Figure 6. At the eastward end, the model and data are closest to a linear fit more than at the northward. It can be observed from table 4 that the correlation coefficient is increasing after calibration and for forecast as well at the east direction. For northward velocity, results before and after calibration were a little different but still highly correlated, the correlation coefficient moved from 0.81 previously to calibration to 0.76 for the calibrated model, and the 0.86 for forecast, it could be model is taking more time to speed up while is forcing with wind, so in this case forecast perform much better after calibrated model also by the time of calibrated model the observations has some outliers during that period of time that we do not removed by smoothing, this could be a reason why after calibration at east direction model correlation coefficients is decreasing. To verify if the model was correctly explaining the dynamic of the bay we revise the vertical depth average velocity within a certain determined time. In this case, it was taken at the 10th hour on the 26th of August during low tide. It can be observed that the current is headed toward the outside of the bay, and that the maximum velocities are observed in the deepest sections. During high tide (hour 17) of the same day, it is observed that the water level along the bay has maximum velocities in the deepest section ( Figures 7a and 7b). The velocity field obtained by modeling the flux of the selected in flowing tide shows that tide forces contributed to the reduction of the water velocity over the shallowest part of the area. It also shows that maintaining higher water velocities through the channels that define the eastern and western limits of the bay. To find out the main direction and maximum velocities at the south station, a comparison model between observation during forecast period is shown in Figure 8. During that period of time, the maximum velocities are between 0.5 and 0.6 ms. This plot also shows two main directions: the first one, while the tide comes in and the second one when tide goes out. Also we observed an increasing in magnitude for low tides. The model was able to capture the same behavior as observed data with some a small over-prediction of 0.1 ms difference. The magnitude of the velocities in the northern direction was greater than in the eastern direction.
Conclusion The model results show that the Delft3D model is capable of reproduce the essential processes in the San Quintin Bay, and can be forced by the tidal model and wind. Delft3D in 2D mode was able to capture changes in water level with a high correlation in the stations where we had observations with very high predictability capabilities. Calibration using OpenDA in a small scale problem, did improve the first model implementation showing that friction and changes in depth were not significant. Conversely, diurnal and semidiurnal components at the boundaries show that a change in their parameter have an important role in the dynamic of the bay. Better performance was shown after the model was calibrated. At north station, it shows a 99 % linear correlation between model and observation. However at the bay entry, measurement for water level prediction capabilities were not as big as expected but still highly correlated. This could be a product of measurement error from the equipment more than the model itself. While tide and wind forcing are the main forces that move the hydrodynamics inside the bay, the depths average velocities for the 3D runs shows that Delft3D and OpenDA calibration in 2D was enough to explain the variability in currents along the bay with a very high correlation during the period of time considered, and there was no need to calibrate any other parameters for the 3D model. In this context predictability capabilities increase over time due the increased speed of the model to capture the dynamics. At the bay entry we found the main direction is NNE during high tide and SSW during low tide with a higher magnitude in this direction. A highlight from this work could be that, for this area, automated calibration has not been applied before as well the forecast capabilities of the model in 3D mode. Implementing Open DA as a tool for calibration in a small scale problem is also a new approach with successful results. In this aspect, the next step will be to consider OpenDA as a tool to deliver real time forecasting capabilities in this region by introducing data assimilation techniques. However new computational challenges need to be addressed to reduce computational costs of these implementations. Knowing the development of these variables over time is very important in the study of growing oysters Crassotrea gigas inside the bay. The 3D output from the model plays a very important role, since it is possible to determine which shallow areas inside the bay has enough current to transport nutrients 7
y component of depth averaged velocity (m/s)
x component of depth averaged velocity (m/s)
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
8
Depth average Velocity EASTWARD (V ) SQ Bay 2010 at SOUTH Station x Observatios Vs Model Output Observed − South Station Model Calibrated Model Forecast
1
0.5
0
−0.5
−1 27−Aug
29−Aug
31−Aug
02−Sep
04−Sep
07−Sep 09−Sep time−−−−>
11−Sep
13−Sep
15−Sep
18−Sep
20−Sep
Depth average Velocity NORTHTWARD (V ) SQ Bay 2010 at SOUTH Station y Observatios Vs Model Output
Observed − South Station Model Calibrated Model Forecast
1
0.5
0
−0.5
−1 27−Aug
29−Aug
31−Aug
02−Sep
04−Sep
07−Sep 09−Sep time−−−−>
11−Sep
13−Sep
15−Sep
18−Sep
20−Sep
Figure 5: Comparison of instantaneous depth average velocity V x Eastward (Top Figure) and Vy Northward (Bottom Figure) at the Bay Entry. 14 days simulation from the calibrated model starting on Aug 27, 2010, and 10 days forecast after calibration is done (black line)
ORIGINAL Model Eastward Velocity (Ux) R=0.80648 RMSE= 0.058906
ORIGINAL Model Northward Velocity (Uy) R=0.81851 RMSE= 0.1751 0.5
Uy Model
Ux Model
0.5
0
−0.5 −0.2
0
−0.5 −0.15
−0.1
−0.05
0
Ux Observed
0.05
0.1
0.15
0.2
−0.8
−0.6
CALIBRATED Model Eastward Velocity (Ux) R=0.87777 RMSE= 0.10146
0
−0.5
0.2
0.4
0.6
0.8
0.6
0.8
0.6
0.8
−0.5 −0.15
−0.1
−0.05
0
Ux Observed
0.05
0.1
0.15
0.2
−0.8
−0.6
−0.4
−0.2
0
Uy Observed
0.2
0.4
FORECAST Model Northward Velocity (Uy) R=0.86441 RMSE= 0.21025 0.5
Uy Model
0.5
Ux Model
0
Uy Observed
0
FORECAST Model Eastward Velocity (Ux) R=0.88622 RMSE= 0.11834
0
−0.5 −0.2
−0.2
0.5
Uy Model
Ux Model
0.5
−0.2
−0.4
CALIBRATED Model Northward Velocity (Uy) R=0.76106 RMSE= 0.18155
0
−0.5 −0.15
−0.1
−0.05
0
Ux Observed
0.05
0.1
0.15
0.2
−0.8
−0.6
−0.4
−0.2
0
Uy Observed
0.2
0.4
Figure 6: Linear regression between observations and 3D Model, for vertical depth average of Vx and Vy, before calibration, calibrated model and forecast. 8
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
for healthy cultivation of the oysters. Additionally, it is important to note that the quality of the baroclinic structure strongly depends upon the quality of the boundary conditions for temperature and salinity in this type of scenario. A future work will study temperature variation along the bay to determine the best cultivation zones for these organisms.
9
[19] M. L. Raltson, R. I. Jennrich, DUD, A derivative-free algorithm for nonlinear least squares., Technometrics 20 (1) (1978) 7–14.
Acknowledgements Special thanks to the Navy Secretary of Mexico (SEMAR), for the use of the currents data at the entrance of the bay, as well as Dr. Rafael Blanco from UABC-Mexico and Jody Payne for their invaluable support during this research. References [1] E. S. Ghada, M. Verlaan, S. Hummel, A. Weerts, J. Dhondia, OpenDA open source generic data assimilation environment and its application in process models., Geophysical Research Abstracts 12: EGU2010-9346-2. URL http://adsabs.harvard.edu/abs/2010EGUGA..12.9346E [2] G. Y. El Sarafy, H. Gerritsen, S. Hummel, A. H. Weerts, A. E. Mynett, M. Tanaka, Application of data assimilation in portable operational forecasting systems the DATools assimilation environment, Ocean Dynamics 57 (4-5) (2007) 485–499. [3] A. H. Weerts, G. Y. El Serafy, S. Hummel, J. Dhondia, H. Gerritsen, Application of generic data assimilation tools (DATools) for flood forecasting purposes, Computers Geosciences 36 (4) (2010) 453–463. [4] A. Kurniawan, H. Gerritsen, S. K. Ooi, S. Hummel, Sensitivity analysis of the tidal representation in Singapore Regional Waters in a data assimilation environment, Ocean Dynamics 61 (2011) 112–1236. [5] F. Zijl, M. Verlaan, H. Gerritsen, Improved water-level forecasting for the northwest european shelf and north sea through direct modelling of tide, surge and non-linear interaction, Ocean Dynamics 63 (7) (2013) 823–847. doi:10.1007/s10236-013-0624-2. URL http://dx.doi.org/10.1007/s10236-013-0624-2 [6] P. Chavez, M. Acosta, Desarrollo gonadal de crassostrea gigas en Bahia San Quintin, Baja California, Mexico. , Ciencias Marinas 20 (2) (1995) 225–242. [7] T. Imai, S. Sakai, Study of breeding of Japanese oyster, Crassostrea gigas., Agric. Res. 12 (1961) 125–171. [8] J. D. Wang, A. F. Blumberg, H. L. Butler, ASC, P. Hamilton, Transport Prediction in Partially Stratified Tidal Waters, Journal of Hydraulic Engineering 116 (1990) 380–396. [9] I. Del Valle Lucero, H. R. Cabrera Muro, Aplicacion de un modelo numerico unidimensional a bahia San Quintin, B.C., Ciencias Marinas 7 (1) (1981) 1–15. [10] O. E. Delgado-Gonzalez, J. A. Jimenez, J. L. Ferman-Almada, G.-E. Z, F. Marvan-Gargollo, A. Mejia-Trejo, Depth and hydrodynamics as tools to select aquaculture areas in the M Depth and hydrodynamics as tools to select aquaculture areas in the coastal zone, Ciencias Marinas 36 (3) (2010) 249–265. [11] I. Ramirez, R. Blanco, R. Vaquez, G. Ramirez, The simulation of the circulation of San Quintin Bay., Submited to Ocean Dynamics. [12] J. Jimenez, A. Mejia-Trejo, Investigacion numerica de la dinamica de la circulacion en bahia de san quintin, b. c.: procesos de transporte, Vol. 12, XII Congreso Nacional de Oceanograf´ıa, 2000. [13] I. Ramirez, D. Vincelli, Mediciones oceanograficas en bahia san quintn, b.c, mexico., CICESE, TECH. [14] Kriging [online]. [15] J. I. Gonzalez. MAR V1.0 para Windows 98/2000/XP/Vista/Windows7 [online] (Sep. 2013). [16] H. Delft, User manual Delft3D-Flow, DELTARES (Sep. 2011). [17] G. D. Egbert, A. F. Bennett, ‘TOPEX/ POSEIDON tides estimated using a global inverse model.’, J. Geophys. Res (1994) 821–824. [18] RIKZ, User’s Guide TIDEMAT tidal analysis, ALKYON Hydraulic Consultancy and Research, PO Box 248 8300 AE, Emmeloord The Netherlands (Jan. 2004).
9
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
10
Velocities Measured Bay Entry 90
1
120
60
0.5
150
30
180
0
330
210
240
300 270
(a) during high tide
Velocities Model Bay Entry 90
1
120
60
0.5
150
30
180
0
330
210
240
300 270
(b) during low tide
Figure 7: Depth Average Velocity, Magnitude (m/s)
Figure 8: Classical compass for current, measured data (Top Figure) and output from 3D model (Bottom Figure) during forecast period
10
S0377-0427(14)00236-2 http://dx.doi.org/10.1016/j.cam.2014.05.003 CAM 9671
To appear in:
Journal of Computational and Applied Mathematics
Received date: 1 October 2013 Revised date: 28 April 2014 Please cite this article as: M. Garcia, I. Ramirez, M. Verlaan, J. Castillo, Application of a three-dimensional hydrodynamic model for San Quintin Bay, B.C., Mexico. Validation and calibration using OpenDA, Journal of Computational and Applied Mathematics (2014), http://dx.doi.org/10.1016/j.cam.2014.05.003 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
mariangelapril2014.tex Click here to view linked References
Journal of Computational and Applied Mathematics 00 (2014) 1–9
Comp. and Applied Math.
4th International Congress on Computational Engineering and Sciences
Application of a Three-Dimensional Hydrodynamic Model for San Quintin Bay, B.C., Mexico. Validation and Calibration Using OpenDA. Mariangel Garciaa,1, Isabel Ramirezb,2, Martin Verlaanc,3, Jose Castilloa,4 a Computational b Centro
Sciences Research Center, San Diego State University, San Diego, C.A, U.S de Investigacion y Educacion Superior de Ensenada, B.C, Mexico c Deltares, Delft, The Netherlands
Abstract A 3D hydrodynamic model (Delf3D) was applied for San Quintin Bay (SQB), B.C., Mexico. Calibration and validation were conducted using measured bathymetry, wind, water surface elevation and velocities fields. The calibration period took place during the winter of 2010. Model predictions were evaluated graphically and statistically against field observations in order to quantify the accuracy of the predictions and evaluate the success of the model calibration. A 2D Calibration of tide constituents at the boundaries, seabed roughness and depths was done to evaluate the improvement in the velocities for the 3D Model. Comparisons between the model predictions and the field observations of the water surface elevation at interior stations indicated that the model had been successfully calibrated. It also indicated that the model predictions were highly correlated with the observed water surface elevations and velocities. This connection between the observed and simulated values was detected through graphical comparisons and root-mean-square errors and correlation coefficient. The objective of this study was to demonstrate that OpenDA could be used effectively to rapidly calibrate a hydrodynamic model for smaller sub-meso scales like the size of SQ Bay. Keywords: Calibration, OpenDA, Data Assimilation, Delft3D, San Quitin Bay, Hydrodynamic.
1. Introduction
improving the starting position of the forcing for a forecast, so that the estimates are different in each cycle. OpenDA software was designed to quickly and efficiently facilitate the use of a tool which enables the reduction of uncertainty of a given model. The generic OpenDA data assimilation environment [1] provides both a range of filtering routines, as well as uncertainty and sensitivity analysis routines. It also supports assimilation of recent observations to improve model forecasts, and allows the user to carry out sensitivity analysis and simultaneous parameter optimization of model parameters in a given model application. It is an updated, open-source version which includes the extended functionalities of the existing DATools system, which has been successfully used in data assimilation of current and salinity profiles [2], for flood forecasting purposes [3] as well as for tidal sensitivity analysis[4]. However OpenDA has been applied before to calibration of coastal seas, but not yet to the submeso scales. A sensitivity analysis of the tidal representation in Singapore Regional Waters was done using 2D WAQUA by [5] but the boundary conditions were not calibrated. In this paper we are using Delft3D for the hydrodynamics model focus in a local scale problem where the open-boundary plays a very important role as well as the dynamics at the entrance.
Numerical ocean models have proven to be essential for the creation of forecast estimates, as well as weather forecasts for coastal areas. However, model errors are currently inevitable, due to uncertainties in both the initial and boundary conditions, as well as in the oversimplified or entirely neglected physical processes and mathematical approximations used in the systems equations. Furthermore, despite numerical weather prediction applications, available coastal ocean observation data is relatively scarce, and often does not capture complex coastal dynamics. Measurement error is another source of variability that is often not well understood, making it difficult to characterize them. Two methods of reducing the uncertainty within the models are presented here. The first method are calibration, with a goal of tuning a set of parameters that are fixed in time, for example: depth, amplitude, phase of tidal constituents and roughness. The second method is data assimilation, with a goal of 1 [email protected] 2 [email protected] 3 [email protected] 4 [email protected]
1
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
The region of study selected was San Quintin Bay (SQB), located South of Baja California., Mexico. SQB forms an interesting ecosystem, which lends itself for the development of aquaculture. In particular, excellent results have been obtained in the growth and survival of the Japanese oyster Crassotrea gigas [6]. One of the most important biological processes for cultivating this organism is its reproductive cycle, which is principally regulated by the amount and quality of its food, as well as by water temperature [7]. The hydrodynamic must be known with accuracy in order to ensure efficient ecological monitoring of the bay. The hydrodynamics of SQB has been studied by [9], [10], and [11]. One of the first results found was by [9] who applied a one-dimensional numerical model to reproduce the hydrodynamics of the lagoon. Observations of tides, currents, temperatures, salinity and winds were used in that research, however the numerical model that was used proved to be of little predictive value for this location. Similar values calculated and observed values of water surface elevations and associated velocities were obtained by applying a two-dimensional hydrodynamic model, which determined that the dynamics of SQB are mainly controlled by the tide [12]. A principal component analysis was used by [11] to show that the effect of the tide and wind on the current could be separate in the data, as well as in the model results. A 3D simulation using ELCOM was done by forcing the boundary using predicted tides with satisfactory results [11], however for this area automated calibration has not been applied before as well the prediction capabilities of the model. This paper presents the details regarding how OpenDA is used to calibrate the parameters at any depth or friction, as well as the components of the tides in SQB. One of the objectives of this work is to show when the numerical model is unable to capture the physics of its region of study. In such cases, the correlation is low, and OpenDA plays an important role in the models calibration. The performance of the model is also quantified in terms of root-mean square errors (RMSE), and the correlation coefficients (r) follow the recommendations of the ASCE Task Committee [8].
2
Figure 1: Three-second resolution bathymetry interpolation of San Quintin Bay, B.C, Mexico figure 2 and the data available at each station is described as follows: • The first station which is named ID: North is located north of the bay, and comes from a predicted tide from MARV0.9 [15]. Since there are few stations to measure, this data is used as a observed data to calibrate Delft3D model for SQ. • The second station which is named ID: Met, is located in the north-west section of the bay, in the N 30.43358 W 115.98988 position. On this site, a weather station was available providing the atmospheric data such as air temperature, wind velocity components, atmospheric pressure, humidity, and rainfall gathered during the 2010 campaign. From this station wind fields were used to force the 3D Model. • The last station located at the bay entrance is named ID: Mouth and the data is provided by the Mexican Navy Secretary. The data contains hydrostatic pressure used for calibration and currents to validate the 3D run. From hydrostatic pressure, it obtained water levels at this location.
2. Study Region and Monitoring Data The bathymetry far from the coast were taken from ETOPO2, while the bathymetry near the shore were taken during campaign of June 2005 [13]. The topographic component and coastal lines were recorded from the digital fusion of level curves, at scale 1:50.000, as generated by the National Geophysical Data Center (NGDC). To optimize our data the Kriging, an interpolation method was applied[14] . The final topographic result is the generation of a barometric digital model with a special resolution of 3 seconds in the latitude and longitude in the horizontal, and 15 sigma-levels in the vertical. The results of the interpolation can be seen in figure 2. To validate the results, time series of water level and currents were used. The curren data were provided by the Navy Secretary from Mexico (SEMAR). This data was taken from the wet season of 2010 starting on August 27th to September 20th with 30 minutes intervals. Their locations are shown in
Before we could do any calibration, the data was treated to remove the noise and outliers due to equipment error. In order to smooth the velocities fields series at the ID:Mouth location, a moving average of 5 steps was applied. For hydrostatic pressure the seasonality in the data was removed to obtain water level. 3. Model Set up 3.1. Hydrodynamics Model The numerical model used in the present study is Delft3D, which has a long history of successful application to oceanic, coastal, and estuarine waters. This model is a multi-dimensional 2
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
3
Figure 2: A zoom in of the study region, including location points for the observed data stations: North, Mouth, and Met. San Quintin Bay, Mexico. (2D or 3D) hydrodynamic (and transport) simulation program which calculates non-steady flow that result from tidal and meteorological forcing on a boundary fitted grid. In 3D simulations, the vertical grid is defined following the sigma coordinates approach [16]. The numerical model solves a coupled system of differential, prognostic equations which describe the conservation of mass, momentum, heat, and salinity. Delft3DFLOW solves the Navier Stokes equations for an incompressible fluid under shallow water, as well as the Boussinesq assumptions. In the vertical momentum equation, the vertical accelerations are neglected, which then leads to the hydro-static pressure equation. In 3D models, the vertical velocities are computed from the continuity equation. This set of partial differential equations, in combination with an appropriate set of both initial and boundary conditions, is solved on a finite difference grid. (A detailed description of the model equation and its numerical algorithm can be found in [16]). Delft Dashboard was used to set up the numerical domain, initial and boundary conditions, and time discretization. Delft Dashboard is a standalone, Matlab-based graphical user interface [16]. The interface is coupled with the Delft3D modeling suite (Deltares), which allows users to carry out hydrodynamic, wave, morphodynamic, and water quality computations.
position from which the tide is observed. The general formula for the astronomical tide [16] is:
H (t) = A0 +
k X i=1
Ai Fi cos(ωi t + (V0 + u)i − Gi )
(1)
in which: H(t) is the water level at time t, A0 is the mean water level over a specific period number of relevant constituents, i is the index of a constituent, Ai the local tidal amplitude of a constituent, Fi is the nodal amplitude factor, ωi represent the angular velocity, (V0 + u)i are the astronomical arguments and Gi is the improved kappa number (= local phase lag). Fi and (V0 + u)i are time-dependent factors which, together with ωi , can easily be calculated and are generally tabulated in the various tidal year books. V0 is the phase correction factor which relates the local time frame of the observations to an internationally agreed upon celestial time frame. V0 is the frequency dependent and Fi and ui are the slowly varying amplitude and phase corrections (also frequency dependent). For this study, the astronomical tide, H(t), due to the eight primary harmonic constituents (A and G) of the semidiurnals (M2, S 2, N2, K2) and the diurnal (K1, O1, P1, Q1) was obtained from a global model of ocean tides. This global inverse tide model (TPX0 7.2) was developed by Oregon State University, and is used to initiate Delft3D [17]. On the other hand, a classical tidal harmonic analysis including error estimates using [18] was done to the pressure time series at the bay entry in order to find the more predominant components of the tide in the region (results are shown in Figure 3). An estimation of the tidal amplitude and phase with
3.2. Model Forcing Function Boundary conditions can specified in terms of astronomical components. The observed tidal motion can be described in terms of a series of simple harmonic constituent motions, each with its own characteristic frequency (angular velocity). The amplitudes A and phases G of the constituents vary with the 3
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
4
starting point, in such a way that they produce the corresponding results within state variables (the procedure in empirical models). When working with an empirical model, the values of the parameters must be calibrated, beginning with a sample of the input values, and output of the model and a purpose function (which value should be minimized). In the present study, OpenDA′ s semi-automated parameter estimation method, called Does Not Use Derivatives (DUD), is utilized. DUD is a derivative-free algorithm for nonlinear least squares [19]. It evaluates and optimizes uncertain model parameters by minimizing a generalized form of a least squares or goodness of fit (GoF) criterion which is formulated in the time domain. DUD is an optimization algorithm which does not use any derivative of the function being evaluated. It can be seen as a Gauss-Newton method, in the sense that it transforms the nonlinear least square problem into a well-known linear square problem. The difference is that, instead of approximating the nonlinear function by its tangent function, DUD uses an affine function for the linearization. A general idea of the DUD implementation is given by the following: Figure 3: Spectrum analysis for water level at the bay entry. Phase (Top) and Amplitude (Bottom) estimation of the NeapSpring significant tidal constituents.
• Start running first guess and modify each parameter;
95 % interval confidence shows that the predominant components are K1, 01, M2 and S 2 . Table 1 shows the parameter values phase and amplitude estimated from observation, as well as those used as an open boundary condition for the SQB model. This information is used to compare the improvement of the calibration within these parameters.
• If this is an improvement, update linearization with new point;
• Linearize the model around these values and solve linear problem;
• Or, do a line-search (only until there is improvement). In tidal modeling, the water level is the key model variable [4] and the built-in GoF formulation is now configured to read as:
3.3. Optimization and Calibration
Parameter M2 S2 K2 K1 O1 M4 MS 4
MAT Tide Observer Amplitude Phase 0.538 134.436 0.3147 131.844 0.0836 172.0200 0.2297 214.343 0.2025 190.2559 0.0087 292.2344 0.0046 327.7652
max max X max X 2 1X obs sim (tn ) /σ2Hobs (2) (tn ) − Hr,s,n Hr,s,n 2 r=1 s=1 n=1
R
Once a model has been chosen, it must be applied to a specific problem. First, the values of the parameters required for the model must be obtained. These can be measured directly within the field (the procedure in a physical model), or can be obtained via the use of optimization techniques. In the latter case, we begin with the unknown values of the variables at the
Go f =
S
N
in which H(t) is the water level measured at time t, sim refers to results obtained from model simulations over the simulations period [0, T ], obs are observed values, Nmax is the number of time steps in the time series, S max is the number of stations in region r, Rmax indicates the blocks for which observations are included, while σ2Hobs denotes the uncertainties assigned to the observations (here: tidal prediction values). σ2Hobs is set at 0.05 m. After calibrating the model, the parameter values should have a physical meaning, however if this is not the case, it could be that the model has predictive influence for the set of data utilized in the calibration.
TPXO 7.2 Amplitude Phase 0.485 129.599 0.215 127.896 0.062 121.610 0.308 203.941 0.196 189.309 0.003 272.662 0.003 301.000
4. Numerical Results
Table 1: Comparison of significant tidal constituents from observation at the Bay Entry and model tidal TPXO 7.2 at the boundary.
The model grid was Cartesian, with 15 layers of sigma coordinates, and a 3-second resolution in the horizontal ( 90 m). Since the main force in the bay is coming from the tide one major boundary-forcing function is applied at the west side of the 4
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
domain, and the global inverse tide model (TPX0 7.2) is used to initiate Delft3D. Since the entrance to the bay is protected by a sand bar of 3 km long, the open boundary was selected to be far from the entrance to allow the tide to approach the lagoon in a natural way. To begin the computations, it is necessary to specify the initial conditions for the water level elevation, velocity and wind force. The initial water surface was assumed to be horizontal, and the velocity components were set to zero throughout the model domain. To run the model proceeded in the following manner:
5
parameter the standard deviation chosen was 0.002. For depths, the standard deviation chosen was 0.05 m, and the new value calculated for the parameter will be added in the entire domain. For tide constituents, the standard deviation chosen for the amplitude was 0.10 m, and the phase was 10 degrees. It is important to note from table 2 that even small changes in the tidal constituents at the boundary can improve the model considerably. The maximum variation in amplitude were 5 degrees for O1 and 0.29 m in phase for M2 For roughness and depths, improvement was small; maximum of 1.5 % in the two stations (Table 3). Since the decrease of uncertainty in the model resulting from the calibration of the friction and depth parameters was not as considerable as expected and the semidiurnal tide is the main forcing of the hydrodynamics on the bay in a two-week period, it is necessary to extend the simulation to a longer period, we then proceeded to calibrate the semidiurnal and diurnal tidal constituents to determine if there was any significant improvement. In fact, there was a noteworthy reduction in the MSRE (see table 3), for station ID: Mouth (up to 12 %) and station ID: North showed an improvement of up to 101 %. Friction was also added to this experiment by changing the roughness coefficient, but the variability did not decrease compared with the last experiment. In both cases, 7 and 14 days simulations, it was observed that the calibration was better adjusted in the north station than in the entrance station of the bay. This may be the result of a bias between observation and the model. A final test was implemented to remove the bias, the bias is computed explicitly and removed from the cost function, result in table 3 shown an improvement of 26.83 % at the mouth station and a 169.48 % at the north station. It is important to note from table 2 that even small changes in the tidal constituents at the boundary can improve the model considerably. The maximum variation in amplitude were 5 degrees for and 0.29 m in phase for After the model was satisfactorily calibrated by removing the bias, a final forecast of 10 days. After the model was satisfactorily calibrated by removing the bias, a final forecast of 10 days (Figure 4) shows how well the model captures the dynamics for Water Level. At North Station the model did very well, but for south station it looks as if there is still some bias between the model and observation. It could be the manner in which the Water Level was taken from hydro-static pressure at ID:Mouth station, more than the calibration itself. To validate the model, the correlation coefficient was calculated for each station. The original model, after calibration and forecast against observations, (Table 4) shows a
1. Run 2D mode: Water Level • 7 days initial guess run, check friction and depths for calibration. • Extend to 14 days calibration for tide constituents at the boundaries. • Chose the best calibrated model.
• 10 days forecast with the best calibrated model. • Validate calibration.
2. Run the 3D mode: Velocities fields • 7 days Initial guess adding wind.
• 14 days calibrated model using the parameters obtained in the 2D run. • 10 days forecast with the calibrated model.
• Compare Initial run with calibrated model and forecast. 4.1. Water Level For water level calibration, the stratification can be neglected at present and the model run in 2D mode. The parameters to consider are the friction, bathymetry, and tidal boundary conditions. The temporary calibration was done as follows: to save computational time, 7 days of simulation was run first, and the friction and depth parameters were calibrated using the data from north of the bay and at the bay entrance. The 7 days of simulation lasted 30 minutes, and it was very simple to run this first set of tests. Also, to minimize uncertainty in the optimization process, the calibration was begun one day after the model was started. The calibration using OpenDA is based on an iterative process of optimization until the parameters of defined tolerance are completed. In our study, the criteria were as follows: • stop criterion 1, imain > maxit = 10 • stop criterion 2, diff < abstol = 1.0 • stop criterion 3, relDiff < reltol = 0.01
Parameter M2 S2 K1 O1
• stop criterion 4, linearized cost relative error < 0.01 To apply the DUD algorithm, first it is necessary to predetermine an observation standard deviation to the GoF function. For calibration of roughness OpenDA changes the Manning coefficient values in the regions specified and parameters are modified with a percentage that is constant within a region. For this
TPXO 7.2 Amplitude Phase 0.485 129.599 0.215 127.896 0.308 203.941 0.196 189.309
After Calibration Amplitude Phase 0.456 124.792 0.200 124.907 0.316 203.348 0.195 183.758
Table 2: Parameters changes before and after calibration at the boundaries. 5
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
% ObsID
Mouth.waterlevel
North.waterlevel
Simulation 2010 Campaign (2D Mode) Parameter Before Calibration (14 days) Depth (7 days ) Roughness-u/roughness-v (7 Days) M2, S2, K1, 01 (14 Days) Roughness, M2, S2, K1, 01 (14 Days) Roughness, M2, S2, K1, 01, Bias Removed (14 Days) Forecast (10 Days) Before Calibration (14 days) Depth (7 Days) Roughness-u/roughness-v (7 Days) M2, S2, K1, 01 (14 Days) Roughness, M2, S2, K1, 01 (14 Days) Roughness, M2, S2, K1,01, Bias Removed, (14 Days) Forecast (10 Days)
STD 0.1681 0.0385 0.02962 0.1446 0.1446 0.1449 0.0567 0.1394 0.0465 0.0437 0.0557 0.0554 0.0539 0.1227
6
RMSE 0.1837 0.15630 0.1549 0.16327 0.16329 0.1449 0.1931 0.1451 0.088350 0.08487 0.0720 0.0720 0.0538 0.1451
IMP % n/a 0.308153 1.168351 12.54576 12.53115 26.8300 n/a n/a 1.7113 5.8739 101.4792 101.7253 169.4833 n/a
Table 3: Sensitivity results of calibrating parameters depth, roughness-u/roughness-v, tidal constituents M2, S2, K1, 01. Observation Vs Model Output. Water Level SQ Bay 2010 at NORTH Station Observatios Vs Model Output
Observed − North Station Model Calibrated Model Forecast
1.5
WL (m) North
1 0.5 0 −0.5 −1 −1.5 27−Aug
29−Aug
31−Aug
02−Sep
04−Sep
07−Sep 09−Sep time−−−−>
11−Sep
13−Sep
Water Level SQ Bay 2010 at SOUTH Station Observatios Vs Model Output
15−Sep
18−Sep
20−Sep
Observed − South Station Model Calibrated Model Forecast
1.5
WL (m) Bay Entry
1 0.5 0 −0.5 −1 −1.5 27−Aug
29−Aug
31−Aug
02−Sep
04−Sep
07−Sep 09−Sep time−−−−>
11−Sep
13−Sep
15−Sep
18−Sep
20−Sep
Figure 4: Comparison of instantaneous water levels at the North station (Top Figure) and Bay Entry (Bottom Figure). 14 days simulation from the calibrated model starting on Aug 27, 2010 , and 10 days forecast after calibration is done (black line) very high correlation, increasing after calibrations and for prediction as well at the North station. The model gets better until 0.99 correlation, almost linear for forecast.
previous research [11] wind was found to be a predominant component in the dynamics of the bay. Wind from ID: Met station was set up in the model. The observed currents data was taken from station ID: Mouth measured at different depths. A vertical depth average will be used to validate the model against observation. Since seven days of simulation in 3D mode takes 8 hours of real-time in a single machine, we validate the first run only with 7 days for validation before any calibration is done, then run the model previously calibrated until September 20th. Figure 5 shows a comparison of instantaneous depth aver-
4.2. Velocity Field In order to create the 3D model for our study, 7 days of simulation were run with boundary conditions determined prior to calibration. As the major dynamic within the bay is a reflected depth between 1 and 15 meters, the 15 layers of sigma coordinates were set up every 1 m (for the first 10 meters). From 6
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
Model Set Up 2D Mode 3D Mode
Before Calibration r 0.9559 0.9681 0.8064 0.8185
Station Mouth.waterlevel North.waterlevel V x Depth average Mouth.velocity Vy Depth average Mouth.velocity
Calibrated Model r 0.9617 0.9759 0.8777 0.7610
7
Forecast r 0.9595 0.9925 0.8862 0.86441
Table 4: Correlation Coefficient r for Model Output vs Observation. Before calibration, calibrated model and forecast for water level and depth average velocities fields age velocity V x Eastward and Vy Northward with depth average currents observations, and also the forecast performed. To validate the model, a linear fit was performed for the components of velocity in the eastern and northern directions, Figure 6. At the eastward end, the model and data are closest to a linear fit more than at the northward. It can be observed from table 4 that the correlation coefficient is increasing after calibration and for forecast as well at the east direction. For northward velocity, results before and after calibration were a little different but still highly correlated, the correlation coefficient moved from 0.81 previously to calibration to 0.76 for the calibrated model, and the 0.86 for forecast, it could be model is taking more time to speed up while is forcing with wind, so in this case forecast perform much better after calibrated model also by the time of calibrated model the observations has some outliers during that period of time that we do not removed by smoothing, this could be a reason why after calibration at east direction model correlation coefficients is decreasing. To verify if the model was correctly explaining the dynamic of the bay we revise the vertical depth average velocity within a certain determined time. In this case, it was taken at the 10th hour on the 26th of August during low tide. It can be observed that the current is headed toward the outside of the bay, and that the maximum velocities are observed in the deepest sections. During high tide (hour 17) of the same day, it is observed that the water level along the bay has maximum velocities in the deepest section ( Figures 7a and 7b). The velocity field obtained by modeling the flux of the selected in flowing tide shows that tide forces contributed to the reduction of the water velocity over the shallowest part of the area. It also shows that maintaining higher water velocities through the channels that define the eastern and western limits of the bay. To find out the main direction and maximum velocities at the south station, a comparison model between observation during forecast period is shown in Figure 8. During that period of time, the maximum velocities are between 0.5 and 0.6 ms. This plot also shows two main directions: the first one, while the tide comes in and the second one when tide goes out. Also we observed an increasing in magnitude for low tides. The model was able to capture the same behavior as observed data with some a small over-prediction of 0.1 ms difference. The magnitude of the velocities in the northern direction was greater than in the eastern direction.
Conclusion The model results show that the Delft3D model is capable of reproduce the essential processes in the San Quintin Bay, and can be forced by the tidal model and wind. Delft3D in 2D mode was able to capture changes in water level with a high correlation in the stations where we had observations with very high predictability capabilities. Calibration using OpenDA in a small scale problem, did improve the first model implementation showing that friction and changes in depth were not significant. Conversely, diurnal and semidiurnal components at the boundaries show that a change in their parameter have an important role in the dynamic of the bay. Better performance was shown after the model was calibrated. At north station, it shows a 99 % linear correlation between model and observation. However at the bay entry, measurement for water level prediction capabilities were not as big as expected but still highly correlated. This could be a product of measurement error from the equipment more than the model itself. While tide and wind forcing are the main forces that move the hydrodynamics inside the bay, the depths average velocities for the 3D runs shows that Delft3D and OpenDA calibration in 2D was enough to explain the variability in currents along the bay with a very high correlation during the period of time considered, and there was no need to calibrate any other parameters for the 3D model. In this context predictability capabilities increase over time due the increased speed of the model to capture the dynamics. At the bay entry we found the main direction is NNE during high tide and SSW during low tide with a higher magnitude in this direction. A highlight from this work could be that, for this area, automated calibration has not been applied before as well the forecast capabilities of the model in 3D mode. Implementing Open DA as a tool for calibration in a small scale problem is also a new approach with successful results. In this aspect, the next step will be to consider OpenDA as a tool to deliver real time forecasting capabilities in this region by introducing data assimilation techniques. However new computational challenges need to be addressed to reduce computational costs of these implementations. Knowing the development of these variables over time is very important in the study of growing oysters Crassotrea gigas inside the bay. The 3D output from the model plays a very important role, since it is possible to determine which shallow areas inside the bay has enough current to transport nutrients 7
y component of depth averaged velocity (m/s)
x component of depth averaged velocity (m/s)
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
8
Depth average Velocity EASTWARD (V ) SQ Bay 2010 at SOUTH Station x Observatios Vs Model Output Observed − South Station Model Calibrated Model Forecast
1
0.5
0
−0.5
−1 27−Aug
29−Aug
31−Aug
02−Sep
04−Sep
07−Sep 09−Sep time−−−−>
11−Sep
13−Sep
15−Sep
18−Sep
20−Sep
Depth average Velocity NORTHTWARD (V ) SQ Bay 2010 at SOUTH Station y Observatios Vs Model Output
Observed − South Station Model Calibrated Model Forecast
1
0.5
0
−0.5
−1 27−Aug
29−Aug
31−Aug
02−Sep
04−Sep
07−Sep 09−Sep time−−−−>
11−Sep
13−Sep
15−Sep
18−Sep
20−Sep
Figure 5: Comparison of instantaneous depth average velocity V x Eastward (Top Figure) and Vy Northward (Bottom Figure) at the Bay Entry. 14 days simulation from the calibrated model starting on Aug 27, 2010, and 10 days forecast after calibration is done (black line)
ORIGINAL Model Eastward Velocity (Ux) R=0.80648 RMSE= 0.058906
ORIGINAL Model Northward Velocity (Uy) R=0.81851 RMSE= 0.1751 0.5
Uy Model
Ux Model
0.5
0
−0.5 −0.2
0
−0.5 −0.15
−0.1
−0.05
0
Ux Observed
0.05
0.1
0.15
0.2
−0.8
−0.6
CALIBRATED Model Eastward Velocity (Ux) R=0.87777 RMSE= 0.10146
0
−0.5
0.2
0.4
0.6
0.8
0.6
0.8
0.6
0.8
−0.5 −0.15
−0.1
−0.05
0
Ux Observed
0.05
0.1
0.15
0.2
−0.8
−0.6
−0.4
−0.2
0
Uy Observed
0.2
0.4
FORECAST Model Northward Velocity (Uy) R=0.86441 RMSE= 0.21025 0.5
Uy Model
0.5
Ux Model
0
Uy Observed
0
FORECAST Model Eastward Velocity (Ux) R=0.88622 RMSE= 0.11834
0
−0.5 −0.2
−0.2
0.5
Uy Model
Ux Model
0.5
−0.2
−0.4
CALIBRATED Model Northward Velocity (Uy) R=0.76106 RMSE= 0.18155
0
−0.5 −0.15
−0.1
−0.05
0
Ux Observed
0.05
0.1
0.15
0.2
−0.8
−0.6
−0.4
−0.2
0
Uy Observed
0.2
0.4
Figure 6: Linear regression between observations and 3D Model, for vertical depth average of Vx and Vy, before calibration, calibrated model and forecast. 8
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
for healthy cultivation of the oysters. Additionally, it is important to note that the quality of the baroclinic structure strongly depends upon the quality of the boundary conditions for temperature and salinity in this type of scenario. A future work will study temperature variation along the bay to determine the best cultivation zones for these organisms.
9
[19] M. L. Raltson, R. I. Jennrich, DUD, A derivative-free algorithm for nonlinear least squares., Technometrics 20 (1) (1978) 7–14.
Acknowledgements Special thanks to the Navy Secretary of Mexico (SEMAR), for the use of the currents data at the entrance of the bay, as well as Dr. Rafael Blanco from UABC-Mexico and Jody Payne for their invaluable support during this research. References [1] E. S. Ghada, M. Verlaan, S. Hummel, A. Weerts, J. Dhondia, OpenDA open source generic data assimilation environment and its application in process models., Geophysical Research Abstracts 12: EGU2010-9346-2. URL http://adsabs.harvard.edu/abs/2010EGUGA..12.9346E [2] G. Y. El Sarafy, H. Gerritsen, S. Hummel, A. H. Weerts, A. E. Mynett, M. Tanaka, Application of data assimilation in portable operational forecasting systems the DATools assimilation environment, Ocean Dynamics 57 (4-5) (2007) 485–499. [3] A. H. Weerts, G. Y. El Serafy, S. Hummel, J. Dhondia, H. Gerritsen, Application of generic data assimilation tools (DATools) for flood forecasting purposes, Computers Geosciences 36 (4) (2010) 453–463. [4] A. Kurniawan, H. Gerritsen, S. K. Ooi, S. Hummel, Sensitivity analysis of the tidal representation in Singapore Regional Waters in a data assimilation environment, Ocean Dynamics 61 (2011) 112–1236. [5] F. Zijl, M. Verlaan, H. Gerritsen, Improved water-level forecasting for the northwest european shelf and north sea through direct modelling of tide, surge and non-linear interaction, Ocean Dynamics 63 (7) (2013) 823–847. doi:10.1007/s10236-013-0624-2. URL http://dx.doi.org/10.1007/s10236-013-0624-2 [6] P. Chavez, M. Acosta, Desarrollo gonadal de crassostrea gigas en Bahia San Quintin, Baja California, Mexico. , Ciencias Marinas 20 (2) (1995) 225–242. [7] T. Imai, S. Sakai, Study of breeding of Japanese oyster, Crassostrea gigas., Agric. Res. 12 (1961) 125–171. [8] J. D. Wang, A. F. Blumberg, H. L. Butler, ASC, P. Hamilton, Transport Prediction in Partially Stratified Tidal Waters, Journal of Hydraulic Engineering 116 (1990) 380–396. [9] I. Del Valle Lucero, H. R. Cabrera Muro, Aplicacion de un modelo numerico unidimensional a bahia San Quintin, B.C., Ciencias Marinas 7 (1) (1981) 1–15. [10] O. E. Delgado-Gonzalez, J. A. Jimenez, J. L. Ferman-Almada, G.-E. Z, F. Marvan-Gargollo, A. Mejia-Trejo, Depth and hydrodynamics as tools to select aquaculture areas in the M Depth and hydrodynamics as tools to select aquaculture areas in the coastal zone, Ciencias Marinas 36 (3) (2010) 249–265. [11] I. Ramirez, R. Blanco, R. Vaquez, G. Ramirez, The simulation of the circulation of San Quintin Bay., Submited to Ocean Dynamics. [12] J. Jimenez, A. Mejia-Trejo, Investigacion numerica de la dinamica de la circulacion en bahia de san quintin, b. c.: procesos de transporte, Vol. 12, XII Congreso Nacional de Oceanograf´ıa, 2000. [13] I. Ramirez, D. Vincelli, Mediciones oceanograficas en bahia san quintn, b.c, mexico., CICESE, TECH. [14] Kriging [online]. [15] J. I. Gonzalez. MAR V1.0 para Windows 98/2000/XP/Vista/Windows7 [online] (Sep. 2013). [16] H. Delft, User manual Delft3D-Flow, DELTARES (Sep. 2011). [17] G. D. Egbert, A. F. Bennett, ‘TOPEX/ POSEIDON tides estimated using a global inverse model.’, J. Geophys. Res (1994) 821–824. [18] RIKZ, User’s Guide TIDEMAT tidal analysis, ALKYON Hydraulic Consultancy and Research, PO Box 248 8300 AE, Emmeloord The Netherlands (Jan. 2004).
9
M. Garcia et al. / Journal of Computational and Applied Mathematics 00 (2014) 1–9
10
Velocities Measured Bay Entry 90
1
120
60
0.5
150
30
180
0
330
210
240
300 270
(a) during high tide
Velocities Model Bay Entry 90
1
120
60
0.5
150
30
180
0
330
210
240
300 270
(b) during low tide
Figure 7: Depth Average Velocity, Magnitude (m/s)
Figure 8: Classical compass for current, measured data (Top Figure) and output from 3D model (Bottom Figure) during forecast period
10