Vertical structure of bottom Ekman tidal flows: Observations, Vertical structure of bottom Ekman tidal flows: Observations, theory, and modeling from the northern Adriatic theory, and modeling from the northern Adriatic

[ 1 ] From September 2002 to May 2003, fifteen bottom-mounted, acoustic Doppler current profilers measured currents of the northern Adriatic basin. Tidal fluctuations at all seven of the major Adriatic frequencies were synthesized from a response tidal analysis of these measurements. Most observed tidal current ellipses were nearly reversing, but near the bottom, tidal current ellipses all shortened and broadened, semidiurnal currents led upper water column currents, and diurnal tidal current ellipse orientations rotated counterclockwise toward the bottom. Theoretical solutions for a tidally forced, bottom Ekman layer with vertical eddy viscosity of the form A z = b z + k were least squares fit to the observations. Average values were b = 3 (cid:1) 10 (cid:2) 4 m/s and k = 5 (cid:1) 10 (cid:2) 4 m 2 /s. The value of k was important in matching tidal orientation and phase changes, and a nonzero b was important in matching tidal amplitude changes. The Navy Coastal Ocean Model (NCOM) and the Quoddy model were also compared to the observations. The average RMS errors for the bottom Ekman layer were 0.22 cm/s for the best fit theory, 0.35 cm/s for NCOM, and 0.36 cm/s for Quoddy. A z structures from NCOM and Quoddy show that time variation in A z is relatively unimportant for Adriatic tides. The bottom shear stresses from theory were larger in magnitude than those from the bottom drag formulations in NCOM and Quoddy.


Introduction
[2] In his pioneering paper on the dynamics of tides on the north Siberian Shelf, Sverdrup [1927] first solved the bottom Ekman theory problem with tidal forcing and found agreement with observed tidal structure changes near the bottom.Since that time, the topic has been further explored and advanced by many investigators.A nonexhaustive list includes the mainly theoretical works of Fjeldstad [1929], Prandle [1982], Soulsby [1983], and Yasuda [1987], and the works of Munk et al. [1970], Kundu et al. [1981], Maas and van Haren [1987], Lueck and Lu [1997], Ullman and Wilson [1998], Tsimplis [2000], Werner et al. [2003a], and Davies et al. [2004], all of which compared theories of vertical current variation to current measurements.
[3] However, despite this body of work, quantitative evaluation of the theoretical vertical changes in tidal structure has been limited by lack of appropriate measurements.Both the studies of Kundu et al. [1981] and Maas and van Haren [1987] verified the major characteristics of tidal Ekman theory, but certain aspects of their fits were less than satisfactory (e.g., the veering of ellipses in the study by Kundu et al. [1981] and the diurnal constituent comparisons in the study by Maas and van Haren [1987]) and they faced major measurement limitations (i.e., vertical current measurements of limited duration from anchored ships by Kundu et al. [1981] and point current samples at limited depths by Maas and van Haren [1987]).Both Lueck and Lu [1997] and Ullman and Wilson [1998] found good agreement between logarithmic layer theory and measured velocity profiles over short durations but conducted individual fits in time with bottom drag coefficients varying over tidal cycles at sites where Coriolis accelerations could be neglected rather than comparisons to tidal Ekman theory.Recently, acoustic Doppler current profiler (ADCP) measurements have been compared to tidal Ekman theory by Tsimplis [2000] and Davies et al. [2004], but neither of these studies had measurements within 10 m of the seabed and therefore missed a region of strong change in tidal characteristics.Of the structure that was observed, Tsimplis [2000] explained 90% of the variance of along-strait velocity and phase using a combination of frictional theory and internal mode theory to account for internal tides.The 3-D model with a quadratic friction law used by Davies et al. [2004] agreed well with available measurements for semidiurnal tides, but less well for diurnal tides, likely because of inaccuracies in the measurements.Werner et al. [2003a] compared tidal velocity observations at one location to 1-D models using either a combination of linear and constant eddy viscosities or using a Mellor-Yamada level 2.5 turbulence closure scheme [Mellor and Yamada, 1982].They found good agreement between the data and both types of models during times of low stratification for M 2 tides, the dominate constituent of the region and the only one they utilized.
[4] Recent mooring measurements from the northern Adriatic present a new opportunity to evaluate tidal variation in the bottom boundary layer.Velocities were measured by ADCPs at 15 different locations spread throughout the basin for more than six months over the winter [Book et al., 2007b] together with pressure measurements from wave/tide gauges (WTGs) at most sites.Depth cell sizes were 1 m or smaller, measurements covered the entire water column with the exception of the bottom blanking and surface contamination zones, and measurements at all sites were made within 3 m of the bottom or closer.In addition, results from a 3-D finiteelement model of the Adriatic dedicated to tides [Janekovic ánd Kuzmic ´, 2005] and results from a 3-D finite-difference model of the Adriatic with tides [Martin et al., 2006] are available for comparison with these measurements and with theory.
[5] The Adriatic Sea is an arm of the Mediterranean Sea.It may be represented roughly as an 800-km-long, 150-kmwide channel, oriented southeast-northwest, open at the southeast end (Strait of Otranto), with the bottom sloping upward toward the closed northwest end.The northern Adriatic (defined here to occupy the region northwest of Ancona and Zadar) is the final 200 km of this ''channel,'' where depths gradually slope from 70 m in the southeast to less than 10 m in the northwest (Figure 1).The southwest (Italian) side of the sea is characterized by a mild bathymetry slope, the Po River Delta and associated Po River plume, and a boundary current, the Western Adriatic Current, with a typical strength of 10 cm/s.In contrast, the northeast (Croatian) side of the sea is characterized by nearly vertical dropoffs at the coast, numerous deep bays and channels, and diffuse and varying currents.Most of these bays and channels on the northeast side are nearly isolated from the main Adriatic, but Kvarner Bay opens up to the Adriatic through a 30-km-wide passage.Figure 1 shows some of the main features of the northern Adriatic.
[6] There has been considerable theoretical and practical research on the tides of the Adriatic (see Cushman-Roisin et al. [2001, chapter 7] for a review of work prior to 2001), but, as in many coastal areas, direct measurements of tidal currents have been limited by technological and fishing pressure restrictions.Malac ˇic ˇet al. [2000] extended earlier semidiurnal applications of the theory of Taylor [1921] by Hendershott and Speranza [1971] and Mosetti [1986] with a general theory of gravity and topographic waves to explain the dynamics of both the semidiurnal and diurnal Figure 1.Bathymetry of the north Adriatic and tide ellipses from vertically averaged currents.M 2 ellipses are drawn in magenta, K 1 ellipses are drawn in red, S 2 ellipses are drawn in yellow, O 1 ellipses are drawn in blue, P 1 ellipses are drawn in green, N 2 ellipses are drawn in cyan, and K 2 ellipses are drawn in black.The velocity scale is given in the bottom left corner.Place names used in this paper are labeled, with cities indicated by black dots.Moorings are located at the centers of the ellipses and named according to the label given to each ellipse set.tides of the northern Adriatic.Recent work [Cushman-Roisin and Naimie, 2002;Janekovic ´et al., 2003; Janekovic ánd M. Kuzmic ´, 2005;Martin et al., 2006] has focused on using 3-D, high-resolution, numerical models with realistic topography and Mellor-Yamada turbulence closure schemes to simulate the Adriatic tides, and then to validate the components of these simulations with available observations.In addition to their focus on depth-averaged tides, Janekovic ánd Kuzmic ´[2005] also compared simulated tidal vertical structure with ADCP observations at one station and found good agreement.Differences were attributed to imperfectly represented vertical mixing.The interaction of stratification and tides in the Adriatic, with implications for the vertical structure of tidal currents, is also presently a topic of investigation, with Chavanne et al. [2007] citing Po River stratification as a possible explanation for model/data discrepancies near the Italian coast and Mihanovic ´et al. [2006] publishing a dedicated study of internal tides in the Adriatic.
[7] This paper reexamines time-dependent, tidally forced, bottom Ekman structure using the data from 15 ADCPs and results from the Navy Coastal Ocean Model (NCOM) and the Quoddy ocean model.Section 2 describes the measurements, section 3 describes the models, section 4 describes the bottom Ekman tidal theory used in this paper, and section 5 compares vertical tide structures from observations, theory, and models.Finally, discussions and conclusions are presented in sections 6 and 7.

Measurements
[8] From September 2002 to May 2003, an array of RD Instruments (RDI) Workhorse Sentinel broadband ADCPs was deployed in the northern Adriatic as part of a Joint Research Project (JRP) between the U.S. Naval Research Laboratory (NRL) and the NATO Undersea Research Centre (NURC).The JRP moorings consisted of 14 trawl-resistant bottom-mounted ADCPs [Perkins et al., 2000] distributed along portions of 4 mooring sections.An additional upward looking ADCP was mounted near the base of a meteorological tower as described by Cavaleri [2000].These mooring positions are shown in Figure 1 and listed in Table 1 with their mean sea level depths.The full mooring sections were populated by both the JRP moorings and moorings from several international partners collaborating on the study of the northern Adriatic [Lee et al., 2005].In addition to the ADCP measurements of currents throughout the water column, bottom pressure (by ADCP or wave/tide gauge) was also measured at each site.Book et al. [2007a] provide further details of the mooring instrumentation and Book et al. [2007b] show monthly mean and storm-driven currents observed at these sites.

Tidal Analysis
[9] The JRP ADCPs were set to measure the currents using 62 -90 s bursts of 1 Hz pings every 15 min, except for the first half of the VR1 deployment, which used 16 min bursts of 0.5 Hz pings every hour.Quality control steps to exclude bad data included an objectively determined velocity error cutoff (velocity errors estimated from independent measures of vertical velocity), exclusion of ensembles with more than 40% (20% for surface measurements) of the data marked bad by internal RDI checks, and additional tests described by Book et al. [2007a].The surface echo interference zone was truncated above a time-varying level determined from a time series of sea surface height constructed from the pressure and acoustic-backscatter-intensity measurements.Linear compass drifts (less than 4°) in some records were verified to be false trends by tidal analysis and corrected by small, gradual rotations of current vectors.Despite the lack of any physical evidence of instrument malfunction, the orientation disagreement between the observed tidal ellipses at station VR5 with neighboring observed tidal ellipses and three independent modeling simulations strongly suggests a compass error is present in the VR5 data.Therefore, the currents at site VR5 were rotated 28°clockwise to align with modeled strong-constraint, variational, data assimilation predictions.
[10] Tidal analyses of the ADCP data were done individually for all depth cells using the Response Method of Munk and Cartwright [1966] on the 15-min, current, ensemble time series (time values assigned to the center time of the measurement bursts).Gaps were introduced into the ADCP records by the quality control steps summarized in the previous paragraph, so a method was used for 2-D interpolation and extrapolation to replace missing values in the ADCP records.Sensitivity tests show that the near-surface tidal solutions (primarily only the ADCP depth cell nearest the surface) were influenced by interpolation and extrapolation method choices, but solutions at all other depths were very insensitive to this as they had few gaps.
[11] After the missing values were filled in, the Response Method was used to solve for estimates of pure tidal time series and for approximations to the harmonic tidal coefficients for the O 1 , P 1 , K 1 , N 2 , M 2 , S 2 , and K 2 constituents.Tides were not calculated from surface ADCP depth cells with 50% or more data marked bad by the quality control procedures.Slightly different procedures, described by Book et al. [2007a], were used to calculate tidal currents at site VR1 to account for the change in ADCP settings midway through the deployment.
[12] For each depth level, error estimates were obtained by creating an ensemble of 40 normally distributed, random noise time series, each having variance equal to the ADCP measurement error variance.Then, the velocity time series of the tide and of the tide residual were obtained by response analysis.The residual was 2-h low-pass filtered and then added to the tide and to one of the random time series.The procedure was repeated for each random time series in the ensemble.95% confidence limits for the tidal ellipse parameters were assigned to 1.96 times the value of the RMS errors of the parameters from the ensemble.The RMS errors were all smaller than the parametric bootstrap error estimates for a standard harmonic analysis of the data using a colored bivariate noise model [Pawlowicz et al., 2002].

Tidal Observations
[13] Figure 1 shows the tidal ellipses of the vertically averaged currents.The eccentricity of the ellipses is high (i.e., nearly reversing tidal currents) at all sites except for KB1 in Kvarner Bay and VR1 and VR2 in the northwest corner.Except for KB1, the ellipse major axes are approximately aligned with the Adriatic axis.There is an increase in tidal currents at the Istrian coast; site VR6 had the strongest tidal current with an M 2 semi -major axis of 10 cm/s.M 2 tidal currents are strongest at all sites (average semi-major axis 7 cm/s), but S 2 and K 1 tidal currents also play prominent roles (average semi-major axes 4 cm/s and 3 cm/s, respectively).O 1 , P 1 , N 2 , and K 2 are all much weaker, with average semi -major axes of 1 cm/s.
[14] The vertical structures of the tidal currents for M 2 and K 1 are shown in Figures 2 and 3, respectively.For graphical clarity, the vertical structures of the tidal currents for sites KB1 and VR1 -2 are not displayed together with the other sites because of their different character (i.e., eccentricity, strength, and orientation).A gradual change in tidal characteristics is resolved in the ADCP measurements from depth cells about 20 m above the bottom to the deepest measured depth cell.For M 2 (Figure 2), the semi -major axes decrease toward the bottom, indicating weakening of the tidal current speeds, the semi -minor axes increase toward the bottom, indicating broadening of the tidal current ellipses and counterclockwise rotation of the current vectors around these ellipses, and the phases decrease toward the bottom, indicating upper currents lagging bottom currents.Changes in the ellipse orientation are small, with some clockwise rotation of the ellipse toward the bottom.This structure agrees with predictions from tidal-forced, bottom Ekman theory [e.g., see Soulsby, 1983], as will be further detailed in section 4 of this paper.The vertical structures of S 2 currents (not shown) have the same character as those of M 2 .
[15] The K 1 tidal currents (Figure 3) also have gradually decreasing semi -major axes and gradually increasing semi-minor axes toward the bottom.However, although the ellipse shape changes are similar to those for the semidiurnal tides, the ellipse orientation and phase changes differ.The K 1 tidal ellipses rotate counterclockwise toward the bottom, with larger rotations than the clockwise rotations of the semidiurnal constituents.The K 1 tidal current phase changes with depth are much weaker than the semidiurnal phase changes.At most sites, the K 1 phases slightly increase toward the bottom (bottom currents lag upper currents).These differences in character between diurnal and semidiurnal tidal currents near the bottom were theoretically demonstrated by Kundu et al. [1981] and are due to the fact that the semidiurnal and diurnal frequencies are respectively faster and slower than the inertial frequency at latitudes greater than $30°.[16] Figure 4 shows SS6 tidal current ellipses as representatives of typical tidal current structure.At middepths, the diurnal and semidiurnal current ellipses align in direction and remain relatively unchanged over a large depth range.Approaching the bottom, all the ellipses shorten and broaden, the diurnal ellipses veer drastically counterclockwise together, and the semidiurnal ellipses veer slightly clockwise.Phase changes are difficult to see in this type of graphic, but semidiurnal currents strongly lead near the bottom, while diurnal currents weakly lag.
[17] Some departures from these general trends should be noted.The K 1 tidal structure is somewhat different at sites (not shown) KB1, VR1, and VR2.At site SS2, there is a K 1 phase drop of 14°from 15 to 23 m above the bottom.This vertical anomaly could possibly be caused by the Po River plume, following the hypothesis of Chavanne et al. [2007] that the plume caused an abrupt K 1 phase drop in surface tidal currents in a horizontal band all along the Italian coast as measured by their high-frequency radars.Site VR4 shows an increase in the K 1 phase toward the surface.

NCOM
[18] The Navy Coastal Ocean Model (NCOM) is a 3-D, finite-difference numerical model based on the primitive equations and the hydrostatic, Boussinesq, and incompressible assumptions.NCOM specifies vertical mixing in terms of vertical eddy coefficients that are calculated according to either the Mellor-Yamada 2.0 or 2.5 turbulence closure scheme.The NCOM runs for this study used the 2.0 scheme, which is described by Mellor and Yamada [1974].Implementation of this scheme in NCOM and details of the model are described by Martin [2000].
[19] NCOM was implemented for the entire Adriatic Sea on a 1.02-km horizontal grid.A difference from previous descriptions of NCOM is that the version of NCOM used here was modified to allow the use of generalized sigma coordinates, where the fractional sigma layer thickness can vary horizontally as well as vertically and sigma layers can be masked to land as the bottom shallows.The vertical grid used for this study was set up with the following properties: 40 layers in deep water, with a gradual reduction to 7 layers in shallow water, logarithmic expansion of the grid away from the surface and bottom to provide increased resolution in the surface and bottom boundary layers (with expansion factors of 1.14 and 1.25, respectively), fairly consistent resolution near the surface everywhere and near the bottom in water shallower than about 100 m (with layer thicknesses at the surface and bottom of about 1.0 and 0.26 m, respectively), nearly horizontal layers in the upper half of the water column, and less slope of the layers than with a regular sigma coordinate grid in the lower half of the water column (except near the bottom).
[20] The NCOM vertical grid spacing for the bottom 24 m of the water column was very similar at all 12 JRP mooring locations where the bottom depths were greater than 25 m.There were between 16 and 18 layers in the bottom 24 m, with vertical grid spacing in this region ranging from 0.2 m to 3.4 m.The distance between the bottom and the center of the deepest sigma layer varied by less than 2 cm from site to site and averaged 0.12 m. [21] Bottom stress, t0 , was calculated from bottom velocity, ũb , using the quadratic law, A logarithmic velocity profile in the bottom layer was assumed in NCOM, and, thus, C d was adjusted for slight spatial variations of the bottom layer thickness according to the equation, where k = 0.4 is von Ka ´rma ´n's constant and Dz b is the bottom layer thickness.The bottom roughness length, z 0 , was kept constant in NCOM at a value of 0.003 m.The implementation of equation ( 2) produced values of C d at the simulated JRP mooring locations from 0.011 to 0.012.
[22] NCOM was previously run for the Adriatic [Martin et al., 2006] to investigate the total circulation of the Sea with as realistic a simulation as possible.Therefore, accurate simulation of the tides was an important component of this goal as the fractional variance of measured northern Adriatic currents that could be explained by tides ranged from 10% to 64% [see Martin et al., 2006, Table 2].Tidal forcing was from tidal sea surface height and depth-averaged velocities that were prescribed at the open boundary in the northern Ionian Sea.These values were taken from the Oregon State University Mediterranean (O 1 , K 1 , M 2 , and S 2 ) and global (Q 1 , P 1 , N 2 , and K 2 ) tidal databases.Tidal potential forcing was used in the interior of the model for these eight constituents, with sensitivity studies showing that about 12% of the M 2 and 7% of the K 1 tidal elevations in the northern Adriatic could be explained by this direct astronomical forcing [Martin et al., 2006].The Adriatic NCOM version used in this paper included the total circulation as well as the tides, just as in the previous study.
[23] NCOM tidal solutions were extracted by harmonic tidal analysis (Q 1 , O 1 , P 1 , K 1 , N 2 , M 2 , S 2 , and K 2 ) using model results from 1 September 2002 to 29 April 2003.Quantitative evaluation of NCOM-simulated tidal sea surface height for the Adriatic was previously done by comparison with data from 27 International Hydrographic Organization stations [see Martin et al., 2006, Table 1].The RMS error was lower than 1.6 cm in amplitude for all constituents and was 9 and 7 degrees in phase for M 2 and K 1 , respectively.Although this previous study was done with an Adriatic version of NCOM with much less vertical resolution near the bottom than the simulation used here, the overall tidal accuracy outside the bottom zone is expected to be similar.Martin et al. [2006] also showed good qualitative agreement between the JRP mooring tidal currents and the NCOM tidal currents but no quantitative comparison was presented.

Quoddy
[24] ''Quoddy'' is a finite-element numerical model based on the 3-D nonlinear shallow water equations and using the hydrostatic, Boussinesq, and incompressible assumptions.Quoddy specifies vertical mixing using a Mellor-Yamada 2.5 turbulence closure scheme [Mellor and Yamada, 1982] with the improvements of Galperin et al. [1988].Implementation of this scheme inside Quoddy and details of the model are fully described by Lynch et al. [1996].
[25] Quoddy was implemented for simulation of the Adriatic tides using a finite-element mesh with typical node spacing of 500 m in coastal areas and 44 km in deep areas.Quoddy used sigma vertical coordinates with 21 vertical layers.Quoddy used a sinusoidal spacing of vertical levels, with highest resolution in both the surface and bottom layers and lower resolution at middepth.At the JRP locations, between 9 and 15 of these levels were located in the bottom 24 m, and vertical node spacings were between 1.0 m and 4.6 m there.Quoddy uses the quadratic law for calculation of bottom stress, but with a constant bottom drag coefficient C d = 0.003 because the bottommost node is always 1 m off the bottom everywhere.Evaluating equation ( 2) with Dz b = 2 m gives an NCOM equivalent C d at 1 m of 0.0047.The Quoddy value for C d was found through a series of numerical tuning experiments.A complete description of the Adriatic Quoddy setup is presented by Janekovic ´and Kuzmic ´ [2005].
[26] Unlike NCOM, Quoddy was not forced from tidal database values at the boundary.Instead, Quoddy was iteratively coupled to a linear, 3-D, finite-element model and inverse system, ''Truxton/Fundy'' [Lynch and Naimie, 1993], to determine sea level boundary conditions in the Strait of Otranto.The Truxton/Fundy data assimilation system used data from six coastal tide gauge stations to produce optimized boundary conditions.These boundary conditions were then used in a Quoddy run to produce residual errors at the stations, which were then in turn used in Truxton/Fundy as data to produce an update to the boundary conditions.The procedure was iterated to obtain boundary conditions specifically optimized for Quoddy.Truxton/ Fundy used l = 1 Á 10 À3 m/s as the linear frictional parameter and a constant vertical viscosity of 0.04 m 2 /s.Boundary conditions for the O 1 , P 1 , K 1 , N 2 , M 2 , S 2 , and K 2 tides were each solved for separately and then Quoddy was run separately for each constituent and once with all seven constituents together.Janekovic ´et al. [2003] conducted numerical experiments with Quoddy in the Adriatic using direct astronomical forcing.They found much less effect than Martin et al. [2006] with less than 1% contribution to M 2 elevation amplitude and 6% contribution to K 1 elevation amplitude at the northwest end of the Adriatic.Therefore, direct astronomical forcing was not used in later Quoddy runs.Janekovic ´and Kuzmic ´[2005] provides further details on the data assimilation and tide forcing procedures.
[27] Quoddy tidal solutions were extracted by harmonic tidal analysis (O 1 , P 1 , K 1 , N 2 , M 2 , S 2 , and K 2 ) using the seven-constituent model results from 1 February 1982 to 16 May 1982, with inference techniques to resolve the K 1 /P 1 and S 2 /K 2 constituent pairs.Janekovic ´and Kuzmic ´ [2005] previously validated Quoddy tidal simulations using coastal tide gauge station data, rotary current meter data at eight sites, and an ADCP record.In comparisons with data from 31 coastal tide gauge stations, the RMS error was lower than 0.8 cm in amplitude for all constituents and was 10 and 6 degrees, respectively, in phase for M 2 and K 1 [see  Janekovic ´and Kuzmic ´, 2005, Table 3].Comparisons with tidal velocity measurements also showed good agreement, especially with respect to simulation of tidal ellipse orientations.The simulated tidal current variability with depth qualitatively matched the major characteristics observed in ADCP data taken at the single location.

Theory
[28] The governing momentum equation for tidal flow in an unstratified environment over a localized area is where ũ is the horizontal velocity vector, t is time, f is the Coriolis parameter (f) times the vertical unit vector, g is gravitational acceleration, h is sea surface elevation, z is vertical height above the bottom, and A z is a coefficient of eddy viscosity.Here we have made the hydrostatic assumption, neglected horizontal advective fluxes and horizontal diffusion of momentum, and assumed that direct astronomical forcing of the tides is negligible compared to the co-oscillating tide because of the localization of the area considered.Because of the strength of the Adriatic tidal currents and their horizontal spatial scales, the horizontal advective fluxes are estimated to be <1% of @ũ @t and the horizontal diffusion of momentum fluxes are estimated at 0.1% of @ @z (A z @ũ @z ).For the northern Adriatic from late September to early May, stratification is generally weak [see Jeffries and Lee, 2007, Figure 2] and a conservative estimate of the Burger number over the entire water column for our application is 0.01.Therefore, especially in the bottom Ekman layer, stratification should play a very minor role in the momentum balance.
[29] The complication that remains is specifying the vertical structure of the eddy viscosity.The easiest choice is to assume that it is constant through a bottom boundary layer.With this choice, solutions for tidal problems with friction were found by Sverdrup [1927], Munk et al. [1970], Kundu et al. [1981], Prandle [1982], and others.If the tidal current above the boundary layer (where friction is unimportant) for a particular constituent is ũ = <{a exp(Àiwt) î + g exp(Àiwt) ĵ}, then the tidal solution in the boundary layer is ũ where a and g are complex constants, w is the tidal constituent angular frequency, the positive signs in the first exponential terms are for the case w > f (semidiurnal tides of the Adriatic), and the negative signs are for the case w < f (diurnal tides of the Adriatic).The boundary conditions that were used are a no-slip condition at the bottom and a no-stress condition at 1. Equations ( 4) and ( 5) are only valid for the Northern Hemisphere.Further sign changes are needed to transform them for the Southern Hemisphere.
[30] A best fit to these equations was sought for the ADCP data.A z was treated as an unknown and a and g were determined from tidal velocity means from ADCP depth cells more than 24 m above the bottom, excluding the depth cell nearest the surface.Equivalent Ekman depths for equations ( 4) and ( 5) differ for the clockwise (CW) and counterclockwise (CCW) rotary components of the tidal flows (see section 6), with the CW rotary depth, , always greater than both the steady flow and CCW rotary Ekman depths.Using A z = 8 Á 10 À4 m 2 /s (see section 5.1), the maximum H E for the JRP locations was 24 m and occurred at site SS2 for the CW rotary component of the K 1 tidal flows.This result, together with the observed structure of the tides in Figures 2 and 3, suggests that frictional effects do not extend significantly beyond 24 m above the bottom.The near-surface depth cell was also excluded as the tide solution was sensitive to interpolation methodology (see section 2.1).Tides for sites SS2, VR1, and VR2 were not fit to the infinite-depth equations ( 4) and ( 5) because the Ekman depth extended over the entire measured water column at these shallow sites.A z was found for each site and for the M 2 , S 2 , and K 1 tides separately by using a range of possible values from 5 Á 10 À5 m 2 /s to 2 Á 10 À3 m 2 /s in steps of 1 Á 10 À5 m 2 /s and computing a cost function as the sum of the squared errors between the data and the results of equations ( 4) -( 5) for all depths with ADCP measurements.
[31] A more complex expression for A z is to use a linear approximation and set A z = bz + k, where b is a constant coefficient with units of velocity and k is a constant kinematic viscosity.Prandle [1982] solved equation (3) applied to tides for this case.Following Prandle [1982] and forms given by Boas [1983], versions of equations ( 4) and (5) were derived for the case of linear A z and boundary conditions of no slip at the bottom and no stress at 1.They are where ker and kei are Kelvin functions of order zero.For terms in equations ( 6) -( 8) with alternative signs, the positive signs are for w < f and the negative signs are for w > f.This set of equations is only valid for the Northern Hemisphere.
[32] As was done for equations (4) and ( 5), a best fit was sought between the ADCP data and equations ( 6) -( 11).An identical procedure to the one described for those equations was used for these, with b varied from 4 Á 10 À5 m/s to 1.6 Á 10 À3 m/s in steps of 4 Á 10 À5 m/s and k varied from 1 Á 10 À5 m 2 /s to 1 Á 10 À3 m 2 /s in steps of 1 Á 10 À5 m 2 /s.b and k were optimized jointly.
[33] Equations ( 6) and ( 7), using best fit values for b and k, do not converge asymptotically to a and g for increasing z as quickly as equations ( 4) and (5) using best fit values for A z .This means that the influence of bottom friction can be responsible for slight curvatures in tidal structure far from the boundary and an infinite-depth approximation is less accurate for linear A z models than for constant A z models.Practically, the tidal ellipse parameters calculated using linear A z fits have broad curvatures above the Ekman layer such that their values slightly differ from ''1'' values throughout the entire water column at the JRP sites.We used an iterative approach to cope with this problem and to reduce the sensitivity of the b and k fits to the choices for a and g.Once the b and k best fit values were found, a and g were recalculated so that the solutions of equations ( 6) and ( 7) would pass through the means from the ADCP depth cells more than 24 m off the bottom (again excluding the depth cell nearest the surface) at the midpoint of these depth cells instead of at 1.Then, new best fit values for b and k were found using these new a and g values.For consistency of comparison, an analogous iterative procedure was also used in constant A z fits, even though iteration had less impact on these solutions.
[34] Bottom stress was calculated numerically for applications of equations ( 6) -(11) using a forward finite difference approximation for the vertical derivative of velocity evaluated at the bottom, Numerical tests showed that u 0 and v 0 as defined by equations ( 6) and ( 7) are approximately linear near the bottom for the parameter values we are using, and A z is also nearly equal to k.Therefore, the stress is nearly constant (as in a viscous sublayer), so evaluation of equation ( 12) within this zone is insensitive to the value of Dz.For the calculations of theoretical bottom stress in this paper, a Dz of 1 cm was used.
5. Comparison of Theory, Data, and Models 5.1.Fitting Bottom Ekman Layer Theory to the Data [35] Figure 5a shows the values for A z in equations ( 4) and ( 5) that gave the best least squares fit to the observations.Fits for the K 1 tides at site KB1 were discarded here and in all subsequent fits to theory because of decreasing errors toward zero friction parameter values.That is, the best fit procedures could not match the observed K 1 vertical structures at site KB1, likely because of its low signal level and signal-tonoise ratio (the mid-water-column, semi -major axis was less than 2 cm/s) and the best fits tended toward vertically uniform tides.The results did not vary much from constituent to constituent, which demonstrates the skill of equations ( 4) and ( 5) in simulating the distinctly different vertical structures of the semidiurnal and diurnal tides.The average A z values were 7.1 Á 10 À4 m 2 /s, 7.8 Á 10 À4 m 2 /s, and 8.4 Á 10 À4 m 2 /s for K 1 , M 2 , and S 2 , respectively, and 7.8 Á 10 À4 m 2 /s overall.Variation from site to site was also weak, with marginally higher values at sites SS4 -5.
[36] The thin lines in Figure 5d show the residual errors between the best fit constant A z theory and the data from the bottom 24 m of the water column.The average error was 0.31 cm/s, with higher average error (0.40 cm/s) for the stronger M 2 tides compared to the S 2 (0.25 cm/s) and K 1 (0.29 cm/s) tides.Overall, the errors peaked at 0.61 cm/s for the M 2 tides at site SS10.The K 1 and S 2 best fit solutions both had error peaks at site SS4.
[37] By factoring out b, the linear form for A z becomes b(z + k/b), where k/b is a scale height.This term cannot The square root of the mean of the sum of all the squared differences in u and v between theory and observations over the bottom 24 m and over a tidal period.The thin line is for depth-constant A z theory, and the thick line is for linear A z theory with the same symbols as in Figures 5a -5c for each.simply be neglected because ker and kei approach infinity as z approaches zero, but Prandle [1982] assumes that it is vanishingly small.This is in accord with the classic ''logarithmic layer'' approach of assigning A z = ku * z + n [Wimbush and Munk, 1970], where u * is the friction velocity and is equal to the square root of the bed shear stress divided by density and n = 1.2 Á 10 À6 m 2 /s is the molecular kinematic viscosity.Soulsby [1983] uses this form, additionally neglecting the scale height term, now n/(ku * ), by setting the condition of no slip at z = z 0 (Soulsby [1983] used z 0 = 0.0009 m) instead of z = 0. Therefore, as a first attempt to improve upon the results of optimizing equations ( 4) and ( 5), equations ( 6)-( 11) were used, k was set equal to n, and a best fit to the data was sought by optimizing b = ku * .
[38] This form for A z was expected to be a more realistic approximation than a constant A z .However, when using molecular kinematic viscosity and solving for an optimal u * on the basis of the ADCP data, this did not prove to be true.The RMS error in the bottom 24 m remaining after optimization (not shown) was lower in only 11% of the fits compared to using constant A z .Despite being worse than constant A z fits in a large majority of cases, the A z = ku * z + n fits were better in simulating the measured semi-major axis amplitudes in 74% of the cases.The reason that this result did not translate to overall better performance was that the A z = ku * z + n form better simulated semi -minor axis amplitudes in only 46% of the cases, ellipse orientation in only 6% of the cases, and phase in only 11% of the cases.
[39] Because of this result, the fits were redone allowing both b and k to vary as specified in section 4. For the S 2 fit at site SS4, k was varied to higher values (2 Á 10 À3 m 2 /s) in order to find its best fit value close to the maximum k considered for other sites.Figures 5b and 5c show the resulting optimum values for these parameters.Values for K 1 at sites VR4 -5, for S 2 at site SS6, and for all constituents at site CP2 are not shown because the fit selected the lowest value of b (4 Á 10 À5 m/s) and thus a true minimum was not found.Clearly for these cases, the best fit was approaching the constant A z case, which is mathematically equivalent to using equations ( 6) -( 11) with b = 0.For other cases, the optimum b occupied a fairly small range of values independent of coefficient and site.Values for K 1 and S 2 at site SS4, for K 1 at site SS9, and for M 2 at site SS10 were somewhat larger, with values more than two times the overall average.Average b values were 3.2 Á 10 À4 m/s, 2.7 Á 10 À4 m/s, and 2.7 Á 10 À4 m/s for K 1 , M 2 , and S 2 respectively and 2.9 Á 10 À4 m/s overall.
[40] Optimum k values (Figure 5b) were more than 100 to more than 800 times greater than molecular kinematic viscosity.Optimum k also occupied a fairly small range of values independent of coefficient and site, with a weak trend toward lower values from southwest to northeast along the SS line.Average k values were 4.6 Á 10 À4 m 2 /s, 4.7 Á 10 À4 m 2 /s, and 5.7 Á 10 À4 m 2 /s for K 1 , M 2 , and S 2 respectively and 5.0 Á 10 À4 m 2 /s overall.
[41] The thick lines in Figure 5d show the residual errors between the best fit linear A z theory and data.We have included the errors for fits with the minimum b for comparison to the constant A z residual errors.Linear A z theory has lower errors than constant A z theory in 91% of the cases.The average error decreased by 30% to 0.22 cm/s.Average errors for particular constituents were 0.22 cm/s for K 1 , 0.26 cm/s for M 2 , and 0.19 cm/s for S 2 .Comparing the error for each tidal current parameter, linear A z better simulated the semi -major axis in 91% of the cases, the semi-minor axis in 74% of the cases, the ellipse orientation in 51% of the cases, and the phase in 83% of the cases.Considering only the ellipse orientation for K 1 and only the phase for M 2 and S 2 (i.e., the respective dominant angle changes), linear A z performed better in 83% of the cases.
[42] Figures 6 and 7are representative of the fits.For both diurnal and semidiurnal tides, constant A z theory (black curves) produces a distinct bulge in the semi-major axes and semi -minor axes at the top of the bottom Ekman layer analogous to the increase in current speed produced for a steady, constant A z , bottom Ekman layer.However, such a bulge is generally not seen in the ADCP data and a linear A z (red curves) is better able to match the currents in this area because it also lacks a distinct bulge and instead has a broad maximum.Both constant A z and linear A z theory match well the sharp changes in diurnal ellipse orientation and semidiurnal phase observed in the bottom ADCP depth cells, but constant A z theory generally overpredicts the change around the depth level of the bottom depth cell.Note that it is the constant component of linear A z theory that is primarily responsible for reproducing these angular changes.Although linear A z theory matches the semi -major axis, semi-minor axis, diurnal ellipse orientation, and semidiurnal phase changes more accurately than constant A z theory in most cases, weak changes in the semidiurnal ellipse orientation and in diurnal phase match constant A z theory better in half the cases.
[43] Sensitivity studies were performed to investigate the impact on the best fit values from using different cost function forms.Specifically, 20 m and, at some sites, 30 m were tried as selection heights for the ''Ekman layer'' and calculation of a and g.Also, for these fits (and an additional fit with a 24-m Ekman layer) the cost functions were evaluated only over the bottom layer and not over the entire water column, and no iteration of a and g was performed.
The results were relatively insensitive to these variations in methods and cost function definitions.Average b and k best fit values (not including three outliers) for the 20-m cost function were 3.3 Á 10 À4 m/s and 4.2 Á 10 À4 m 2 /s, respectively.The average constant A z value for this cost function was 8.0 Á 10 À4 m 2 /s.

Tidal Current Vertical Structure Comparison
[44] The degree of mismatch between the observations and model simulations in the bottom Ekman layer is sensitive both to how the model represents the bottom layer and to the model solution errors above the Ekman layer.Separation of these two effects is desirable but cannot be simply achieved.The approach we have taken is to estimate biases in the upper water column and remove these from the entire water column before calculating error statistics in the Ekman layer.This is an imperfect solution since, using linear A z theory for an example, it removes bias from the a and g terms in equations ( 6) and (7) but does not remove the effect of model solution error in B 0 and D 0 caused by incorrect a and g values in equations ( 8) and ( 9).Thus the depth structure of the tides depends on a and g in a complex way.Sometimes removal of the upper water column bias will cause error in the Ekman layer to increase.However, all errors were calculated with and without bias removal and, overall, the errors are lower when the bias is removed.Bias removal increases the NCOM average (over sites) RMS error in the Ekman layer for the K 1 tide, but decreases this RMS error for all the Quoddy tides and for the M 2 and S 2 NCOM tides.Therefore, to better examine how the models represent the Ekman layers, errors with upper water column bias removed are the ones that are used.
[45] Figures 8 -10 show the calculated upper water column biases from NCOM and Quoddy.These were calculated as the difference in ellipse parameters using average a and g values calculated from the models and data over the same depth ranges determined from ADCP depth cells higher than 24 m off the bottom.Quoddy tends to have semi -major axes that are biased too large and removal of these biases greatly helps to lower RMS errors for the Ekman layer.NCOM M 2 and S 2 tidal currents all lead the observed tidal currents except at site KB1.Removal of this phase bias greatly lowers the NCOM Ekman layer errors.The K 1 biases for the semi -major axes and semi -minor axes are especially interesting as both NCOM and Quoddy have, uniquely for these cases, very similar patterns of bias.This suggests the presence of a dynamics or parameter error that is common to both models and affects the K 1 tidal currents without causing large effects in the M 2 or S 2 tidal currents.
[46] Figure 11 shows the total RMS error for each tidal constituent evaluated for the layer within 24 m of the bottom.The bias correction described in the previous paragraphs has been applied.The error from the best fit for linear A z from Figure 5d is also shown for comparison.The average errors for NCOM are 0.37 cm/s for K 1 , 0.36 cm/s for M 2 , and 0.33 cm/s for S 2 (0.35 cm/s overall) and for Quoddy are 0.33 cm/s for K 1 , 0.44 cm/s for M 2 , and 0.30 cm/s for S 2 (0.36 cm/s overall).Results for linear A z theory averaged 38% lower than these values.NCOM and Quoddy have similar error levels, but, near Istria, where M 2 tidal currents are especially strong and the observed tidal structure changes strongly with depth, NCOM matches the currents exceptionally well while Quoddy has only weak tidal current depth changes.Although linear A z theory, NCOM, and Quoddy all match the observed tidal structure in the bottom Ekman layer relatively well, better agreement could be achieved before reaching the error level of the measurements, as directly illustrated in Figures 6 and 7 for site SS6.
[47] Figures 12 -14 show the breakdown of the Ekman layer RMS error according to tidal ellipse parameter.NCOM averaged the lowest errors for M 2 semi-minor ellipse axes, Figure 6.Observations, theoretical fits, and model simulations of K 1 tidal currents at SS6.The four panels show tide ellipse parameters versus depth as individually labeled.The blue curves are measured values with dashed lines indicating 95% confidence limits, the black curves are from the constant A z best fit, the red curves are from the linear A z best fit, the green curves are from NCOM, and the magenta curves are from Quoddy.The model simulations have been interpolated to ADCP-observed depth levels and the upper water column biases have been removed.Note that the x axis scale for ellipse orientation has a larger range than the x axis scale for phase because the dominant diurnal angular response is in ellipse orientation.but linear A z theory averaged the lowest errors for all other constituents and ellipse parameters.Quoddy simulated the M 2 and S 2 semi-major axis changes relatively well, with only 5% higher average errors than linear A z theory but this result strongly depends on the removal of the Quoddy upper water column bias (e.g., see Figure 10).In contrast, the removal of the upper water bias greatly increases the Ekman layer RMS errors for the K 1 semi-major axis at sites SS4 -5 for NCOM and site SS4 for Quoddy (i.e., the three highest errors for this parameter in Figure 12), but average (over stations) K 1 semi -major axis Ekman layer RMS errors are reduced by bias removal.There is a distinct pattern in the degree that linear A z theory averaged lower errors for the different ellipse parameters, with overall 22% and 17% lower errors for semi -major ellipse axes and semi -minor ellipse axes, but 47% and 53% lower errors for ellipse orientation and phase.In particular for the semidiurnal tides, linear A z theory averaged 62% lower errors in phase.

Time and Vertical Structure of Eddy Viscosity
[48] The analytical solutions of section 4 assume very simple vertical structures for A z and also assume that A z is constant in time.However, both NCOM and Quoddy fully simulate A z with complex vertical and temporal structures using Mellor-Yamada 2.0 and 2.5 turbulence closure schemes, respectively.These temporal and depth structures were extracted from the models for analysis and comparison.A z was saved at hourly intervals along with other variables for the main run of NCOM.A z was derived using other saved turbulence parameters for Quoddy from a dedicated 172-day run made without nodal modulation.
[49] Figure 15 shows the site-averaged power spectral density (psd) of A z for the bottom depth levels of NCOM and Quoddy.They were calculated from the A z time series with their means removed, using Welch's averaged periodogram method over block lengths of 1024 h ($43 days), with 50% overlapping Hanning windows and no detrending.The frequency structures of the NCOM and Quoddy results are similar, with psd peaks clustered around particular frequency bands near 0, 1, 2, 3, and 4 cycles per day (cpd).Quoddy often has higher psd in these bands and has further psd peaks at higher frequencies that are not excited in NCOM.In contrast, NCOM has higher psd between bands as would be expected from a model simulating the total circulation rather than only tides.Both in NCOM and Quoddy, strong peaks occur at particular frequencies (marked with dotted lines in Figure 15) within each cluster, and each of these peaks can be related to specific frequencies of tidal constituents or tidal interactions as annotated.
[50] Figure 16 shows the time means of A z in the Ekman bottom layer for NCOM (red) and Quoddy (blue) at select sites.NCOM has higher gradients and higher values of time mean A z than Quoddy throughout this layer.Both NCOM and Quoddy have considerable curvature in their time mean A z depth structures, but the curvatures in the NCOM structures tend to occur further off the bottom or at higher A z values and are therefore not as evident in the SS8, CP3, and  6 but for M 2 tidal currents at SS6.Note that the x axis scale for phase has a larger range than the x axis scale for ellipse orientation because the dominant semidiurnal angular response is in phase.
VR6 panels.The theoretical fits (orange, black, and magenta lines) show considerable spread at these scales from each other and from site to site, but generally have smaller gradients close to the bottom than those of the Quoddy or NCOM A z time means.
[51] Because of the time dependence of A z , there are frictional impacts at all depths not described by these time means.Because the depth structure of A z and the depth structure of @ũ @z both change with time and with respect to each other, their interaction in equation (3) will vary in time and we cannot separate their frictional impacts from each other.One way the effect of A z can be roughly estimated is to calculate a gain in the semi -major ellipse axis at tidal frequencies due to the action of A z .This was done by dividing the semi -major axis values of harmonically analyzed kinematic shear stress, (A z @ũ @z ), by the semi -major axis values of harmonically analyzed velocity gradient, @ũ @z .The results for K 1 (green) and M 2 (cyan), and for NCOM (solid) and Quoddy (dashed) are shown in Figure 16.For Quoddy, these gains are always higher than the Quoddy time mean A z values and often approach the NCOM time mean values and structure at depths close to the bottom.In contrast, the NCOM gains are only slightly higher than the NCOM time mean values for most sites near the bottom, and they abruptly shift farther up the water column to lower and nearly depthconstant values.Above the shift, the NCOM gains are noisy because of less well determined semi -major axis values of @ũ @z .For graphical clarity, gains where the estimated error were !50% of the estimated value for either axis are not plotted in Figure 16.
[52] Although the ''impact'' of A z alone can only be approximated, the combined action of A z and @ũ @z , i.e. the frictional kinematic shear stress, for tidal momentum at specific frequencies can be exactly calculated and analyzed.Here we exploit the fact that the Fourier Transform of the multiplication of two time-dependent signals is the convolution integral of the Fourier Transforms of the individual signals.Therefore, for a particular tidal frequency of interest, the Fourier coefficient of the frictional kinematic shear stress, tk , will be where F represents Fourier transformation.
[53] By this expression, if either A z or @ũ @z have dominant spectral peaks, then only a few Fourier coefficients would determine the frictional kinematic shear stress that impacts a particular tidal constituent's momentum.The convolution integral shifting property causes the time mean of A z to interact with the ±w peaks of @ũ @z , and unless @ũ @z has high psd at other frequencies, which when shifted correspond to other peaks in A z , then the frictional kinematic shear stress for that tidal constituent will be mainly determined by the   time mean of A z .This is the case for near bottom frictional kinematic shear stresses in NCOM and Quoddy and is why the gains track the time means in Figure 16.
[54] Figure 17 shows this in a different form.Here, the cumulative sums of the discrete version of equation ( 13) for the direction of strongest tk are plotted for NCOM and Quoddy at M 2 and K 1 frequencies and at 1.5 and 15 m height above the bottom.Jumps in the cumulative sums at particular frequencies indicate their relative importance in determining frictional kinematic shear stress.The largest jumps are all at zero frequency.Because NCOM simulates the strong Western Adriatic Current flow at SS4 and the time mean of @ũ @z interacts in the convolution integral with the tidal peaks of A z , there are also large jumps at 2 cpd for M 2 and 1 cpd for K 1 .Quoddy completely lacks these jumps because it does not simulate the mean currents.In the upper Ekman layer, the Quoddy kinematic shear stresses are weaker but the frequency contribution structure is similar to that near the bottom.This is not true for NCOM, where the higher A z psd between bands at lower frequencies often interacts with various shears to interfere destructively and thereby reduce the frictional kinematic shear stress that would have otherwise been established by a time constant A z .This likely explains the abrupt shift and noise noted earlier in the NCOM gains (Figure 16).

Bottom Stress Comparison
[55] Bottom stress at the mooring sites can be calculated from each of the models and from the best fit linear A z theory results.Figure 18 shows examples of this calculation for the M 2 and K 1 tides at moorings SS8 and VR6.The theoretical kinematic bottom shear stresses (solid black curves) were calculated from equation (12) using the best fit b and k values, and thus the time evolution of the shear stresses takes the form of an ellipse that is shaped, rotated, and phased in accord with the near-bottom tidal currents.The black dashed lines show bottom shear stresses that would be calculated using theoretical values at 1-m height in equation ( 1) with C d optimized to produce the same maximum shear stresses as the solid black ellipses.Use of a quadratic drag law on currents tracing an ellipse in time leads to a double lobed bottom shear stress structure; that is, the shear stresses depart from an elliptical shape because of nonlinearity.Also, the orientation and phase of the quadratic drag law stresses will match the currents at the level where they are evaluated and miss any changes taking place deeper in the water column.The values for C d needed to match the maximum shear stress magnitudes of equation ( 12), averaged over all the theoretically fitted sites, were the seemingly high values of 0.022 and 0.051 for M 2 and K 1 , respectively.
[56] Also shown in Figure 18 are the bottom shear stresses from three models.Generally, NCOM and Quoddy produce stresses that are similar to each other, but smaller than the theoretical result.Both have double lobed shapes from the use of the quadratic drag law.Except for the K 1 shear stresses at site SS4, NCOM has larger maximum M 2 and K 1 bottom shear stresses than Quoddy at all sites.However, in 19 out of these 22 cases, the percentage increase is less than what would be expected simply because of the effectively 58% larger (1-m) bottom drag coefficient in NCOM.Thus this suggests that NCOM tidal currents are somewhat weaker at that level than Quoddy tidal currents.The percentages by which NCOM bottom shear stresses exceed those of Quoddy generally increase toward the northwest end of the basin, with averages for the various mooring lines of 22% for SS, 39% for CP, and 54% for VR.Note that for both NCOM and Quoddy, the actual bottom shear stresses associated with the tides will be higher than the values in Figure 18 since the double lobed ellipses were calculated using single tidal constituent velocities only and the nonlinearity of equation (1) will produce cross terms with other tides for Quoddy and with other tides and other currents for NCOM.
[57] The blue ellipses in Figure 18 are from a depthaveraged, linear, shallow water equation tide model using strong-constraint variational data assimilation of the tidal observations in this paper.The model is described by Griffin and Thompson [1996] and the setup for the northern Adriatic is described in Book [2007].The linear friction parameter was varied to improve agreement with the data and 5 Á 10 À4 m/s was determined as the optimal value for this model.In the linear, depth-averaged model equations, bottom stress is simply this friction parameter multiplied by the depth-averaged tidal currents.Because of this, the bottom stresses for a particular tidal constituent trace out a true ellipse.Also, since the bottom stresses are based on depth-averaged currents, they account for very little of the broadening, rotating, and phase shifting that takes place in near-bottom tidal currents.Therefore, they are much too weak in the cross-axis direction of the Adriatic.Contrastingly, the optimized stresses for this model seem much too large in the along-axis direction compared to other bottom shear stress estimates.

Discussion
[58] Our observations of near-bottom tides show common general characteristics relative to mid-water-column tidal currents: (1) tidal ellipses shorten and broaden, (2) diurnal tidal current ellipse orientations rotate strongly CCW, and (3) semidiurnal current phases advance strongly.All these can be qualitatively explained by considering the rotary components of the tidal flow separately and examining each of them in a frame of reference that is rotating such that the flow becomes stationary.In these reference frames the steady current Ekman solution above a rigid surface [Kundu, 1990] applies, and the speed and angular changes of steady currents with depth are equivalent to the amplitude and phase changes of the respective rotary components.The steady current Ekman solutions are in Earth coordinates and are already in a rotating coordinate system.Therefore, additional rotation at ±w to match the tidal rotary component rotation only acts to increase or decrease the equivalent ''Coriolis parameter.'' [59] Thus, for the CCW rotating Northern Hemisphere, CCW rotary components have an Ekman depth of , CW rotary components have an Ekman depth of p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , and both rotary components phase angles rotate CCW with depth.However if w > f (semidiurnal tides of the Adriatic), the reference frame itself is rotating CW for CW rotary components and ''Southern Hemisphere'' steady current Ekman solutions must be used, causing their phase angles to rotate CW instead of CCW.
[60] For near-reversing tides in the Adriatic, mid-watercolumn tides have rotary components that are nearly equal in amplitude (required to produce reversing tides) and have distinct times and orientations when the rotary vectors align in direction producing maximum currents and setting the ellipse orientations and phases.The characteristics (1 -3) outlined at the beginning of this section can be explained as follows.Because of their larger Ekman depths, the CW rotary component amplitudes have decay scales which extend further up in the water column than the CCW amplitudes.Consequently (1) the tidal ellipses broaden as they shorten, and tidal current vectors rotate CCW with time near the bottom.For diurnal tides, both CW and CCW rotary component have phase angles that rotate CCW with depth (albeit at slightly different rates) and, because they rotate together, (2) changes in timing of the maximum current (phase) are small but the ellipse orientations rotate CCW along with the rotary component phases.For semidiurnal tides, CCW rotary component phase angles rotate CCW with depth but CW rotary component phase angles rotate CW with depth and therefore (3) the rotating vectors will align at nearly the same orientation but the timing of the alignment (phase) will advance.In short, differing decay scales for rotary component amplitudes cause effect 1, and rotation of rotary component phases either together or oppositely cause effect 2 or 3.
[61] Both constant A z theory and linear A z theory reproduce such depth changes.However, in general, constant A z theory does not match well the observed curvature of the semi-major axis and semi -minor axis profiles, and linear A z theory (with a constant k much larger than molecular viscosity) matches better the observed structure of all the most important tidal parameter changes in the Ekman layer.The average scale height (where bz = k) for optimized b and k values found in this study was 2.9 m.This implies an extended region away from the bottom where eddy viscosity is relatively high and approximately constant.Of course, since the ADCP measurements do not extend down into this layer, the actual structure near the sea bed could be different and A z could have a different form than the one that fits tidal velocities 2.5 m and higher off the bottom.The VR4 mooring was configured with a higher-frequency ADCP and it measured velocities as close as 1 m off the bottom as part of its bottom depth cell, yet error levels were not anomalous (Figure 5d) compared to other sites with measurements farther off the bottom.[62] Bottom roughness is one possible explanation for high eddy viscosity near the bottom.If the roughness elements extend beyond a hypothetical viscous sublayer, molecular viscosity is unimportant even near the bed [Soulsby, 1983].Of course, the moorings themselves could act as roughness elements.As described by Perkins et al. [2000], BARNY mounts are shaped like large barnacles with a 2 m diameter circular footprint and 0.5 m maximum height.The effect of a BARNY mount on the tidal currents is unknown, but some perturbation of the flow field near the sea bed should be expected.
[63] Werner et al. [2003b] used Benthic Acoustic Stress Sensor tripod deployments at 76 m depth on the southern flank of Georges Bank to determine M 2 tidal friction shear velocities (u * ) through logarithmic fits of observed nearbottom velocities.They used standard logarithmic layer theory to derive bottom stress, a bottom drag coefficient, and bottom roughness.M 2 bottom shear stresses estimated from their experiment averaged 1.6 Á 10 À4 m 2 /s 2 [Werner et al., 2003a], 11 times larger than the average from the northern Adriatic best fit linear A z theory results.This factor is similar to the average squared ratio, 15, of M 2 current semi -major axis amplitudes in their study and this one, so, as expected, u * scales roughly as the current speed.Inserting our time mean u * values in equation (C2) from Werner et al. [2003a], gives a logarithmic layer thickness of $2.0 m, which is below the depth range observed by most ADCPs used in our study.For unstratified conditions, Werner et al. [2003a] found good agreement between measured tides above the logarithmic layer and a constant A z model (linear A z below the log layer).But their model did not produce a distinct bulge in semi -major axes and semi -minor axes as the Adriatic fits did, probably because they used A z values ($0.03 m 2 /s) approximately 40 times the northern Adriatic best fit values, with an associated Ekman depth >100 m spreading the bulge over a much larger depth range.However, the use of such large constant A z values produces poor fits to the observed Adriatic tides.
[64] When upper water column biases are removed from NCOM and Quoddy, (particularly the amplitude bias in Quoddy and the phase bias in NCOM) they too match the observed currents in the bottom Ekman layer well.The Mellor-Yamada level 2.0 and 2.5 turbulence closure schemes implemented in these models seem to reproduce vertical eddy coefficients that are well suited for mimicking realistic tidal changes near the bottom.Although the errors for use of linear A z theory are smaller, it is an optimized best fit from site to site and constituent to constituent which is reasonable because of the potential for spatial changes in bottom roughness, but is a luxury not enjoyed by NCOM or Quoddy.However the fact that optimized linear A z theory consistently gives much better results, specifically for semidiurnal phase and diurnal ellipse orientation, suggests a dynamical difference.Linear A z theory, with use of a relatively high value for k, seems to be particularly adept at matching such angular changes.The sizes of the theoretical bottom shear stress ellipses in Figure 18 are primarily dictated by the values of k and if these extrapolated results are valid, they suggest that the bottom shear stresses of the two models are too low.
[65] One possible explanation is that use of equations ( 1) and (2) assumes a logarithmic velocity profile.Clearly equations ( 6) -( 11) depart from this in agreement with the general finding of Soulsby and Dyer [1981], who suggest that in an accelerating (but nonrotating) tidal flow the nearbed velocity profile departs from the usual logarithmic form.In contrast, Lueck and Lu [1997] in their tidal channel measurements found that the velocity profile above 3 m off the bottom in the along channel direction matched a logarithmic form well and that departures from their fits were not consistent in this depth zone with Soulsby and Dyer [1981] acceleration corrections.But, they also found that the velocity profiles in the across channel direction matched a linear profile instead of a logarithmic one.Figure 18 from our results shows that bottom stress estimates from the quadratic law are particularly low relative to estimates from equation ( 12) in the semi -minor ellipse direction.Another effect to consider is form drag [Chriss and Caldwell, 1982] which Lueck and Lu [1997], Ullman and Wilson [1998], and Werner et al. [2003b] all used to explain their findings of larger roughness lengths or drag coefficients than typical values.However, it is less clear how this could explain our general finding of implied higher bottom stress, as the sites near Italy had very muddy bottoms and typical bed forms responsible for form drag are not expected to be present.
[66] It is interesting to note that despite a wide spread in A z gradients between the various model and theoretical solutions shown in Figure 16, all the solutions seem to converge to close to the same A z values in the bottom 1 m.If the vertical derivative is carried through the friction term in equation (3), then, for a given depth level, it is the numerical value of A z that multiplies @ 2 ũ @z 2 and the slope of A z that multiplies @ũ @z , which together with the convergence of solutions to a particular A z value, suggests the relative importance of @ 2 ũ @z 2 for determining the tidal structure in this region.If the velocity structure above this level is roughly Figure 17.Cumulative sum versus frequency for the convolution integral between the Discrete Fourier Transform of A z and the Discrete Fourier Transform of @ũ @z at sites SS4 and SS9.The colors indicate NCOM M 2 (red), NCOM K 1 (cyan), Quoddy M 2 (magenta), and Quoddy K 1 (blue).The solid lines are from the level closest to 1.5 m off the bottom (circles mark the zero frequency values), and the dashed lines are from the level closest to 15 m off the bottom (crosses mark the zero frequency values).In each case, the component of tk with the strongest stress is the one shown.logarithmic, then @ 2 ũ @z 2 and @ũ @z should have opposite signs, and therefore higher numerical values of A z in Quoddy and NCOM could offset the effects of larger A z slopes and thus explain how the diverging A z depth profiles result in similar tidal solutions.
[67] The analysis of the time structure of A z from Quoddy and NCOM and the fact that theoretical fits with timeconstant A z match the observations best show that time variation of A z is not very important for Adriatic tidal momentum balances.As the CW rotary component of the Adriatic tidal flows decays slower toward the bottom than the mean flow, tidal components become more dominant in the spectra of A z and @ũ @z , and the interaction of the time mean of A z and the tidal components of @ũ @z predominately determine the frictional tidal shear stresses in this region where the stresses tend to have maximums.For models without mean or low-frequency flows, such as this application of Quoddy, the role of the time mean of A z is even greater in determining the frictional shear stresses for tides.
[68] Shear stress in the cross-axis direction of the Adriatic is particularly important for diurnal tides.Diurnal tidal currents near the bottom rotate with depth more toward this direction than semidiurnal tidal currents, and the diurnal tidal waves propagate in this direction from the northeast to the southwest coasts unlike the semidiurnal tide waves, which propagate along the axis of the sea [Malac ˇic ˇet al., 2000].These facts may account for the matching bias structure for K 1 amplitudes in NCOM and Quoddy.Both models use spatially uniform bottom roughness and therefore may miss the potentially significant effect on K 1 of a northeast to southwest bottom drag difference caused by varying sediment types (sand to mud).The trend for the SS line in Figure 8 (top left) is in the right sense for this where a moderate and uniform model C d would not damp enough energy near the northeast coast but damp too much near the southwest coast.Lower bottom shear stress in the semiminor axis direction may also have an effect.

Conclusions
[69] A large observational database unmatched in previous studies was compiled and used with a suite of mathematical and numerical models in an effort to explore Adriatic tidal dynamics.Fifteen bottom-mounted ADCPs deployed for more than six months in the northern Adriatic during a time period with generally little stratification were able to resolve strong tidal current structure changes at 1-m intervals from Figure 18.Kinematic bottom shear stress from theory and models at sites SS8 and VR6 for (left) M 2 tides and (right) K 1 tides.Red is from Quoddy, green is from NCOM, blue is from a model with tuned linearized friction, and black is from the best fit linear A z theory.The black solid lines are from equation (12), while the black dashed lines use quadratic drag law from the 1-m height theoretical results with C d chosen to yield the same maximum t 0 /r magnitude as the solid ellipses.Dots indicate relative phasing.The direction of rotation of t 0 /r with time is from the dots to the gaps.heights of 3 m off the bottom upward to heights where the tidal current structure was nearly vertically uniform.Tidal currents in much of the water column were near-reversing tides, with M 2 currents having the largest amplitudes and S 2 and K 1 currents both having significant energy.Near the bottom, tidal current ellipses were all shortened and broadened, semidiurnal currents led upper water column currents, and diurnal tidal current ellipse orientations rotated CCW descending through the boundary layer.
[70] Such changes are in accord with tidal bottom Ekman layer theory and match solutions using a linear form for A z that are obtained here for conditions of no slip at the bottom and no stress at 1. Linear A z solutions have 30% smaller errors than constant A z solutions, but only if the constant term in the linear form is far above molecular values.NCOM and Quoddy simulations were also compared to observations.Accounting for upper water column biases, the model simulations in the bottom Ekman layer accurately reproduce the major tidal structure changes with depth, although not as well as optimized theory.
[71] Analysis of the time structure of A z solved by using Mellor-Yamada turbulence closure in NCOM and Quoddy show that the time mean of A z predominately controls the frictional tidal shear stresses.Comparison of the depth structure of A z between the model solutions and the theoretical fits shows a large spread of values in the bottom 24-m layer, but also shows convergence toward similar values at depths less than 1 m off the bottom.The optimized theory fits show that the constant part of A z was most important in matching the tidal orientation and phase changes and a nonzero A z slope was important in matching the tidal amplitude changes.
[72] A comparison of kinematic bottom shear stresses from the optimized theory, from Quoddy, from NCOM, and from a data assimilation model using an optimized linear friction parameter shows a wide range of values produced by these different models in simulating the northern Adriatic tides.Quoddy tends to have the smallest stresses, but NCOM stresses are similar.Theoretical stresses are larger, especially in the cross-axis direction.The model using linearized friction has the largest stresses in the along-axis direction but weakest in the cross-axis direction.For near-reversing tidal currents, using a quadratic drag law or a linear friction parameter can produce low bottom shear stresses in the semi -minor axis direction.
[73] Further work is needed to explore the impact of potentially underestimating this directional component of stress in numerical models of the Adriatic and elsewhere.Also, the northern Adriatic with its general lack of stratification during half of the year and significant semidiurnal and diurnal tides would be a good location to test the validity of equations ( 6) -( 11) near the bed from examination of appropriate current measurements from long-term, bottom tripod deployments.

Figure 2 .
Figure2.Vertical structure of M 2 tidal currents for all sites except KB1 and VR1-2.Conventions for tidal current parameters are those ofForeman [1978].

Figure 3 .
Figure 3.As in Figure 2 but for K 1 tidal currents.

Figure 4 .
Figure 4. Measured tidal current ellipses (cm/s) at station SS6 for four different depths.M 2 ellipses are drawn in magenta, K 1 ellipses are drawn in red, S 2 ellipses are drawn in yellow, O 1 ellipses are drawn in blue, P 1 ellipses are drawn in green, N 2 ellipses are drawn in cyan, and K 2 ellipses are drawn in black.Dots indicate relative phasing.The directions of rotation of the ellipses with time are from the dots around the ellipses to the gaps.Distance from the bottom to the midpoint of the ADCP depth cells are 60, 41, 22, and 3 m as labeled.

Figure 5 .
Figure 5. Best fit values and RMS errors for theoretical solutions.(a) Best fit, depth-constant A z values with circles for K 1 , dots for M 2 , and asterisks for S 2 .Best fit (b) k and (c) b values for linear A z theory with triangles for K 1 , diamonds for M 2 , and stars for S 2 .(d)The square root of the mean of the sum of all the squared differences in u and v between theory and observations over the bottom 24 m and over a tidal period.The thin line is for depth-constant A z theory, and the thick line is for linear A z theory with the same symbols as in Figures5a -5cfor each.

Figure 7 .
Figure7.As in Figure6but for M 2 tidal currents at SS6.Note that the x axis scale for phase has a larger range than the x axis scale for ellipse orientation because the dominant semidiurnal angular response is in phase.

Figure 8 .
Figure8.Estimated model K 1 bias for currents from heights more than 24 m off the bottom.Squares are for NCOM, and triangles are for Quoddy.Positive ellipse orientation bias indicates that the model tidal ellipse is rotated counterclockwise with respect to the observed tidal ellipse.Positive phase bias indicates that the model tidal currents lag the observed tidal currents.Mooring names are shortened for ease of display.

Figure 9 .
Figure 9.As in Figure 8 but for M 2 tidal currents.

Figure 11 .
Figure 11.Total RMS error for all tidal currents from heights 24 m or closer to the bottom.Currents were corrected for upper water column bias (Figures 8 -10) before computation of RMS.RMS errors for (top) K 1 , (middle) M 2 , and (bottom) S 2 are shown.Squares are from NCOM, triangles are from Quoddy, and asterisks are from best fit linear A z theory.

Figure 12 .
Figure12.K 1 tidal ellipse parameter RMS errors for all currents from heights 24 m or closer to the bottom.Currents were corrected for upper water column bias (Figure8) before computation of RMS.Squares are from NCOM, triangles are from Quoddy, and asterisks are from best fit linear A z theory.

Figure 14 .
Figure 14.As in Figure 12 but for S 2 tidal currents.The off-scale RMS errors for KB1 ellipse orientation and phase are 16.5°and 15.4°, respectively.

Figure 15 .
Figure15.Average power spectral density of A z for all 12 JRP mooring locations where the bottom depths were greater than 25 m.The thick lines are calculated from the Quoddy model (at 1 m off the bottom), and the thin lines are calculated from the NCOM model (at $0.26 m off the bottom).The dotted lines mark specific frequencies of tidal importance as annotated.

Figure 16 .
Figure16.Vertical structure of various estimates of A z at sites SS5, SS8, CP3, and VR6.The theoretical linear fits to observations are shown for the K 1 (orange), M 2 (black), and S 2 (magenta) constituents.The time mean values for NCOM (red) and Quoddy (blue) are drawn with dots to indicate the depth levels where A z is defined.The estimated gains at the K 1 (green) and M 2 (cyan) tidal frequencies are drawn as solid lines for NCOM and as dashed lines for Quoddy.Gains with estimated errors that are greater or equal to 50% of the gain value are not drawn.

Table 1 .
Mooring Positions and Depths