Nonlinear Processes in Geophysics Lagrangian structure of flows in the Chesapeake Bay: challenges and perspectives on the analysis of estuarine flows

In this work we discuss applications of Lagrangian techniques to study transport properties of flows generated by shallow water models of estuarine flows. We focus on the flow in the Chesapeake Bay generated by Quoddy (seeLynch and Werner , 1991), a finite-element (shallow water) model adopted to the bay by Gross et al. (2001). The main goal of this analysis is to outline the potential benefits of using Lagrangian tools for both understanding transport properties of such flows, and for validating the model output and identifying model deficiencies. We argue that the currently available 2-D Lagrangian tools, including the stable and unstable manifolds of hyperbolic trajectories and techniques exploiting 2-D finite-time Lyapunov exponent fields, are of limited use in the case of partially mixed estuarine flows. A further development and efficient implementation of three-dimensional Lagrangian techniques, as well as improvements in the shallow-water modelling of 3-D velocity fields, are required for reliable transport analysis in such flows. Some aspects of the 3-D trajectory structure in the Chesapeake Bay, based on the Quoddy output, are also discussed.


Introduction
The Chesapeake Bay is the largest estuary in the United States, stretching for about 320 km from the mouth of the Susquehanna River (Havre de Grace, Maryland) to its entrance at the Atlantic Ocean (Cape Charles, Virginia). About 15.1 million people live in its watershed and the bay provides habitat for more than 3600 species of plants and animals. Despite its strategic location and commercial importance, relatively little is known about the flow dynamics in the bay and Correspondence to: M. Branicki (m.branicki@bristol.ac.uk) its coupling to external forcing. The Chesapeake Bay poses a challenging place to develop a predictive model largely due to its elongated and and narrow shape, its complex coastline and bathymetry, and a large network of tributaries (see Fig. 1). Tidal forcing in the bay is relatively weak, compared to other large estuaries, with tidal range rarely exceeding 1 m (Browne and Fisher, 1988). However, strong storm winds in the autumn can occasionally significantly influence the dynamics in the bay and destratify the entire water column (Goodrich et al., 1987;. Problems with a reliable prediction of storm surge in the bay are also well documented (e.g. aftermath of hurricane Isabel 1 in 2003; see Shen and Wang, 2005;Shen et al., 2006). Moreover, the salinity and temperature distributions, which represent two important variables regulating the bay's biological activity, are notoriously difficult to model and are sensitive to the turbulence closure schemes used (see Xu et al., 2002;Guo and Valle-Levinson, 2007;Li et al., 2005;Gross et al., 2001). Marine biologists, in particular, need the estimated salinity distribution to map the habitat of species and the location of the oxygen-depleted, "dead zones" and to design monitoring protocols more efficiently.
A number of models have been developed, or configured, over the past 20 years in order to simulate the flow in the Chesapeake Bay (e.g. CH3D, Johnson et al., 1993;Johnson, 2000, Princeton Ocean Model, POM, Guo andValle-Levinson, 2007;Regional Ocean Modelling System, ROMS, Li et al., 2005;Quoddy Gross et al., 2001). Validation and skill assessment of these models are generally carried out via comparisons of hindcast simulations with time series measurements of high resolution hydrographic data, shipboard measurements and Doppler radar imagery. While such data provide very important reference when analyzing the reliability of model predictions of the (Eulerian) velocity fields and temperature and salinity distributions, they do M. Branicki and R. Malek-Madani:   one of several ocean models configured to simulate flows in the bay and we plan to compare the Lagrangian flow structures predicted by these models in the near future.
The Chesapeake Bay estuary consists of a main stem connected to a number of tributaries arrayed along its axis (see figure 1a). With the exception of a deeper paleochannel, running through the bay in the north-south direction with depths greater than 15m (maximum depth 63 m), the bay is quite shallow, with an average depth of 8.42 m. The rather complicated bathymetry and the 'shallow-water' nature of the flow dynamics in the bay poses a number of modelling challenges. One approach to tackle this problem is to use a topographically conformal vertical coordinate system, the socalled sigma-coordinate system (see, for example, Lynch and Werner (1991); Haney (1991) for more details). In a sigmacoordinate model the number of vertical levels in the water column is the same everywhere in the domain irrespective of the depth. Thus, in comparison to z-level models, such an approach allows for a computationally efficient resolution of the bottom boundary layer and a better numerical treatment of domains with complicated bathymetry.
The sigma-coordinate models are not free of deficiencies. Most notably, these models are susceptible to large roundoff errors which arise when calculating the pressure gradients between grid points characterized by sharp topographic changes (e.g. Haney (1991)). However, from the point of view of Lagrangian analysis, discretization-induced inconsistencies in the imposed boundary conditions, and problems associated with preservation of conserved quantities (e.g. mass and momentum) are particularly acute. Moreover, in contrast to the z-level models, the layer thickness in sigma-coordinate models varies widely between grid points (see figure 2) and the geometry of the computational sigma levels changes in time, complicating the subsequent interpolation procedures. We therefore briefly outline the main characteristics of the Quoddy model and discuss some of the limitations in using the output of this and similar models for the subsequent Lagrangian analysis. Some more technical aspects related to the implementation of the Lagrangian techniques to the discrete model output are discussed in Appendix A.
Quoddy is a 3D, finite-element, sigma-coordinate, model which solves the shallow water equation for the surface elevation and the Navier-Stokes equation, under the Bousinnesq and hydrostatic approximations, for the velocity components. It uses the Galerkin method admitting nodal quadrature on linear six-node elements so that the use of quadrature points coincident with the element nodes exactly integrates its volume. The model is 'driven' by externally imposed oceanic water level, the wind stresses across the bay and freshwater discharge from eleven rivers into the bay. The river discharge data is obtained from U.S. Geological Survey.
The governing equations for the velocity components 2 in Quoddy are given by (see Lynch and Werner (1991)) and 2 The actual equations solved by Quoddy in the Chesapeake Bay are formulated in spherical coordinates and solved in the spherical, sigma-coordinate system but we do not need this complication here. not facilitate a meaningful way of assessing transport properties of flows generated by different models. When a velocity field changes in time, its streamlines are no longer barriers to transport, even in the idealized case of an inviscid fluid flow. In the time-dependent setting the role of transport barriers is taken over by more complex structures which can be revealed via the so-called Lagrangian transport analysis rooted in the dynamical systems theory (e.g. Wiggins, 1992Wiggins, , 2005Haller, 2000;Haller and Yuan, 2000;Jones and Winkler, 2002). In quasi-turbulent flows the structures revealed by the analysis in the Lagrangian, trajectory-based framework provide a template for time-dependent transport on sufficiently large spatial scales. Many of these Lagrangian structures were successfully identified in geophysical flows as oceanic currents, fronts and eddies (e.g. Samelson and Rogerson et al., 1999;Mancho et al., 2008;Lekien et al., 2005a;Olascoaga et al., 2006;Mathur et al., 2007). These evolving structures are difficult to cross for fluid particles and they can be rather persistent. Thus, water masses separated by such Lagrangian barriers can retain their different physical and chemical properties, and nutrient content for prolonged periods of time.
Complex nonlinear models used for Eulerian ocean forecasting, including those used in modelling the dynamics in the Chesapeake Bay, are often sensitive to small variations in model attributes (e.g. different turbulence closure or slightly different boundary conditions). This can produce significantly different but equally plausible velocity fields and salinity distributions (e.g. Gross et al., 2001;Li et al., 2005). In situations when the available experimental data is sparse it is often difficult to single out the most accurate and robust model. Moreover, it is well-known in the Dynamical Systems theory that even "nearly identical" velocity fields (i.e. velocity fields which are close to each other in an appropriate norm) can produce drastically different fluid particle trajectories and affect the Lagrangian transport properties of the resulting flow 2 .
In this work we consider some issues associated with the application of Lagrangian diagnostic tools to the study of transport in the Chesapeake Bay. We certainly do not attempt to give a complete "Lagrangian picture" of the flow dynamics in the bay. Rather, our objective is to illustrate how these tools can help uncover new information on spatiotemporal flow geometry affecting Lagrangian transport in the complex geometry of the Chesapeake Bay. We also point out the limitations associated with application of such methods to shallow estuarine flows, highlighting the need for further developments. While the currently available Lagrangian tools are not ideal, we argue that techniques which emphasize the spatio-temporal geometry of considered flows, rather than their Eulerian characteristics, can provide an additional method of validation of numerical models.
In the presented analysis we mainly focus on the mouth of the bay where the most energetic dynamics takes place due to the exchange of waters between the bay and the open ocean. This area is strongly influenced by the tidal and wind activity, leading to highly intermittent flows. In particular, estuarine outflow plumes are important coastal phenomena which provide a mechanism for transporting organic and inorganic matter contained in inland water masses onto the continental shelf. Although such outflows are important to fisheries, pollution dispersion and naval operations, their variable spatio-temporal nature makes them difficult to monitor using traditional ship-based surveys, especially due to their large along-shore to cross-shore aspect ratios (e.g. Ruzecki, 1981). Remote sensing, including satellite imagery, becomes increasingly more reliable in detecting salinity fronts associated with the occurrence of such plumes (e.g. Dzwonkowski and Yan, 2005). However, validation and calibration of such methods is often problematic due to the lack of detailed in situ data. It is hoped that a synergetic approach, combining Lagrangian diagnostics of numerical models with available experimental data will lead to improvements in both the sensing techniques and the understanding of the fundamental dynamical problems.
The paper is structured as follows: in Sect. 2 we introduce the numerical model, Quoddy, which is used to obtain the velocity field in the Chesapeake Bay. We take the opportunity there to point out some common model limitations which affect the Lagrangian transport analysis in numerically generated, discrete velocity fields. Of particular importance is the nature of the boundary conditions imposed in the numerical model and the validity of the Lagrangian analysis near the boundary. In Sect. 3 we briefly introduce some currently available dynamical systems tools which can be used to study the structure of aperiodically time-dependent fluid flows in the Lagrangian framework. In Sect. 4 we discuss issues associated with applicability of the 2-D Lagrangian transport analysis in shallow, partially mixed estuaries based on simulated flows in the Chesapeake Bay. The examples illustrated in Sect. 4 point to a number of possible improvements of the currently existing techniques and provide directions for future research which we summarize in Sect. 5.

The QUODDY model
In this work we focus on the output from the Quoddy model, developed by Lynch and Werner (1991), which was adapted to the Chesapeake Bay by Gross et al. (2001). Quoddy is just one of several ocean models configured to simulate flows in the bay and we plan to compare the Lagrangian flow structures predicted by these models in the near future.
The Chesapeake Bay estuary consists of a main stem connected to a number of tributaries arrayed along its axis (see Fig. 1a). With the exception of a deeper paleochannel, running through the bay in the north-south direction with depths greater than 15 m (maximum depth 63 m), the bay is quite shallow, with an average depth of 8.42 m. The rather complicated bathymetry and the "shallow-water" nature of the flow dynamics in the bay poses a number of modelling challenges. One approach to tackle this problem is to use a topographically conformal vertical coordinate system, the so-called sigma-coordinate system (see, for example, Lynch and Werner, 1991;Haney, 1991 for more details). In a sigmacoordinate model the number of vertical levels in the water column is the same everywhere in the domain irrespective of the depth. Thus, in comparison to z-level models, such an approach allows for a computationally efficient resolution of the bottom boundary layer and a better numerical treatment of domains with complicated bathymetry.
The sigma-coordinate models are not free of deficiencies. Most notably, these models are susceptible to large roundoff errors which arise when calculating the pressure gradients between grid points characterized by sharp topographic changes (e.g. Haney, 1991). However, from the point of view of Lagrangian analysis, discretization-induced inconsistencies in the imposed boundary conditions, and problems associated with preservation of conserved quantities (e.g. mass and momentum) are particularly acute. Moreover, in contrast to the z-level models, the layer thickness in sigma-coordinate models varies widely between grid points (see Fig. 2) and the geometry of the computational sigma levels changes in time, complicating the subsequent interpolation procedures. We therefore briefly outline the main characteristics of the Quoddy model and discuss some of the limitations in using the output of this and similar models for the subsequent Lagrangian analysis. Some more technical aspects related to the implementation of the Lagrangian techniques to the discrete model output are discussed in Appendix A.
Quoddy is a 3-D, finite-element, sigma-coordinate, model which solves the shallow water equation for the surface elevation and the Navier-Stokes equation, under the Bousinnesq and hydrostatic approximations, for the velocity components. It uses the Galerkin method admitting nodal quadrature on linear six-node elements so that the use of quadrature points coincident with the element nodes exactly integrates its volume. The model is "driven" by externally imposed oceanic water level, the wind stresses across the bay and freshwater discharge from eleven rivers into the bay. The river discharge data is obtained from US Geological Survey.
The governing equations for the velocity components 3 in Quoddy are given by (see Lynch and Werner, 1991) ∂v and where v v v(x x x,t)∈IR 3 is the fluid velocity at a point x x x∈IR 3 , ξ (x,y,t) is the free surface elevation, g is the acceleration of gravity, ∇ h is the "horizontal" gradient operator, f f f is the vertical component of the Coriolis force, and R is the baroclinic pressure gradient given by where ρ 0 (x x x) is the hydrostatic fluid density and ρ (x x x,t) represents the departure from the hydrostatic density. The vertical mixing coefficient K m (x x x,t) in (1) is obtained from the second-order turbulence closure models due to Mellor & Yamada Mellor and Yamada (1982) or the Munk and Anderson (1948) closure model (see Gross et al., 2001 for more details). Note that the horizontal mixing coefficient is neglected in (1) which essentially implies an inviscid horizontal flow (see Sect. A3 for a discussion of how such boundary conditions affect trajectory calculations). The free surface elevation, ξ(x,y,t), in (1) satisfies the shallow water equations in the form (with overbars denoting the vertical average) where H (x,y,t) = ξ(x,y,t) + B(x,y) is the fluid depth, B(x,y) denotes the bottom topography andv v v h (x,y,t) represents the vertically averaged horizontal velocity components. The presence of the atmospheric shear stress and the bottom velocity function (v b ) in (5) arise from the imposed boundary conditions on the free surface and the bottom (see Lynch and Werner, 1991), in the form It is important to note that another set of "kinematic" boundary conditions, imposed on the free surface and the bottom, is implicitly utilized when deriving (4) from (2); where w denotes the vertical component of the velocity field. The implementation of Quoddy in the Chesapeake Bay solves (1) for the horizontal velocity components, (2) for the vertical velocity component, and (4)-(5) for the free surface elevation in spherical sigma coordinates. The 3-D mesh has 15 sigma levels with 9700 triangular elements in each sigma level (see Fig. 1b). The computational sigma levels are topologically identical. The time-dependent transformation from the sigma coordinates is given within the model domain, D⊂IR 3 , by where In the following discussion and in all of the figures we will refer, somewhat inaccurately, to the northward coordinate, θ, as the latitude (denoted by "lat") and to the eastward coordinate, ϕ, as longitude (denoted by "lon"). The instantaneous geometry of a sigma level in spherical coordinates is simply defined as the image of T t (·,·,σ ) given by the graph An example of instantaneous sigma level geometry in spherical coordinates is shown in Fig. 2. Note, in particular, that 1 coincides with the free surface (cyan) and 0 (brown) coincides with the bottom topography of the domain.
A reliable determination of the vertical velocity components from shallow-water ocean models is well-known to be problematic since it is inferred from the continuity Eq. (2) and not from (1). In Quoddy the continuity equation is integrated in the vertical from bottom up assuming no normal flow across the bottom topography (see 7bb). Since no other conditions are imposed, this generally leads to a nonvanishing flux across the free surface. A simple verification of this fact can be performed in the following way: the flux across any computational sigma layer, σ i , σ i ∈Z * , is proportional to the velocity component across the sigma layer, w σ i . These "vertical" velocity components in the sigma coordinates can be easily obtained from the Quoddy output at every mesh node (ϕ n ,θ n ,σ i ) as where u,v,w denote the velocity components in the spherical coordinate system, i.e. v v v = ue e e ϕ + ve e e θ + we e e r , and is the "vertical" velocity component of the mesh node (ϕ n ,θ n ,σ i ) (the mesh nodes can only move along e e e r so that v v v msh n,i = w msh n,i e e e r ).
In Fig. 3 we show an example of w σ at the free surface, the stable and unstable manifolds of those hyperbolic trajectories. Just as in the steady flow configuration, these manifolds are flow invariant, which means that trajectories cannot "leave" these manifolds or cross them. Moreover, trajectories which are confined to these manifolds and are sufficiently close to the "central" hyperbolic trajectory approach (if on stable manifold) and move away (if on unstable manifold) from the hyperbolic trajectory at exponential rates. When viewed at a sequence of observational times, the instantaneous geometry of these manifolds is given by evolving material curves (i.e. curves which are made up of the same fluid or tracer parcels) in the flow domain. Contrary to the steady 2-D flows, stable and unstable manifolds of hyperbolic trajectories in unsteady 2-D flows can intersect along multiple other hyperbolic trajectories 4 . At any fixed time, the invariant manifold structure is represented by a network of material "stable" and "unstable" curves which intersect at isolated points and form invariant lobes. Motion of these lobes provides the main mechanism for mediating Lagrangian transport between different regions and the rules governing these motions constitute the core of the so called lobe dynamics (see, for example, Wiggins, 1992).
A key to understanding the Lagrangian transport properties of a flow via the invariant manifold analysis lies in identification of those hyperbolic trajectories whose stable and unstable manifolds form homoclinic/heteroclinic tangles extending in the phase space throughout the region of interest. Such Distinguished Hyperbolic Trajectories (DHTs) play a role of organizing centers in the Lagrangian flow structure. The invariant lobes formed by their intersecting manifolds and evolving between the DHTs provide a mechanism for Lagrangian transport. In 2-D time-periodic configuration the DHTs are usually identified as low-order fixed points of the corresponding Poincaré map (low-period hyperbolic cycles). Identification of DHTs in complex aperiodically timedependent geophysical flows is generally more complicated and it may be not unique, especially in flows undergoing transitions or flows defined for finite time (see . Although the relevant theory and iterative algorithms for identifying the DHTs in arbitrary unsteady flows were developed in Ide et al. (2002);Ju et al. (2003); Mancho et al. (2004), the choice of the initial guess required in these algorithms is often not straightforward. We refer the reader to the above references for more information and only briefly note here that such a guess does not have to be a trajectory and it is often chosen to be a path of hyperbolic instantaneous stagnation points (i.e. x x x isp (t) ∈ IR n such that v v v x x x isp (t),t =0 and the eigenvalues of ∇v v v x x x isp (t),t have non-zero real parts for all t∈Ĩ ). If the DHT can be identified, its stable and unstable manifolds can be computed using algorithms developed in Mancho et al. (2004Mancho et al. ( , 2006. It is important to note here that the methods for computing invariant manifolds are not restricted to the computations of stable and unstable manifolds of hyperbolic trajectories; they are also useful in identifying the evolving geometry of any material contour in the flow. When a suitable DHT is difficult to identify, computations of evolving geometry of appropriately chosen material segments can still provide important information about the Lagrangian flow structure (e.g. Miller et al., 1997Miller et al., , 2002. We will illustrate such an application in Sect. 4.2. Another technique used in the finite-time Lagrangian transport analysis is based on the so-called finite-time Lyapunov exponent fields (FTLE). The finite-time Lyapunov exponents are measures of (finite-time) growth rates of infinitesimal perturbations of the linearized dynamics about a trajectory of the considered velocity field. Any trajectory in an autonomous, n-dimensional, continuous time dynamical system (and a smooth velocity field can be cast in such a form, possibly in an extended phase space 5 ) has n finitetime Lyapunov exponents. One of these exponents which is associated with a direction tangent to the trajectory is always zero. The values of the n−1 finite-time Lyapunov exponents associated with the remaining directions determine the finite-time stability properties of the trajectory. For each fixed t m within the considered time interval I , forward FTLE fields at t = t m , denoted here by λ T (x,y,t m ), are obtained at a point (x,y) by estimating the (locally) maximum separation at t = t m + T of nearby trajectories at t = t m (the maximum trajectory separation is related to the maximum finite-time Lyapunov exponent; for details see, for example, Shadden et al., 2005). By performing such a computation for an ordered sequence of "observation times", {t m } m∈Z , t m ∈I , one can examine the spatial evolution of the structures exhibited by the forward FTLE fields in time. Clearly, the backward FTLE fields can also be computed by reversing the direction of time. A more precise quantification of such structures have led to the notion of Lagrangian Coherent Structures (LCS) (see Haller and Yuan, 2000;Shadden et al., 2005) which are defined as the ridges of the FTLE fields. Numerous groups have computed the FTLE fields over the years in the context of fluid transport (e.g. Pierrehumbert, 1991;Pierrehumbert and Yang, 1993;von Hardenberg et al., 2000;Lermusiaux et al., 2006) and have noted that these fields appear to exhibit a great deal of structure. The FTLE computations often provide a useful diagnostic tool in the Lagrangian transport analysis and they are an attractive choice due to the relative simplicity of the required algorithms. Moreover, they often coincide with the stable and unstable manifolds of distinguished hyperbolic trajectories, at least in the neighborhood of the DHTs. However, this technique is also not free of limitations. When computing an FTLE field at some fixed t m within the considered time interval I , one can consider a whole one-parameter family, λ T (x,y,t m ) T , t m + T ∈I , of FTLE fields. It is often not obvious which length of the integration time interval T should be chosen in such computations, especially when the structure of the resulting FTLE fields varies significantly for different values of T (see  for details). Moreover, since the FTLE fields are essentially "functions" of trajectory separation, the LCS are not, a priori, material curves, and therefore not necessarily barriers to transport. (The stable and unstable manifolds of finite time hyperbolic trajectories are a priori barriers to transport since they are computed as curves of fluid particle trajectories.) However, segments of LCS are often "close" to a transport barrier in the sense that the flux across the curve may be small Shadden et al. (2005), and we use this technique extensively in Sect. 4.
An important and a rather obvious formal limitation to the two-dimensional Lagrangian transport analysis occurs if motions along the third dimension are present. Therefore, strictly speaking, no output of a three-dimensional ocean or atmospheric model should be analyzed using the 2-D Lagrangian techniques. The status quo is a combination of the lack of widely implemented methods to tackle the Lagrangian transport in three dimensions (although some methods have been developed; see Lekien et al., 2007) and the fact that many geophysical flows can be approximated, at least at some time-scales, as quasi two-dimensional flows along isentropic or isopycnal surfaces (see, for example, Samelson and Pedlosky et al., 1997;Pierrehumbert and Yang, 1993;Rogerson et al., 1999;Ngan and Shepherd, 1997;Miller et al., 2002). In the next section, we discuss problems associated with these issues in the context of flows in the Chesapeake Bay. Clearly, further developments are needed in this area in order to establish the range of applicability of the 2-D analysis and to efficiently implement the fully 3-D techniques. However, the limitations pointed out in Sect. 2 suggest that development of more reliable models might be a necessary prerequisite, especially in the coastal regions, before such assessments can be made.

A glimpse of Lagrangian flow structures in the Chesapeake Bay
We discuss here some aspects of the Lagrangian analysis of the flow generated by the Quoddy model in the Chesapeake Bay. We use these results to highlight the potential strengths of the Lagrangian techniques in uncovering the spatio-temporal flow structure of estuarine flows, as well as limitations of such analysis.

Is the 2-D Lagrangian transport analysis in the Chesapeake Bay adequate?
As already mentioned in the previous section, Lagrangian transport analysis of geophysical flows have been largely carried out under various quasi-2-D assumptions based on the flow stratification arguments. While this approach has enjoyed wide-spread applications to deep ocean basins, one might question applications of such techniques to shallow coastal regions or estuaries. Important limitations in the fully 3-D Lagrangian transport analysis of simulated geophysical flows stems from the fact that the vertical velocities produced by the shallow water ocean models are notoriously unreliable, since they are usually determined diagnostically from the continuity equation, rather than from the 3-D Navier-Stokes equation. As already discussed in Sect. 2, the Quoddy model is no exception in this regard. Therefore, while one might question the validity of the two-dimensional analysis of the 3-D Quoddy output, far-reaching conclusions based on the analysis of the vertical velocities generated the model are equally prone to criticism. Nevertheless, we felt it necessary to investigate, based on the 3-D velocity field generated by Quoddy, whether or not any two-dimensional Lagrangian analysis can be justified in the Chesapeake Bay. Computations of the 3-D trajectories in the bay were performed using the interpolation technique described in Appendix A2. This technique is capable of processing the Quoddy 3-D velocity data, computed and stored at the nodes of a time-dependent, semi-structured grid.
In Fig. 4 we show an example of the 3-D trajectory geometry at the mouth of the Chesapeake Bay. In order to better visualize the trajectory geometry, two different projections of the same configuration are shown in the two panels with the point of view indicated by the red arrow in the insets. The trajectories are shown after a 12-h run and they originate in the second sigma level, 2 , along a line which "blocks" the entrance to the bay in the top view (green dots); the magenta dots denote the end locations. The horizontal to vertical axes aspect ratio is deliberately exaggerated in order to aid the visualization. Note the significant difference in the depth range explored by different trajectories which seems to be correlated with the bottom topography (see also Fig. 6).
It is tempting to assume that, at least on sufficiently short time scales, either the computational sigma levels (see Fig. 2) or constant depth levels are approximately flow invariant. Such an invariance would allow, in principle, for quasi-2-D analysis. We illustrate a computational test of such a hypothesis in Fig. 5 where we compare the geometry of trajectories evolving according to Quoddy's fully 3-D velocity field (green) and trajectories confined to evolve within the second sigma layer (red). Both types of trajectories evolve from initial locations within the second sigma level, 2 . The "confinement" of trajectories to the particular sigma level is achieved by simply computing trajectories of the 2-D velocity field in the sigma level and determining the depth from (10). Despite the fact that the average ratio of the magnitude of the vertical to the horizontal velocity components is small (of the order of 10 −4 ), the 3-D trajectories cross a number of sigma levels (and depths) within hours (see also Fig. 6) and the quasi-2-D assumptions are difficult to justify, except for the free surface given by 1 . In Fig. 5b we compare the same trajectories as in Fig. 5a but projected onto the longitude-latitude subspace. It is worth noting that, although there are significant discrepancies in the end locations of the computed trajectories, the overall correlation between the 3-D and the "sigma-constrained" trajectories is surprisingly good.
The information about the three-dimensionality of the flow generated by Quoddy in the Chesapeake Bay can be cast in the form of "depth-range" maps. The computational procedure for generation of such maps consists of the following three simple steps: (1) Choose a grid of initial conditions located at a fixed depth (depth-range maps) or within a chosen sigma level (sigma-range maps), (2) Compute trajectories from the grid of initial conditions for a fixed period of time T from t=t 0 to t=t 0 +T , (3) Identify the maximum and minimum depth (or sigma layer) visited by the respective trajectory during the evolution and assign the depth/sigma range to the corresponding initial condition.
The scalar maps resulting from such a procedure are parameterized by the starting time t 0 , the integration period T and the reference depth at which the grid of initial conditions is located. Thus, for example, for any fixed t 0 and the reference depth one could generate a one-parameter family of depth range maps parameterized by T . This fact obviously leads to non-uniqueness of the computed results, similarly to the case of FTLE maps (see Sect. 3). Nevertheless, such a technique can serve as a useful diagnostic tool when assessing the validity of 2-D analysis in complex geophysical flows in the absence of more rigorous theoretical arguments.
We show two examples of such maps in Fig. 6 which are computed at the mouth of the bay for T =12 h. Figure 6a shows a sigma-range map which indicates the range of sigma levels visited by trajectories originating within the second sigma level. In this example, a significant proportion of the 3-D trajectories in the bay mouth area remains within the first four layers, although there are islands of deeper protrusions which are correlated with the bottom topography. Our extensive calculations (not shown here) suggest that this type of behavior is quite representative of the Lagrangian behavior at the mouth of the bay (see also Figs. 4 and 5). Figure 6b shows a depth-range map which indicates the vertical extent of trajectories originating at the depth of 1 m and evolving for 12 h. In summary, the 2-D Lagrangian transport analysis in the Chesapeake Bay, at least according to the Quoddy model output, does not seem adequate within the main body of the fluid. The shallow bathymetry of the bay combined with the complicated bottom topography not only makes the flow modelling a difficult task but also requires fully 3-D Lagrangian analysis of the resulting flows. We note that there exist techniques extending the invariant manifold and FTLE methods to aperiodically time-dependent, threedimensional fluid flows (see Lekien et al., 2007). However, bearing in mind the issues associated with a non-vanishing flux across the discretized model boundary and the free surface (see Sect. 2), it seems that development of more reliable models is necessary be-fore embarking on efficient implementation of the 3-D Lagrangian methods to numerically generated velocity fields.

Lagrangian geometry of the surface flow in Chesapeake Bay
In this section we use the 2-D Lagrangian tools to study the spatio-temporal geometry of the surface flow in the Chesapeake Bay. Such an approach can be used in transport analysis of buoyant Lagrangian tracer patches (e.g. oil spills, algal blooms) which can be thought of as non-diffusing, passive tracers on the considered spatio-temporal scales. Based on the results presented in the preceding subsection, the 2-D flow analysis in the top computational sigma layer seems to

Fig. 5.
Comparison of the geometry of trajectories evolving according to the fully 3-D Quoddy velocity field (green) and trajectories confined to evolve within the second sigma layer (red). Both types of trajectories evolve from initial locations within the second sigma level, 2 . Clearly, the sigma levels (with the exception of 1 which coincides with the free surface) bear little relation to material surfaces. (b) Comparison of the same trajectories as in (a) but projected onto the (lon, lat) coordinates. Although there are significant discrepancies in the end locations of the computed trajectories, the overall correlation between the 3-D and the sigma confined trajectories is surprisingly good.

M. Branicki and R.
Malek-Madani: Lagrangian structure of flows in the Chesapeake Bay 11 able nature of these events (see, e.g. Ruzecki (1981); Guo and Valle-Levinson (2007)) makes it difficult to monitor using mooring arrays and shipboard surveys which infer the evolution of the salinity front associated with the freshwater plume from density or salinity measurements (Xu et al. (2002); Johnson et al. (2001)). A typical plume event begins with an injection of fresh estuarine waters onto the nearshore continental shelf (Chao and Boicourt (1986)). On the northern hemisphere, the outflow initially makes an anticyclonic turn back towards the right-hand coast, forming a pronounced southward jet. In large estuaries, such as the Chesapeake Bay, such outflows can be considerable and travel along-shore for tens of kilometers so that the Earth's rotation becomes a non-negligible factor in their dynamics (see (Li et al., 2005, Figure 6) for some interesting simulations with and without the Coriolis force which confirm this point).
The ocean color data (OC) from remote satellite sensing is increasingly used to track the evolving structure of estuarial plumes (see, e.g. Dzwonkowski and Yan (2005), see also flow plume. However, due to unavailability of in situ data which would allow for a direct verification of such correlations, these findings remain speculative. The Lagrangian analysis of numerical models of the bay can be useful in calibrating the remote sensing data and help obtain a coherent picture of the plume evolution.
The right column in figure 7 shows the backward FTLE fields (see §3) computed over a 24 hour period and at the same 'observation' times as the corresponding panels in the left column; the instantaneous geometry of the material contour (red curve) considered earlier is superimposed on each of the FTLE fields for comparison. Recall that the ridges (yellow) of the backward FTLE fields are, in principle, hallmarks of strong coherent, attracting structures in the flow. There is a clear correlation between the instantaneous location of the considered material contour and a subset of the strongest ridges in the computed FTLE fields, suggesting that the evolved material segment is indeed shadowing a developing front. be the only case where such an approximation is adequate since, due to the imposed boundary conditions, the free surface in shallow water ocean models is, in principle, a material surface (see (7b)). We note that the Lagrangian transport in surface flows, obtained both from shallow water ocean models and high-frequency radar measurements, was analyzed numerous times in the past using the 2-D assumptions (see, for example, Lekien et al., 2005b;Olascoaga et al., 2006Olascoaga et al., , 2008 and led to some insightful results.
When carrying out such an analysis it is worth bearing in mind certain model-dependent limitations. As already mentioned in Sects. 2 and 4.1, the non-vanishing flux across the top sigma level in Quoddy violates, strictly speaking, the assumed "impermeability" assumptions imposed on the free surface. The flux is small compared to the typical vertical velocities in the Quoddy output but it may nevertheless cause problems in trajectory computations. We recall, however, that the erroneous flux arises as a consequence of the "bottom-up" integration of the continuity equation Nonlin. Processes Geophys., 17, 149-168, 2010 www.nonlin-processes-geophys.net/17/149/2010/ which, as in other shallow water models, is enslaved to the Navier-Stokes equation for the horizontal velocity components. Thus, in our surface flow computations, we confine the computed trajectories to the free surface by determining the vertical trajectory coordinates from (10) rather than from a direct integration of the vertical velocities. The strong dependence of surface flows on the imposed wind stresses requires reliable and detailed 6 experimental data which are often less than ideal. Clearly, the robustness and accuracy of the surface flows simulated by Quoddy cannot be investigated based on the output of a single model with prescribed wind forcing. Therefore, we regard the following analysis as a prelude to a wider study comparing the Lagrangian structures in the output of multiple models. We first consider the surface flow near the mouth of the Chesapeake Bay which, due to its proximity to the open ocean, is characterized by the most energetic dynamics. In the examples shown in Fig. 7 we are interested in the exchange of surface waters across the mouth of the bay. We first analyze the evolution of a material segment which initially "blocks" the entrance to the bay. In the left column of Fig. 7 we show the instantaneous geometry of such a segment (thick red) at different times during its evolution. The timedependent geometry of the material contour is computed using the same techniques as those exploited in the numerical approximation of stable and unstable manifolds introduced in Mancho et al. (2004Mancho et al. ( , 2006, except that the initial segment is not seeded at a hyperbolic trajectory. The computation is performed in an adaptive fashion, employing insertion and redistribution of points on the contour, representing the initial conditions for subsequent trajectories, in order to maintain the required resolution. Consequently, the evolving segments are not "made-up" of the same trajectories at each "observational" time instant. A few examples of the trajectory geometry are shown in blue in the left column panels; the initial locations of these trajectories are marked by the blue dots and the end locations are denoted by the green dots. Clearly, one can obtain a much clearer picture of the tracer redistribution by following the material contour (red) instead of relying on the "spaghetti" plots of trajectories (blue) of individual fluid parcels. An interesting aspect of this computations is associated with the development of a series of "bulges" along the south-western Atlantic coast, below the mouth of the bay, which is reminiscent of the estuarine freshwater outflow plume. The highly variable nature of these events (see, e.g. Ruzecki, 1981;Guo and Valle-Levinson, 2007) makes it difficult to monitor using mooring arrays and shipboard surveys which infer the evolution of the salinity front associated with the freshwater plume from density or salinity measurements (Xu et al., 2002;Johnson et al., 2001). A typical plume event begins with an injection of fresh estuarine waters onto the near-shore continental shelf (Chao and Boicourt, 1986).
On the Northern Hemisphere, the outflow initially makes an anticyclonic turn back towards the right-hand coast, forming a pronounced southward jet. In large estuaries, such as the Chesapeake Bay, these outflows can be considerable and they can travel along-shore for tens of kilometers so that the Earth's rotation becomes a non-negligible factor in their dynamics (see Li et al., 2005, Fig. 6) for some interesting simulations with and without the Coriolis force which confirm this point). The ocean color data (OC) from remote satellite sensing is increasingly used to track the evolving structure of estuarial plumes (see, e.g. Dzwonkowski and Yan, 2005, see also http://oceancolor.gsfc.nasa.gov). Contrary to the more traditional in situ techniques which usually monitor salinity gradients, these techniques rely on differences in optical properties of the inland waters carrying large amounts of dissolved organic materials. As reported in Dzwonkowski and Yan (2005) there is compelling evidence that there is a strong correlation between patterns in the OC observations near the mouth of the Chesapeake Bay and the outflow plume. However, due to unavailability of in situ data which would allow for a direct verification of such correlations, these findings remain speculative. The Lagrangian analysis of numerical models of the bay can be useful in calibrating the remote sensing data and help obtain a coherent picture of the plume evolution.
The right column in Fig. 7 shows the backward FTLE fields (see Sect. 3) computed over a 24-h period and at the same "observation" times as the corresponding panels in the left column; the instantaneous geometry of the material contour (red curve) considered earlier is superimposed on each of the FTLE fields for comparison. Recall that the ridges (yellow) of the backward FTLE fields are, in principle, hallmarks of strong coherent, attracting structures in the flow. There is a clear correlation between the instantaneous location of the considered material contour and a subset of the strongest ridges in the computed FTLE fields, suggesting that the evolved material segment is indeed shadowing a developing front.
The example shown in Fig. 7 illustrates an important symbiotic relationship between the invariant manifold methods and the FLTE methods in the Lagrangian analysis of complex geophysical flows. Note that since the material contour in Fig. 7 was not associated with a stable/unstable manifold of a DHT, any conclusions about transport properties of the surface flow based solely on the evolution of such a segment would be unfounded. Similarly, the analysis based solely on the FLTE maps does not provide sufficient insight into the Lagrangian flow structure in the surface flow. The strong ridges present in the FTLE maps at the mouth of the Chesapeake are rather short (in the arc length sense) and their distribution and connectivity is sensitive to the time interval chosen for the FTLE computation. Moreover, the network of strong ridges revealed in the computed FTLE fields can be quite complex, making it impossible to identify the most important coherent structures responsible for organizing the Fig. 7. Trajectories (blue, left column) and backward FTLE fields (right column) computed for a surface flow at the mouth of the Chesapeake Bay. The red curves denote the geometry of material segments evolving in the flow from the straight line 'blocking' the mouth after 12h (first row), 31h (middle row), and 55h (bottom row). Note that the material segments shown develop a series of bulges along the western shore which are reminiscent of the freshwater outflow plume. In this flow, the location and connectivity of the ridges in the FTLE fields are sensitive to the integration time T (see also §3 and Branicki and Wiggins (2010)) Fig. 7. Trajectories (blue, left column) and backward FTLE fields (right column) computed for a surface flow at the mouth of the Chesapeake Bay. The red curves denote the geometry of material segments evolving in the flow from the straight line "blocking" the mouth after 12 h (first row), 31 h (middle row), and 55 h (bottom row). Note that the material segments shown develop a series of bulges along the western shore which are reminiscent of the freshwater outflow plume. In this flow, the location and connectivity of the ridges in the FTLE fields are sensitive to the integration time T (see also Sect. 3 and Branicki and Wiggins, 2010). Lagrangian transport at the bay's mouth. However, if two approaches are combined, the correlation between some of the strong FTLE ridges and the geometry of the material contour help create a more complete picture.
Determination of the Distinguished Hyperbolic Trajectories and their stable and unstable manifolds in the flow allows, in principle, for identification of dominant structures mediating Lagrangian transport via the turnstile mechanism (see, for example, Samelson and Wiggins, 1992Wiggins, , 2005Mancho et al., 2006, 2008 andSect. 3). However, this proves to be a difficult task in the case of the studied surface flow. This is due to the fact that the paths of the instantaneous hyperbolic stagnation points (ISPs) in the surface flow in the Chesapeake Bay, which are often used as the initial guess in an the iterative DHT-finding algorithm (see Ju et al., 2003;Ide et al., 2002), frequently bifurcate in time and are too short-lived to allow for a computation of sufficiently long segments of the relevant hyperbolic trajectories. The initial guesses for the location of the DHTs in the flow could be constructed, in principle, by following in time intersections between ridges of forward and backward FTLE maps. Such intersections mark instantaneous locations of finite-time hyperbolic structures in the flow and represent likely instantaneous locations of DHTs (see . However, the complexity of the "ridge network" at the mouth of the Chesapeake Bay and the fact that the ridge intersections are often not preserved during the flow evolution, render techniques based on such a strategy prohibitive at this stage. We finish this section with an example of Lagrangian structures which can be uncovered in the surface flow in the midsection of the bay where the flow is less energetic. In Fig. 8a-b we show examples of a backward FTLE field (a) and forward FTLE field (b) computed in the surface flow with T =24 h. The strong ridges (red) in the backward FTLE field (a) are indicative of attracting structures and the ridges in the forward FTLE fields in (b) are indicative of repelling structures in the flow field. In contrast to the situation encountered in the surface flow at the mouth of the bay, and probably owing to a longer characteristic time scales of the flow in this region, the location of the ridges in the FTLE fields does not vary in time as rapidly as at the bay's mouth. In Fig. 8c we show a ridge network generated by superimposing the attracting ridges (red) obtained form the backward FTLE field (a), and the repelling ridges (blue) obtained from the forward FTLE field (b). In this "quieter" part of the bay, the paths of hyperbolic instantaneous stagnation points are longer-lived than at the mouth of the bay and it is possible to compute longer segments of distinguished hyperbolic trajectories. In such cases we found that the stable and unstable manifolds of these DHTs tend to align well with the respective ridges of the FTLE fields. We show one such case in Fig. 8d where the dashed green curve represents the instantaneous geometry of the unstable manifold (grown for about 20 h), and the dashed blue curve is the stable manifold of a hyperbolic trajectory (black dot). It is, perhaps, interesting to note that, although the intersections of strong ridges of forward and backward FTLE fields are likely instantaneous locations of the DHTs, the ridge intersection visible just below the identified DHTs does not produce another DHT. Instead, trajectories originating in the neighborhood of this intersection are first attracted towards the identified DHT and swept away form this region along its unstable (green) manifold.

Concluding remarks and future work
It has been long recognized that the Chesapeake Bay, as many other estuaries, is a difficult system to model due to its complex bathymetry, interaction with the open ocean and the presence of high salinity gradients. Due to the strategic location of the bay and its commercial importance, there have been many efforts aimed at improving the Eulerian model predictions of the 3-D flow structure in the bay. While the accuracy of the developed models and their parameterizations have been steadily improving, relatively little is known about the nature of Lagrangian transport in large estuaries like the Chesapeake Bay. Moreover, it is not known how good models like Quoddy, POM (see Guo and Valle-Levinson, 2007), or ROMS (see Li et al., 2005) are at predicting the Lagrangian transport properties in the bay. The presented work was aimed at initiating the study of such issues, and at providing a critical overview of existing methods in the hope of stimulating further research and improvements.
We used in the analysis the output of Quoddy, a shallow water, finite-element model adopted to the bay by Gross et al. (2001). Lagrangian analysis of the Qouddy output highlighted a number of problematic issues which require further improvements, both in the (Eulerian) modelling and in the techniques for Lagrangian "postprocessing". The discretization process of the imposed boundary conditions on the rigid walls and on the free surface makes the trajectory computations a challenging task, since the domain boundaries are not necessarily "impermeable", as desired. We developed a simple method capable of remedying this deficiency in Sect. A3 without compromising the model data. However, the broader issue of adequacy of the boundary conditions used in such models and the validity of Lagrangian analysis in coastal regions remains unclear. It is important to be aware of the potential for an appearance of artifacts in the Lagrangian flow structure in such situations. The currently available 2-D techniques for studying the Lagrangian transport, based on computation of the stable and unstable manifolds of hyperbolic trajectories or exploiting the 2-D finite-time Lyapunov exponent fields, are of limited use in the case of partially mixed shallow-water estuarine flows. We showed in Sect. 4 that the 2-D Lagrangian transport analysis in the Chesapeake Bay, at least according to the Quoddy output, is not adequate within the main body of the bay's waters and that it should be restricted only to the analysis of the surface flow. Similarly, paths of hyperbolic instantaneous stagnation points in the bay are generally too short-lived to serve as a useful initial guess in the DHT finding algorithm (see Ide et al., 2002;Mancho et al., 2006 and Sect. 3). In the situations where relatively short segments of DHTs can be identified the stable and unstable manifolds of these DHTs tend to align well with the respective ridges of the FTLE fields, as shown in (d); the dashed green curve represents the instantaneous geometry of the unstable manifold, and the dashed blue curve is the stable manifold of a hyperbolic trajectory whose instantaneous location is denoted by the black dot.
Clearly, a development and efficient implementation of threedimensional tools is needed in order to gain a better insight into the spatio-temporal structure of flows in estuaries like the Chesapeake Bay. We note that there exist non-trivial extensions of the the invariant manifold and FTLE methods to aperiodically time-dependent, three-dimensional fluid flows (see Lekien et al., 2007) but the theory of lobe dynamics in this higher-dimensional setting has not been developed to date. Bearing in mind the issues associated with a non-vanishing fluxes across the model domain boundaries, it seems that development of more reliable models is necessary before embarking on the implementation of the 3-D Lagrangian methods to numerically generated velocity fields. The examples of the 2-D Lagrangian analysis of the surface flow in the Chesapeake Bay discussed in Sect. 4.2 showed their usefulness in capturing the spatiotemporal variability of the freshwater outflow events, setting the stage for a systematic validation and comparison of predictions obtained from different models. However, it is also worth remembering that these techniques are not free of deficiencies. Computations of finite-time Lyapunov fields can reveal a network of ridges delineating regions of rapid trajectory separation (in both forward and backward time). These structures are, however, generally non-unique (see  and the detected ridges are often too short (in the arc length sense) for a meaningful transport analysis within the framework of lobe dynamics. Computation of stable and unstable manifolds of relevant hyperbolic trajectories in realistic flows also poses a number of challenges. These are mainly associated with difficulties with constructing the initial guesses for the locations of relevant (i.e. distinguished) hyperbolic trajectories which locally organize the flow structure. The paths of hyperbolic instantaneous stagnation points used frequently in such a case are often too short-lived for this purpose in the analyzed surface flow in the Chesapeake Bay.
As discussed in  based on a series of idealized examples, the intersections of the ridges of the forward and the backward FTLE fields are likely instantaneous locations of the Distinguished Hyperbolic Trajectories (see also Sect. 4.2). A robust, automated technique for detection and continuation in time of such ridge intersections is however lacking. It would be highly desirable to develop a synergetic approach to the analysis of Lagrangian transport by combining the various Lagrangian tools instead of artificially separating them. The deterministic Lagrangian approach to transport in oceanic flows is certainly not the ultimate framework capable of solving all the issues associated with transport in geophysical flows which, after all, involve diffusive as well as advective processes operating on all scales. It is also possible that the great mesoscale complexity of Lagrangian structures transpiring in realistic flows will necessitate a drastically different approach. We believe, however, the ideas and concepts arising through an interplay of Lagrangian analysis, rooted in the dynamical systems theory, and improved Eulerian modelling of complex geophysical flows will lead to an advancement in our understanding of these complex processes.

Interpolation of the QUODDY data
The output velocity data of the Quoddy model is stored at 9700 nodes of a semi-structured grid which is composed of 15 sigma levels in the vertical and unstructured triangular meshes within each of the time-dependent sigma levels (see (10) and Fig. 2). Thus, the computation of fluid particle trajectories requires the use of a suitable interpolation technique. In what follows we first outline the 2-D interpolation procedure in each sigma level. We show subsequently how this technique can be extended to the interpolation in the 3-D physical domain between the time-dependent sigma levels.

A1 2-D interpolation within a sigma level
We employ here a local interpolation technique which is based on a procedure introduced in Franke and Nielson (1980), referred to as the modified Shephard's method. This method was proposed to overcome drawbacks of the original interpolation scheme introduced by Shephard (1968). Further extensions of this method, employing the radial basis functions to interpolate multivariate data sets, can be found in Lazzaro and Montefusco (2002). The modified Shephard's method has many advantages such as numerical efficiency, good reproduction quality and stability. Moreover, the presented method is local, i.e. the interpolated value at a given point depends on a subset of nodal values of nodes which are close to the considered point. Consequently, the method can be rather easily parallelized. Appropriate enhancements are, however, beyond the scope of this paper.
We note that a different 2-D interpolation technique was used in Brasher (2005). The global interpolation employed there used compactly-supported Wendland radial basis functions (see Wendland, 1998).
Below we describe the main steps of our interpolation technique for a scalar field in IR 2 . Interpolation of vector fields in IR 2 is performed in a similar way, independently for each component.
Consider a set of nodes M = {m 1 ,...,m n } ∈ IR 2 connected to each other by edges encoded in the sparse and symmetric connectivity matrix M. For a triangular mesh the connectivity matrix satisfies Moreover, if M ij = 1 and M ik = 1, then necessarily M kj = 1. We will refer to the pair (M, M) as the triangular mesh. Assume now that a function F(x x x, t) : IR 2 × IR → IR is given at some fixed t at the nodes of the mesh, i.e. we know the set of values {F(m i )} m i ∈M . We want to interpolate F at a point x x x = (x, y), given its K neighboring nodes m x 1 , . . . , m x K in the mesh (M, M). We will denote this set of K neighbors by K ⊂ M.
The set of neighbors K is chosen based on the connectivity matrix rather than the Euclidean distances from the considered node. When the domain boundary has a complex, non-convex geometry (as in the case of the Chesapeake Bay), such a procedure avoids choosing near-boundary nodes which lie within a given radius from the considered point in the Euclidean metric but are, say, on the other side of an island. We denote the Cartesian coordinates of the neighbors as m k x = (x k , y k ) ∈ IR 2 and the points corresponding to these mesh nodes in the graph of F as (x k , y k , z k ) ∈ IR 3 , i.e. z k = F(m k ).
Following Franke and Nielson (1980), the interpolating function F (x, y) is given by where the nodal functions, Q k , satisfy Q k (x k , y k ) = z k and fit the values of the remaining points in the least-squares sense. Here, we choose the nodal functions in the form of the n-th order bivariate polynomials given by (In order to simplify the notation, we will skip the second subscript n in (A3) and write Q k in (A2) instead of Q k,n .) The weight functionsW k are defined as Assume now that a function F (x x x,t): IR 2 ×IR→IR is given at some fixed t at the nodes of the mesh, i.e. we know the set of values {F (m i )} m i ∈M . We want to interpolate F at a point x x x = (x,y), given its K neighboring nodes m x 1 ,...,m x K in the mesh (M,M). We will denote this set of K neighbors by K⊂M.
The set of neighbors K is chosen based on the connectivity matrix rather than the Euclidean distances from the considered node. When the domain boundary has a complex, non-convex geometry (as in the case of the Chesapeake Bay), such a procedure avoids choosing near-boundary nodes which lie within a given radius from the considered point in the Euclidean metric but are, say, on the other side of an island. We denote the Cartesian coordinates of the neighbors as m k x = (x k ,y k )∈IR 2 and the points corresponding to these mesh nodes in the graph of F as (x k ,y k ,z k )∈IR 3 , i.e. z k = F(m k ).
Following Franke and Nielson (1980), the interpolating function F (x,y) is given by F (x,y) = K k=1W k (x,y)Q k (x,y), where the nodal functions, Q k , satisfy Q k (x k ,y k ) = z k and fit the values of the remaining points in the least-squares sense.
Here, we choose the nodal functions in the form of the n-th order bivariate polynomials given by Q k (x,y) = z k + n i=0 c i (x − x k ) i (y − y k ) n−i .
(In order to simplify the notation, we will skip the second subscript n in (A3) and write Q k in (A2) instead of Q k,n .) The weight functionsW k are defined as where, denoting r k (x,y) = (x − x k ) 2 + (y − y k ) 2 , W k (x,y) = r W k − r k + r W k r k α , α > 2, and r W k − r k + = r W k − r k , if r k < r W k , 0, if r k r W k .  Even in such a case the interpolated velocity is, in general, not going to be tangent to the boundary. Note that imposing a constraint in the interpolation method which would force the interpolated velocities on each boundary segment to be tangent to that boundary segment introduces further problems at the boundary nodes. Of course, one would ideally like know the continuous boundary geometry, rather than its discretized version. In the absence of such an information, one is forced to either somehow determine a C 1 reconstruction of the 'true boundary' so that the model velocities at the boundary are everywhere tangent to it, or to modify the velocities at the boundary nodes. In the computations presented in the next section we try to remedy this problem by adding a 'ghost' boundary whose nodes are slightly offset outwards with respect to the Quoddy boundary nodes. We then set the velocity to zero at the nodes of the secondary boundary. This procedure obviously leads to a generation of a sharp boundary layer (outside the original computational domain) within which the flow is forced to adjust to the added no-slip boundary. We do not attach significance to the results obtained within the coastal 'buffer zone'. Although it seems very desirable for models to incorporate more realistic boundary conditions, it is rather difficult to imagine an implementation which would realistically treat the complicated shore line of the Chesapeake Bay and the river inflows. Indeed, similar questions can be raised in the case of other ocean models in coastal regions. Therefore, one must be cautious when interpreting the results of Lagrangian analysis in coastal regions. This is especially important in situations associated with identifying points of separation of the certain Lagrangian structures from the coast.