Numerical simulation of the Adriatic Sea principal tidal constituents

Abstract. The primary goal of this study was to incorporate data-derived harmonic constants into a complex dynamic model using a form of variational data assimilation, with a view to improve the prediction of 7 dominant tidal constituents in the Adriatic Sea. Firstly, harmonic-constant data for 6 Adriatic stations were fed into a steady-state, 3-D, forward/inverse model to furnish optimal boundary conditions (OBCs). Calculated OBCs were then used to derive individual constituent responses, as well as to synthesise seven-constituent boundary conditions for the time stepping, 3-D model. A separate set of 25 stations provided control harmonic constant data. In validating the model output particular attention has been given to the often-ignored tidal currents. To that end 14 current meter data records were processed into tidal current ellipse parameters and used to examine the comparable model output. Comparison to gauge data has shown that the present solution is better than our own previous one, and shows an improvement over recent solutions by other authors. The model accurately reproduces available data with individual station amplitude differences rarely exceeding 1cm, and with the phase error commonly staying well below 10°. For all tidal constituents individual station differences result in RMSE in the 0.33-0.71-cm range for amplitude, and the 5.6°-19.2° range for phase. Semidiurnal currents appear to be modelled better than the diurnal ones (generally over-predicted). High eccentricity of both data and model-derived ellipses often impaired calculating the proper sense of rotation; inclination of the ellipses proved to be the most robust parameter, successfully predicted for most constituents at all depths.


Introduction
Recent years have witnessed a renewed interest in the Adriatic tides. Janeković et al. (2003) used a data assimilation Correspondence to: I. Janeković (ivica@irb.hr) approach, combining a 2-D assimilative and 3-D forward model, to obtain dynamically consistent optimal boundary conditions for the two major tidal constituents (M2 and K1). The authors found that the data assimilation contributed to a significant reduction in the model-to-gauge data discrepancy; the mismatch was more pronounced at the shallow north, with M2 amplitudes generally under-predicted, and K1 overpredicted. Cushman-Roisin and Naimie (2002) applied a 3-D finite-element model to the Adriatic, simulating 4 major Adriatic constituents, separately and in concert, as well as their residual effect. Comparison to available coastal station elevations, in particular, demonstrated their model skill. In a study based on basin-wide, vessel-mounted ADCP and complementary, moored, current-meter measurements, Ursella and Gacić (2001) also addressed the Adriatic tidal structure. The two kinds of measurements were combined to address the tidal pattern changes in relation to stratification. Malacic et al. (2000) investigated the tides in the Adriatic Sea with a very high resolution 2-D model focused on the northern part, and a 3-D model of the whole basin. They reinforced the traditional interpretation of the semidiurnal constituents as an incoming Kelvin wave propagating along the eastern coast and then reflected at the basin's northern end. However, in their interpretation the diurnal constituents are topographic waves travelling across the basin, from the eastern to the western shore. Nevertheless, both species are seen as members of the same family of linear waves, existing under the combined influence of gravity and topography.
The studies mentioned build on a century long tradition, addressed in more detail in Janeković et al. (2003). From early studies (Defant, 1914;Sterneck, 1919) to later modelling works (e.g. Hendetshott and Speranza, 1971;Accerboni and Manca, 1973) the Adriatic tides have been primarily interpreted as co-oscillations with the Mediterranean Sea. Mediterranean modelling studies incorporating the Adriatic (e.g. Canceill et al., 1993;Tsimplis et al., 1995) seem to reinforce that view. It is worth mentioning here that the Adriatic Sea is among very few Mediterranean Sea basins which exhibit moderate tides (Mediterranean tidal signal is predominantly weak). This is particularly true for the shallow, confined shelf of the northern Adriatic, where observed M2 and K1 amplitudes approach 30 cm and 20 cm, respectively. The goal of the present paper has been to incorporate dataderived tide constituent constants into a complex dynamic model using a variant of incremental variational data assimilation. To reach it we have consolidated the data assimilative procedure (leading to more accurate predictions), improved the model's spatial resolution (along islands -torn eastern Adriatic coast, in particular), and extended the spectrum of the considered tidal constituents to the seven most pronounced ones. In validating the model output particular attention has been given to the often-ignored tidal currents, and the more demanding requirements that their reproduction incurs. To that end the data sets and processing are briefly presented in Sect. 2. The models and assimilation methodology are the subject of the Sect. 3, whereas results are presented and discussed in Sect. 4. The work is summarised in the final section.

Data
In this study we rely on published tidal constituent constants data, as well as on our own analysis of the time series of sea level data from 6 gauge stations and current meter data records from 9 mooring stations. Time series of available data, both at sea level and for currents, were analysed with the same set of Matlab routines (Pawlowicz et al., 2002), minimising the differences which may arise due to using different types of tidal analysis software. For shorter time series we have used inference from stations which have long enough time series, in order to resolve tidal harmonics K1/P1 and S2/K2; for details, see Pawlowicz et al. (2002) and Foreman (1977).
Another type of data is original tide-gauge time series for 6 stations: Bakar, Rovinj, Ancona, Split, Dubrovnik and Trieste (marked with large, bold numbers on Fig. 1 and bold names with an asterisk in Table 1). For the Bakar gauge station we analysed 1 year of hourly data (1 January 1985-31 December 1985, producing a set of 59 tidal constituents; for Rovinj station 2 years of hourly data (1 January 1992-31 December 1993), extracting 68 tidal constituents; for Ancona station 212 days (10 January 1982-17 August 1982), resolving 59 constituents; for Split station 2 years of hourly data (1 January 1986-31 December 1987), resolving 68 constituents; for Dubrovnik station 2 years of hourly data (1 January 1986-31 December 1987), resolving 68 constituents and for Trieste 5 years (1 January 1989-31 December 1994), determining 68 constituents. For those data sets we also computed the appropriate error estimates at a 95% confidence interval. Only seven major tidal constituents were further used in the assimilation part of the study, with the exception of the Trieste station, which was left out for independent comparison. The choice was made keeping in mind that tidal signal is the strongest in the shallow, northernmost part of the Adriatic, and that we have had at our disposal a reliable fiveyear long series. Closer inspection of the literature sources revealed notable discrepancies between the different sources (particularly for weaker tidal constituents). For example, we have found eight different published values for the most energetic tidal constituent (M2) in Trieste; the differences are illustrated in Fig. 2.

Currents
The first set of current data was obtained from 8 Aanderaa RCM-4 current meter stations (marked A to H on Fig. 1) with deployment details presented in Table 2. The data were recorded on digital tape and were made available for reanalysis, making a database of 13-time series, mostly measured in wintertime or early spring when barotropic conditions were present. Unfortunately, the data appear to be of uneven quality, affecting the results of tidal analysis (especially for phases, because some data in the files appear to be missing). Beside those data sets we used recently recorded ADCP current meter data (marked with letter J in Fig. 1 and Table 1. Tide-gauge stations with harmonic constant values used in the study; bold names with the asterisk mark the stations for which we calculated the constants. Amplitudes are in cm and phases are degrees relative to 15 E. Table 2), programmed to sample the water column with a cell size of 2 m and a time interval of 10 min, with 50 pings per ensemble. That data were analyzed in the same manner as the RCM-4 data, but providing a much higher vertical resolution (19 useful vertical layers). It is important to point out that the obtained results for currents were used only to verify the flow field obtained from model simulations, not in the data assimilation procedure.

Models and assimilation technique
Although the Adriatic Sea has only one short open boundary across the Straits of Otranto, it is not a trivial one, as there has been no published pressure gauge data covering it. Attempts to estimate the amplitude and the phase are limited to the available coastal stations, and some assumed spatial distribution in between. That method can result in reasonable estimates of tidal variability within the basin but it does not always provide a dynamically balanced solution for the open boundary. Recent work by Janeković et al. (2003) showed that by using incremental data assimilation (DA) and a combination of 2-D and 3-D models, one can successfully obtain the OBCs for the 3-D model. In that approach the 2-D model with its adjoint was used in the assimilation part, and the 3-D model was used as the forward counterpart. In this study we decided to use another type of DA based on the use of a frequency-domain, 3-D, finite element, linear model and its inverse, to obtain OBCs for each tidal constituent separately. The method is similar to the one described in Janeković et al. (2003); its central point is the use of the two models iteratively. An assimilative part (3-D linear model with its inverse) provides OBCs a for 3-D, prognostic, nonlinear model which generate error estimates at assimilated stations. These error estimates are then used to improve the assimilative part and provide updated OBSs for the forward model. This methodology was successfully applied to Georges Bank by assimilating ADCP observations (Lynch et al., 1998). Both (inverse and forward) assimilative models use the same mesh resolution and bathymetry as the main 3-D, nonlinear, time stepping, finite element model. Boundary conditions obtained that way are taken to be optimal and were combined together, to synthesize the sea level used in forcing the 3-D, nonlinear model.

Three-dimensional data assimilative model
The model used in the DA part of this study is the finite element model of Lynch et al. (1993) ("Truxton/Fundy"). This model solves the 3-D, linearised, hydrodynamic equations for harmonic-in-time motions. The model is formulated in terms of the complex amplitudes Z(x, y) and V (x, y, z) of the sea level elevation ξ and velocity v, at the given frequency ω for each of the seven tidal constituents. The time-domain response is then: Table 2. Current meter stations from which data were used in this study with deployment details.
For the Dirichlet boundary value problem described here, the domain response V and Z is linear in the complex elevations Z bn at the open boundary nodes bn=1 . . . N , for a given frequency ω where [A] is matrix of the domain system response (representation of model equations), {Z bn } is the vector of the OBCs and {Z} is the modelled complex sea level amplitude due to the response of the system to boundary forcing. The goal is to best fit the model to observations in the leastsquares sense, so we define the error vector as a mismatch between modelled and observed complex amplitude values, as: There are many strategies to deal with solving this type of equation, in order to obtain the vector of boundary conditions {Z bn } for which error {ε} has a minimum value. Depending on the condition of the system matrix [A] (singular or not), one can use single value decomposition (SVD) or weighted least-squares (WLS). The latter method uses an additional weighting function on the solution vector; restrictions on the size (w 0 term) and the slope (w 1 term) of the obtained open boundary elevation, as well on the error vector (Lynch et al., 1998): with: Integration is made at the open boundary only (ds is segment), where solution vector {Z bn } is defined. We have used the WLS method for observations of sea level (tidal constants) at 5 tide gauge stations (Bakar, Rovinj, Ancona, Split and Dubrovnik) where we were able to analyse and obtain harmonic constants for seven major tidal constituents. In order to obtain satisfactory results, we had to add to the previous five one more gauge station  in Table 1), due to otherwise poor control of the system. In the DA model we used the value of k=0.001 m/s for the coefficient of the linearized bottom stress and constant vertical viscosity N=0.04 m 2 /s, based on our experience, which gives the best model-to-data fit. A similar type of inversion (SVD) and the Truxton/Fundy model has been used and in detail explored by Xu et al. (2001).
A finite element grid used in the study consists of 23 055 nodes and 37 200 elements (Fig. 3). The size of the triangle areas varies from 0.02 km 2 to 757 km 2 , with the length of the nodal distances varies from about 500 m in coastal areas, up to 44 km in the deep parts of the domain. With this mesh we have been able to include in the simulations 77 major islands and recognise realistic topography and lateral geometry better than any previous Adriatic tidal model. Representing local features (narrow channels, many islands) correctly is an important requirement in realistic tidal modelling. To illustrate the capability of the new grid and the richness of Croatian coastal features, we show a magnification of the mesh for Kvarner Bay in the upper right part of Fig. 3.

Three-dimensional forward model
The model is of a finite element variety ("Quoddy"), based on the 3-D, nonlinear, shallow water equations (Lynch et al., 1996). The 2.5-level turbulence-closure scheme of Mellor and Yamada (1982) is used in the model with the improvements of Galperin et al. (1988). The horizontal diffusion parameterisation scheme is that of Smagorinsky (1963). The free-slip condition is imposed along the coast. The details of the model solution on a finite element grid are described in Lynch and Werner (1991). Bottom stress is estimated from the classical quadratic law as a function of bottom velocity. We use a value of 0.003 for the drag coefficient, based on a set of numerical experiments targeted to produce the most appropriate value. The well-known "sigma" coordinate system is used in the vertical, with 21 non-uniformly placed nodes whose sinusoidal vertical spacing provides an increased resolution in the surface and bottom layers. This high resolution grid is necessary for adequate simulations of the turbulent dissipation and details of the currents, especially along the complex Croatian coast. The model is forced only by time varying, sea level boundary conditions along 40 • N parallel, synthesised from the output of all seven individual constituents generated by the DA model. The same high resolution, finite element mesh is used in the forward model as in the DA model.

Results and discussion
It has been shown in a previous study (Janeković et al., 2003) that the direct astronomical forcing is not very important for the Adriatic tides. The direct astronomical influence can account, at most, for about 1% of the observed M2 and about 6.2% of the K1 amplitude, leaving co-oscillations of the Adriatic and Ionian Sea through their common boundary as the main driving force. The above-described DA procedure was used to derive the sea level at the open boundary.

Optimal boundary conditions
Using the DA and 6 tide-gauge data sets (five for which we had longer hourly time series, and one with only the published constants, but close to the open boundary), we have obtained the sea level OBCs for each of the seven tidal constituents (Fig. 4). Assimilating values from the shallower northern part of the Adriatic Sea (tidal amplitudes have higher values there, for both semidiurnal and diurnal tidal constituents), we have been able to infer the values at the open boundary for all tidal constituents, including relatively weak ones. Obtained open boundary amplitude structure exhibits relatively small values, with only 3 tidal constituents (M2, S2 and K1) surpassing the 2-cm mark (Fig. 4, top). The phase structure suggests two groups: diurnal with angle values between 50 • and 65 • , and semidiurnal with phases ranging between 95 • and 115 • , for particular tidal constituents (Fig. 4, bottom). The OBCs obtained for the two most energetic tidal constituents (M2 and K1) are in good agreement with the values obtained in our previous study (Janeković et al., 2003). Small differences may originate from the different methods employed in the previous DA process (incremental, with a combination of 2-D finite differences and 3-D finite element modelling). Presently derived OBCs were used to run single constituent cases, as well as for a combined, transient seven-constituent case.

Sea level response
The Adriatic basin response to individual constituent forcing has been studied before (including by the present authors - Janeković et al., 2003). In this study we have focused on the physically, more realistic, combined response of the seven major constituents. To maintain correspondence with previous works we also performed some single-constituent runs. Co-tidal charts for the S2, N2 and K2 (not shown) tidal constituents generally repeat the pattern of the dominant M2 tidal constituent shown in Fig. 5. All semidiurnal tidal constituents picture a cyclonically rotating amphidromic system, appearing in earlier studies, like that of Polli (1959), or Accerboni and Manca (1973), but with a somewhat different positioned amphidromic point. The Polli's solution is centred on the Sibenik-Ancona line (our stations 15 and 17), whereas Accerboni and Manca have it closer to the Italian coast. Hendershott and Speranza (1971) were apparently the first to interpret the system in terms of incoming and reflected Kelvin waves.
The O1 and P1 follow the pattern of the K1 tidal constituent shown in Fig. 6. These constituents do not display an amphidromic point in the Adriatic. Instead, the waves appear to travel in the northeast to southwest direction, with the amplitude rising northward along the central axis -behaviour interpreted by Malacic et al. (2000) in a topographic wave framework. In both groups the amplitudes reflect the weights that a particular tidal constituent carries. The relative importance of individual diurnal and semidiurnal tidal  constituents can be compactly expressed through the Form Number defined as: where H X is the respective amplitude value. Although not strictly valid in every particular moment, the number provides a measure of diurnal-semidiurnal interplay; the palette in Fig. 7 is designed to reflect traditional F categories. One can observe that there is neither a diurnal nor a semidiurnal part of the Adriatic. The response is mixed, mainly semidiurnal, with a circular area of primarily diurnal response centred at the amphidromic point. The model solution confirms, with dynamic consistence and spatial detail, predominantly qualitative results obtained previously and reported (e.g. Defant, 1960). In order to obtain more a realistic seven tidal constituents response, the time-dependent model was forced in time for three months with free surface OBCs (synthesised from the single constituent results). After a one-month spin-up time, two more months of hourly values were calculated and stored at selected locations, and later harmonically analysed. Due to the relatively short modelled time series of the sea level, we used the previously described inference technique to resolve close (in frequency domain) tidal constituents (for details, see Data Section). The success of this labour is summarised in Table 3, where model-to-data differences are tabulated for 7 constituents at 31 locations. Comparison of modelled amplitudes with those data-derived at gauge-station locations generated small basin-wide differences. The model accurately reproduces available data with amplitude differences rarely exceeding 1 cm (e.g. at Venezia-Lido: Station 4, or in the Zadar Canal: Stations 11 and 12), but the phase errors diverge from observations considerably; the difference for N2, for example, approaches 60 • at Svetac (Station 21). The individual station differences result in RMS constituent errors in the 0.33-0.71 cm range for the amplitude, and 5.6 • -19.2 • for the phase. One should note, however, that the errors are smaller than those we reported previously for M2 and K1, or those that Cushman-Roisin and Naimie (2002) reported for M2, S2, K1, and O1.
The Zadar Canal is worth a further comment. It is a long, narrow, and shallow structure, horizontally resolved in the model. However, although finer model resolution helped to improve the semidiurnal phase delays (along-channel travelling waves), the prediction still falls short of closely reproducing the observed phase values (Fig. 8). Although still less than perfect, the solution clearly illustrates an improvement possible with the use of high resolution, unstructured mesh.
One should also bear in mind that the uncertainty in the data alone could be rather high. For example, analysing time series of five-year, hourly, sea level data from Trieste (St. 1), we obtained the phase error of 4.6, 5.1, 1.4, 1.9, 0.3, 0.5 and 1.4 deg for the O1, P1, K1, N2, M2, S2, K2 constituents, respectively. For the 200-day long series collected at Ancona (St. 17) the result is considerably worse: 14.4, 12.8, 4.9, 7.8, 1.5, 2.9 and 11.8 deg for the O1, P1, K1, N2, M2, S2, K2 tidal constituents, respectively. Nevertheless, our solutions, qualitatively similar to older and more recent ones -quoted earlier in this section -are quantitatively closer to, and better supported by, the empirical data. Observed discrepancies may be attributed to various procedural inadequacies, ignoring problems like meteorological influence on empirical derivation of harmonic constants. Crisciani et al. (1995), for example, suggest that part of the harmonic constant's variability they calculated can be attributed to a meteorological regime at the time of the observations (atmospheric pressure and wind).
A different view of the same differences is offered by inspecting Fig. 9. One can readily observe a lineup of the semidiurnal amplitudes along the equality line (Fig. 9A), and pronounced scatter of diurnal phases (Fig. 9D). Diurnal amplitudes (Fig. 9B) are also relatively lined up, with the comparable model overshoots and undershoots of K1.
The phases of less pronounced semidiurnal tidal constituents (K2 and N2) have proved to be more difficult to predict than the more energetic ones (M2, S2, Fig. 9C). As a further test of the model skill we addressed the question of how well does a numerical seven tidal constituent solution reproduce observed sea level data at a particular gauge station. To that end we selected the Trieste station (No. 1), not used in the assimilation procedure, yet providing long time series of high quality. The 5-year series of hourly data was analysed, and from the 68-harmonic output the major seven constituents were extracted. The same program (see Data Section) was used to analyse the two-month long series of hourly values derived from the model at the Trieste location, and to calculate tidal constants for the same dominant constituents. Again, due to a relatively short model time series and the Rayleigh criteria used the inference method was used to resolve S2/K2 and K1/P1 tidal constituents. The constituent's sets were then used to synthesise signals for the first two weeks of 1981. The outcome is presented in Fig. 10. One can readily observe not only the correctly reproduced mixed, predominantly semidiurnal, nature of the response at Trieste, but also an excellent agreement in numerical details.

Current response
The sea level alone is not sufficient to judge a model's skill in reproducing tidal dynamics in a particular basin. Tidal currents provide important additional information, particularly when a fully three-dimensional structure is provided. But current vectors have a higher level of irreducible variability than the scalar field of the sea level (Godin, 1983), Table 3. Amplitude and phase differences between observation and model results for sea level at 31 gauge stations. Differences in amplitudes are in cm and phases are in degrees. The last three lines represent the basic statistical parameters: standard deviation (STD), root mean square error (RMS) and absolute mean differences (|M|). Table 3. Amplitude and phase differences between observation and model results for sea level at 31 gauge stations. Differences in amplitudes are in cm and in phases in degrees. The last three lines represent basic statistical parameters: standard deviation (STD), root mean square error (RMSE) and absolute mean differences (|M|).   and generally exhibit lower horizontal and vertical coherence. Consequently, the current records are used less often than the sea level in tidal model validation studies. In our previous study (Janeković et al., 2003) we used 6 current meter records from 2 mooring locations along the West Istria coast to verify our M2 and K1 predictions. An excellent agreement in orientation was reported, with maximum error in the semi-major and semi-minor axes of 1.8 cm and 1.5 cm, respectively. Cushman-Roisin and Naimie (2002) used previously reported Italian data for the Gulf of Trieste (M2, S2, K1, O1) and Emilia-Romagna coastal area (M2), to assess the skill of their model to predict the vertically averaged currents, and generally found good agreement in both amplitude and direction. In the present study we have found that although the four semidiurnal and three diurnal tidal constituents we simulated have different current magnitudes and phases, they exhibit similar intra-group behaviour, patterned after the M2 and K1 responses, respectively. Vertically averaged M2 currents (not shown) are most pronounced in the Northern Adriatic, less articulated in the Middle Adriatic, and barely traceable in most of the Southern Adriatic. The vertically averaged ellipses in the coastal stretch along the  West Istria coast, in particular almost collapse to a straight segment, exhibiting the largest major axes. A similar broad description is valid for the K1, with the additional effect of its cross-basin amplification roughly tracing the 200-m isobath. Generally, the vertically averaged current response reinforces the one reported by Cushman-Roisin and Naimie (2002).  The problem of the three-dimensional response is much more demanding than the vertically averaged flow, invoking questions of vertical shear and turbulence parameterisation. To further aggravate the problem the density of current data stations and continuity in data collection are still inferior to the respective gauge data achievements. In an effort to address our model skill in predicting the three-dimensional current response, we looked at 14 current observation records collected at 9 stations in different field programs. Eight Aanderaa stations provide 13 records, and one ADCP gives 19 depth cell records.
Observations were made in winter and early spring when the water column was well mixed, in keeping with the barotropic conditions assumed in the model. The mooring deployments of these records span different time periods, as summarised in the Data Section. As opposed to tide-gauge Table 4. Tidal-ellipse parameters for the seven major tidal constituents from observed (obs) and modelled (mod) data at the station J and two depths. Inclination angles are in degrees and semi-axes in cm/s. Table 4. Tidal-ellipse parameters for the seven major tidal constituents from observed (obs) and modelled (mod) data at the station J and two depths. Inclination angles are in degrees and semi-axes in cm/s. records the processed current observations were typically limited to a month or two, not allowing for a separation of major tidal constituents (e.g. K1 and P1) without using the inference method. The currents were not assimilated in any of the simulations. Examination of tidal ellipses for both tidal constituents reveals a decrease in the current magnitude with depth, with M2 ellipses more polarised at the surface. Both constituents exhibit excellent model-to-data agreement in ellipse orientation at all stations and at both depths. The magnitude of the major axes, the M2 in particular, is also well predicted. The minor axis is more of a problem, often because of the smallness of its magnitude, which can be appreciated better by examining numerical values. Table 4 lists the major and minor axes magnitudes, together with directional orientation of the major axis and the S/N ratio. The error bounds are also provided for the first three parameters. The same information is reported for model-and data-derived quantities. To keep the table down to a manageable size, the above is reported for the seven constituents, but at only one station (J). The table reveals several notable points. One readily observes narrower model error margins, which can be attributed to a much higher model S/N ratio. The prominence of the three tidal constituents (M2, S2 and K1) in terms of the observed (and modelled) semi-major axis is clearly visible. One also notes a larger observational error for the K1 constituent (0.9 cm/s), as well as its over-prediction (5.9/4.3 cm). For the semidiurnal one the error is smaller (0.5 cm/s). All semi-minor axes are rather small with large observational errors which make the results suspect. All tidal constituents have a calculated inclination of the major axis within, or very close to, the observational error of the measured mean. However, for the three dominant constituents the discrepancy is within 2 • , although the observational error for the K1 goes as high as 13.5 • . For the smaller diurnal ones (O1 and P1) it reaches an excessive 41 • , while the actual P1 discrepancy appears to be less than 2 • . The relative ranking of the constituents observed in the near-surface holds for the near-bottom, too, but smaller magnitudes of all semi-major axes are accompanied with observational error magnitudes, as observed near the surface. The orientation of the major axes continues to be a better-predicted parameter near the bottom, with the three dominant tidal constituents exhibiting discrepancies within the 2 • limit. The model's success in predicting the current ellipses of the other less pronounced constituents at this station (J) can be gleaned from Fig. 12.
A still finer zoom on vertical variability is given in Fig. 13, where the ADCP vertical resolution of the M2 and K1 dataderived and modelled parameters (major semi-axis, minor semi-axis, inclination) is provided. Clearly, the model is excellent in predicting the M2 major semi-axis and its variability with depth. A model-derived slight increase in maximum magnitude, followed by a sharp decrease, closely matches the observed behaviour. The modelled minor semi-axis stays close to, but only in the middle of the water column, when taking into account the error bounds. Depth veering of the M2 inclination (Fig. 13C) is correctly modelled in terms of direction and the rate of change, but the model inclination is somewhat outside the error bounds, except in the bottom layer. The K1 constituent pictures a somewhat opposite situation. The minor semi-axis (Fig. 13B) is reproduced remarkably well, although the observational error bounds are larger than the axis value. The predicted major semi-axis appears offset by about 2 cm/s throughout the water column, but the near-bottom decrease is correctly predicted. The dataderived inclination of the K1 ellipse (Fig. 13D) has a large uncertainty (about 10 • throughout the water column), but the clockwise turning of the major axis (by about 12 • ) is correctly reproduced. Where do the differences come from? We believe that the high similarity of the model-and data-derived responses near the bottom testifies to the proper modelling of the bottom friction. Very high spatial resolution of the model grid allows for adequate representation of lateral geometry and bottom topography, although sporadic, local bathymetric features may have remained unresolved, affecting the modelled topographic steering. We suspect that imperfectlyrepresented vertical mixing of momentum by the level 2.5 turbulence closure scheme is largely responsible for the mismatch.

Summary and conclusions
We have used a variant of incremental variational DA to incorporate data-derived tidal constituent constants into a complex dynamic model, with a view to provide an improved prediction of seven dominant tidal constituents in the Adriatic Sea. More specifically, tidal constant data for 6 Adriatic stations were fed into a steady-state, 3-D, forward/inverse model, in order to furnish optimal boundary conditions for the dominant-constituent simulations. Calculated OBCs were used to derive individual tidal constituent responses, as well as to synthesise the seven-constituent boundary condition for the time-stepping, 3-D model. A separate set of 25 stations provided control constituent constant data. As a further test of the time-stepping model skill, 14 current meter data records were processed into tidal current ellipse parameters and used to examine the comparable model output. No current data were used in the assimilation procedure.
Comparison of modelled amplitudes with those dataderived at gauge-station locations generated small basinwide differences. It suggests that the present solution is better than our own previous one (Janeković et al., 2003 -M2 and K1 only), as well as that of a recent contender (Cushman-Roisin and Naimie, 2002 -M2, S2, K1, O1). In the latter case the RMS error in S2 amplitude is the only exception: we obtained 0.5 cm, compared to their 0.4 cm -a discrepancy for which we do not see an obvious explanation. For the remaining 3 tidal constituents, not modelled previously (P1, N2, K2), the model-to-data comparison singled out the phase as the source of the occasional significant RMS error. However, a lack of conclusive processing details about published tidal constituent constants often leaves one suspecting uncertainty in the data, as much as the model inadequacy.
Comparison of the current response has offered more encouraging results, although observational/processing uncertainty plagues those results, too. The signal-to-noise ratio of the observed and modelled current series offers a useful guidance. Semidiurnal currents appear to be modelled better than the diurnal ones (generally over-predicted). High eccentricity of both data-and model-derived ellipses often impair calculating the proper sense of rotation, but the inclination of the ellipses proved to be the most robust parameter, successfully predicted for most constituents at all depths.
A major roadblock to further improvements in the barotropic framework appears to be the quality and quantity of current data records. Major international field programmes recently completed in the Adriatic hold a significant promise in that respect. We hope that our immediate future work will be a contribution to that end, a necessary step before the more demanding baroclinic Adriatic tidal response is tackled.