High-resolution nested model simulations of the climatological circulation

Abstract. As part of the Mediterranean Forecasting System Pilot Project (MFSPP) we have implemented a high-resolution (2 km horizontal grid, 30 sigma levels) version of the Princeton Ocean Model for the southeastern corner of the Mediterranean Sea. The domain extends 200 km offshore and includes the continental shelf and slope, and part of the open sea. The model is nested in an intermediate resolution (5.5 km grid) model that covers the entire Levantine, Ionian, and Aegean Sea. The nesting is one way so that velocity, temperature, and salinity along the boundaries are interpolated from the relevant intermediate model variables. An integral constraint is applied so that the net mass flux across the open boundaries is identical to the net flux in the intermediate model. The model is integrated for three perpetual years with surface forcing specified from monthly mean climatological wind stress and heat fluxes. The model is stable and spins up within the first year to produce a repeating seasonal cycle throughout the three-year integration period. While there is some internal variability evident in the results, it is clear that, due to the relatively small domain, the results are strongly influenced by the imposed lateral boundary conditions. The results closely follow the simulation of the intermediate model. The main improvement is in the simulation over the narrow shelf region, which is not adequately resolved by the coarser grid model. Comparisons with direct current measurements over the shelf and slope show reasonable agreement despite the limitations of the climatological forcing. The model correctly simulates the direction and the typical speeds of the flow over the shelf and slope, but has difficulty properly re-producing the seasonal cycle in the speed. Key words. Oceanography: general (continental shelf processes; numerical modelling; ocean prediction)


Introduction
The continental shelf at the southeastern corner of the Mediterranean Levantine Basin has a fairly simple bathymetric structure with depth contours aligned roughly parallel to the coastline.The shelf along the entire coast of Israel is relatively narrow as compared to the shelf of the Egyptian coast, where the effects of the sediment drift from the Nile littoral cell are most prominent.The shelf width varies from approximately 60 km off the coast of Sinai to 20 km wide in the south of Israel and narrows to about 10 km in the north.The bathymetry contours are shown in Fig. 1.Direct current measurements conducted at various locations along the shelf and slope of Israel over a ten-year period reveal a predominantly bathymetry-following, northward flow, which is part of the basin wide, cyclonic current that follows the coast of the Levantine Basin.In contrast to this, hydrographic data collected during the Marine Climate, MC, cruises (Hecht et al., 1988), the Physical Oceanography of the Eastern Mediterranean, POEM, cruises (Brenner, 1989;The POEM Group, 1992;Oszoy et al., 1993) and the EDDY cruises (Brenner et al., 1991;Brenner, 1993) all indicate that the circulation in the offshore waters (deeper than 1000 m) in the central latitudes of the Levantine Basin is predominantly anticyclonic due to the presence of the recurrent Shikmona gyre system.Properly simulating this transition between the coastal to the deep-water circulation regimes will be a crucial test of the model's capabilities.
In addition to the extensive data collection programs that have been conducted throughout the Mediterranean over the past 15-20 years, various efforts have focused on modelling the circulation of the entire Mediterranean, as well as various sub-basins, with the goal of understanding the general phenomenology and specific process studies of the region.The basin wide models, such as Roussenov et al. (1995), Zavatarelli and Mellor (1995), and Wu and Haines (1998), were low resolution (roughly 25 km horizontal grid spacing and 19 vertical levels or less), designed to study the climatolog- ical general circulation of the entire Mediterranean.Other studies focused on particular process studies and/or the subbasin circulation.For example, Korres et al. (2000a, b) used the same model as Roussenov et al. (1995) but increased the number of vertical levels to 31 and improved the surface forcing methodology, in order to study the response of the general circulation to interannual atmospheric variability.Wu et al. (2000) increased the horizontal resolution of the model of Wu and Haines (1998) to 0.125 • (roughly 13 km), in order to study the process of deep-water formation in the eastern Mediterranean but used a highly simplified surface forcing parameterization.Finally, Lascaratos and Nittis (1998) used an intermediate resolution, eddy-resolving (5.5 km grid) model of the Levantine and southern Aegean Seas to study the process of intermediate water formation.While this list of previous modelling studies is by no means exhaustive, it does provide a fair indication of the wide interest in studying the circulation of the Mediterranean.
Most recently, the Mediterranean Forecasting System Pilot Project (MFSPP) has addressed the development of a prototype operational forecasting system (MFSPP, 2000;Pinardi et al., 2003).The project has included all of the necessary aspects of an operational forecasting system, including data acquisition, data assimilation, near real-time forecasts of the basin scale circulation, and the development of various nested intermediate regional and localized, high-resolution shelf models.It is within this context that we have devel-oped a high-resolution model for the southeastern corner of the Levantine Basin, as shown in Fig. 1.Our main goal here is to implement the model, establish the nesting procedure, and test the model in a perpetual year, climatological simulation.In Sect. 2 we describe the model, nesting procedure, and surface forcing, while in Sect. 3 we present the results from a multiyear simulation.In Sect. 4 we discuss our results within the context of present-day knowledge of the circulation in this region and in Sect. 5 we present a brief summary and conclusions.

Numerical model description
The model that has been chosen for our work is the Princeton Ocean Model (POM).The model has been used extensively to simulate the coastal circulation in various regions around the world, including the Mediterranean Sea (e.g.Zavatarelli and Mellor, 1995), and as such has been described elsewhere.The most detailed, basic description of the model is given in Blumberg and Mellor (1987).Here we only provide a brief description of the model.POM is a three-dimensional, timedependent model based on the primitive equations.It consists of prognostic equations for the two components of the horizontal momentum, potential temperature, salinity, and the free surface.Three diagnostic equations consisting of the hydrostatic equation, an equation of state, and the vertical velocity, which is derived from the mass continuity equation, supplement these.In addition, POM contains an imbedded higher order turbulence closure scheme to simulate the vertical mixing (Mellor and Yamada, 1982).This adds two additional prognostic equations to the model -one for the turbulent kinetic energy and one for the turbulence macroscale.
The equations are solved using finite differencing in both space and time.In the horizontal, the variables are located on the grid following the second order accurate, quadratic conserving Arakawa C scheme.Potential temperature, salinity, and the free surface are located in the center of each grid box, while the horizontal velocity components are staggered one-half grid box in their respective coordinate directions.
In the vertical, POM uses a terrain following sigma coordinate (Phillips, 1957) in which the water column at each grid point is normalized by the total water depth.Thus, sigma varies from 0 at the surface to -1 at the bottom.The sigma coordinate facilitates including the effects of topography and simplifies the representation of the bottom boundary condition.However, this is done at the expense of introducing some potential inaccuracies in the computation of the horizontal pressure gradient, which consists of the difference of two large terms in the sigma system.To minimize this "sigma problem", we follow the procedure of Mellor et al. (1994) in which the domain mean profile of the density is subtracted before computing the horizontal pressure gradient, and the bathymetry is smoothed based on a criterion to avoid hydrostatic inconsistency.
The time integration is done with a split explicit scheme in which the barotropic and baroclinic modes are integrated separately with a leapfrog scheme but with different time steps.Since POM is a free surface model, the barotropic time step is more restrictive, since it is limited by the speed of the external gravity wave.The external step is chosen to be smaller than the CFL requirement, while the internal mode time step is 30-60 times the external step.
The domain covered by our shelf model is shown in Fig. 1.It covers the region from 33. 5-35.43 • E and 31.15-33.7 • N. The horizontal grid resolution is 1/48  (Korres and Lascaratos, 2003).Similarly, in the vertical, we use the same 30 sigma levels as in the intermediate model.This choice of open boundaries and sigma levels is done to facilitate the interpolation between the intermediate model grid and our high-resolution grid.The sigma levels are more closely spaced in the upper ocean so that the top tenth of the water column contains seven levels.The topmost layer is only 0.0025 thick.Below σ = 0.1, the levels are equally spaced with a layer thickness of 0.04.The bathymetry for the model is bilinearly interpolated from the one-minute resolution (1/60 • ) data available from the U.S. Navy.For depths below 100 m the bathymetry is smoothed, as noted above, while for shallower depths we retain the original bathymetry.The minimum depth in our model is 6 m.The model bathymetry, as well as a portion of the horizontal grid, is shown in Fig. 1.

Surface forcing
In contrast to some of the previously mentioned modelling studies in which the surface forcing consisted of wind stress and buoyancy forcing parameterized as the relaxation of the surface temperature and salinity to climatological values, here we force the model with wind stress, heat fluxes, and fresh water flux.The wind stress and heat flux components (shortwave, longwave, latent and sensible) were computed with standard bulk formulae using the ECMWF reanalysis data for the years 1979-1993.The six-hourly values of the total cooling flux (longwave + latent + sensible) and wind stress were computed and then averaged to form 10-day means, while monthly means of the shortwave flux were computed from an astronomical formula with climatological cloud cover.The use of bulk formulae for computing the fluxes can lead to variations of as much as 50% depending upon the details of the formulation.One advantage of the Mediterranean region is that an independent constraint on the net surface heat flux is given by the net heat flux through the Straits of Gibraltar.In order to achieve a balanced annual heat budget, various flux formulations were tested and compared by Garrett et al. (1993) and by Castellari et al. (1998).The final calibrated heat fluxes used here were the same as those used to drive the intermediate model and are described in detail by Korres and Lascaratos (2003).These consisted of values that were computed with the bulk formulae and then adjusted by adding a weak relaxation to climatological surface temperatures, in order to produce the correct heat budget.These fluxes were used in our model with no further adjustment.Finally, for the fresh water flux we computed the difference between the evaporation (computed from the latent heat flux) and the precipitation taken from the monthly mean climatological values of Jaeger (1976).This flux was adjusted by adding a weak relaxation to the climatological surface salinity.All of the above were bilinearly interpolated to the model grid and linearly interpolated in time to the model time step.

Lateral boundary conditions and nesting
As part of MFSPP, our high-resolution shelf model, as with the shelf models for other regions, is the third in a sequence of one-way nested models.The system begins with the full Mediterranean general circulation model (Pinardi et al., 2003), OGCM, which is run on a grid of 1/8 • (roughly 13 km).Within the OGCM, various regional or intermediate models are nested with a grid resolution of around 5 km.In our case the intermediate model of the Levantine, Aegean, and Ionian basins (Korres and Lascaratos, 2003), ALERMO, which is also based on POM, is run with a resolution of 1/20 • .Finally, our model is nested in ALERMO.In each case, the nesting ratio is around 2.5:1.The lateral boundary conditions for our model were extracted from the third year of a perpetual year forcing ALERMO simulation.Tenday mean ALERMO fields were available at 10-day intervals.The ALERMO values were bilinearly interpolated onto our model grid.For grid points located in water depth shallower than 50 m (the minimum depth of ALERMO), the required values were extrapolated.The variables needed to nest our model are: potential temperature, salinity, the two total (barotropic + baroclinic) velocity components (normal and tangential), and the two barotropic velocity components.The free surface is not nested since this would over specify the system.In the following are the details of the nesting procedure for the western and northern open boundaries.
For potential temperature and salinity, at inflow points, the values are specified from ALERMO, while at outflow points the values are extrapolated from the interior solution using an upstream advection equation of the form where T represents the potential temperature or salinity, t is time, U is the normal component of velocity, and x is the normal coordinate.The boundary condition is applied to each vertical layer individually.
Both the normal and tangential components of the total velocity were specified by bilinearly interpolating directly from ALERMO.Thus, where U and V are the normal and tangential components, respectively, and the subscripts POM and COARSE refer to the shelf model and the intermediate model (ALERMO), respectively.As with the potential temperature and salinity, these boundary conditions are applied to each vertical layer individually.
For the barotropic velocity, when the normal component was directly interpolated from ALERMO, we found that the model blew up after a very short time.To circumvent this problem, an extra term representing a weak external gravity wave radiation condition was added to allow the excess energy to slowly leak out of the limited domain.Thus, the normal component was given by where the overbar stands for the barotropic (depth averaged) velocity, provides the relative direction of the outward propagating wave, according to =1 on the northern boundary and -1 on the western boundary, g is gravity, H is the water depth, and ζ is the free surface.The tangential component of the barotropic velocity on the boundary was specified directly from ALERMO so that Finally, in order to ensure long-term stability and to maintain compatibility with the coarse resolution model, an integral constraint was imposed on the normal velocity across each of the open boundaries.The purpose of this constraint is to adjust the net mass flux across the boundary so that the high-resolution flux exactly matches the original coarse resolution flux (i.e.before the spatial interpolation).While there are various ways to impose such a constraint (i.e.different choices of the weighting function), we have chosen the simplest form in which the cross sectional area mean flux across the entire boundary is preserved.Thus, the integral constraint for the total normal velocity takes the form

C H
where l is the horizontal coordinate along the open boundary with end points l 1 and l 2 , σ is the vertical (sigma) coordinate, the superscript corr refers to the adjusted or corrected high-resolution variable, and all other notation is as defined above.The integrals are computed over the coarse (C) and high-resolution (H) cross sections, respectively.The adjustment reflects the difference between the cross sectional area of the lateral boundary in the coarse and high-resolution grids which arises from the spatial interpolation.The adjustment of the high-resolution normal velocity is, therefore, computed as follows.First, we define the operator which is the cross sectional area average of U (using the respective coarse or high-resolution cross section) from which it follows that the integral constraint of Eq. ( 5) can be written as T COARSE = T corr POM .The total transport across the boundary as computed by the high-resolution model before the adjustment is given by H which in general will not be the same as T COARSE due to the slight differences between the boundary cross sectional area before and after the spatial interpolation.Thus, in order to preserve the original coarse grid transport, we correct the high-resolution velocity according to

H
The high-resolution cross sectional area integral of Eq. ( 7) is exactly equal to Eq. ( 5) and thus, we find that the corrected high-resolution velocity will indeed preserve the original total transport across the boundary.

Results
The initial conditions for the results presented here were taken from 10 January of the third year of the perpetual year ALERMO simulation.The required fields were bilinearly interpolated to our model grid and extrapolated for grid points with depths less than 50 m, as noted above.The model was then integrated for 37 months with the perpetual year surface and lateral boundary forcing described in the previous section.Unless otherwise noted, the results presented are from the third year of the integration.We begin with Fig. 2, which shows the evolution of the domain mean potential temperature, salinity, and kinetic energy throughout the entire simulation.These three integral quantities measure the model spin up and indicate whether or not the model has reached an equilibrium state in the form of a repeating annual cycle.The potential temperature and salinity fields spin up rather quickly with no noticeable initial shock to the system that would require a major adjustment.The potential temperature clearly reflects the seasonal cycle of the upper layers with a maximum in late July to early August and a minimum in late February to early March.There is some internal variability, as shown by the lower maximum in the second summer as compared to the higher maxima in years one and three.This is despite the perpetual year surface and boundary forcing.The other feature that stands out is a small relative maximum on the curve in early January of  each year.This is due to the repeating cycle of the perpetual year and the fact that there is apparently a slight discontinuity between the forcing in December (the end of the forcing year) and the forcing in January (the beginning of cycle).
The salinity curve is analogous to the potential temperature.It shows little evidence of any initial shock and spins up to a repeating annual cycle rather quickly.The annual cycle is controlled primarily by processes in the upper thermocline where the seasonal variations are most apparent.The integrated salinity is lowest in the summer, despite the fact that the thin surface layer of Levantine Surface Water (LSW) has the highest salinity values in the entire basin.The low summer integrated salinity is associated with the strongest inflow of relatively fresh Atlantic Water (AW).The maximum integrated salinity occurs in late winter when the upper 200 m of the water column are well mixed and when the relatively saline Levantine Intermediate Water (LIW) forms.As with the potential temperature curve, a discontinuity appears in January due to the reset of the perpetual forcing year, although here it is much less apparent.
Of the three integral quantities, the kinetic energy exhibits the most pronounced adjustment during the first few months of the simulation.While the kinetic energy curve is much noisier than the thermodynamic curves, the dynamics also shows a clear repeating seasonal cycle with a maximum in February and a minimum in July-August.During the initial adjustment process, the maximum of over 9×10 14 J near the beginning of the simulation drops to less than one-half that value (4×10 14 J) in the subsequent winters.The higher frequency variability, with a time scale of about two months, is most likely due to the dynamically active mesoscale field (Robinson et al., 1987;Hecht et al., 1988).Based on all three quantities we can see that the third year of the simulation is certainly representative of the fully adjusted fields.
In Fig. 3 we show snapshots of the shelf model potential temperature and velocity fields at 30 m for 20 February of year 2 (Fig. 3a) and for year 3 (Fig. 3b).Thus, these two plots represent the flow 13.5 and 25.5 months into the integration.In both cases the three dominant features that recur are: (i) an intense jet forming a large cyclonic meander which enters in the southwest and exits near the center of the western boundary (we will refer to this as J1); (ii) a large anticyclonic eddy in the northwest corner of the domain (AE1), and (iii) a small anticyclonic eddy in the southeast corner of the domain (AE2).The meander formed by J1 is the largest feature present and occupies nearly half of the domain.It advects a sharply defined tongue of cold water that is 1-1.5 • cooler than the surrounding environment.Speeds in the jet are as high as 40-45 cm/s.In the northeast, there is a clear tendency to form secondary meanders which lead to the formation of smaller scale eddies with a typical radius of 15-20 km.For comparison, the local Rossby radius for the first baroclinic mode is on the order of 10-12 km (Hecht et al., 1988).The bifurcation of the jet and eddy formation is much more apparent in the year 2 plot, although even in year 3 we can see the early stages of formation.AE1 in the northwest corner is a large, warm core eddy, which clearly originates outside of our model domain.It is controlled primarily by the lateral boundary conditions and, therefore, appears quite similar in both years, although the southeastern flank of this eddy is clearly formed by jet J1.This eddy can be identified as the Shikmona gyre that is known to be a large, recurrent anticyclonic feature, typically centered around 34 • E, 34 • N (Hecht et al., 1988;Brenner et al., 1991;The POEM Group, 1992;Brenner, 1993).Finally, AE3 is a small, warm core, anticyclonic eddy that slides along the southern edge of the eastward flowing portion of jet J1.Its size is comparable to the local Rossby radius of deformation.It is centered over the continental slope, although the eastern edge enters the shelf region.While the overall structure is fairly similar in both years, the differences noted above clearly indicate that there is internal variability and differences from year to year, even though the external forcing is repetitive.
Figure 4 shows snapshots of the potential temperature and velocity at 30 m on 20 August of year 2 (Fig. 4a) and year 3 (Fig. 4b), which are 20.5 and 32.5 months, respectively, into the integration.A large cyclonic feature formed by jet J1 that enters the domain from the southwest corner dominates the circulation.Maximum speeds in the jet are still on the order of 40 cm/s.In contrast to the winter case, it now covers almost the entire domain and forms a large, well-defined eddy.Also, the jet closely follows the coast with a clear bifurcation near 33.2 • N, where most of the flow continues northward and exits the domain while only a small part turns to the west.The other noticeable difference from the winter is the nearly complete absence of any small mesoscale eddies.This further confirms our assertion above that the summer minimum in the mean kinetic energy is associated with the weakening or lack of an energetic mesoscale field.The temperature contrast between the shelf and the open sea is much more pronounced than in winter, with the shelf being as much as 5-6 • warmer over most of the domain and nearly 9 • warmer in the southwest corner.The part of the jet that separates from the coast and turns westward forms a narrow stream of warm water that separates the cool water in the center of the eddy from the cool water to the north.The pictures are similar for both years, although in year 3 the eddy center and associated westward jet along its northern flank are displaced 45 km to the south.
Next, we consider the deeper thermocline circulation.In Fig. 5 we show the potential temperature and velocity on 20 February at 300 m for year 2 (Fig. 5a) and year 3 (Fig. 5b).Many of the prominent near-surface features also appear at this depth, including the large cyclonic meander associated with jet J1 in the center of the domain; the large anticyclonic eddy, AE1 (Shikmona gyre) located in the northwest corner (along the open boundaries); and the smaller anticyclonic eddy, AE2 in the southeast corner of the domain.Jet J1 is not as tightly defined along the northern flank of the meander as in the near-surface layer and as might be expected, the current speeds are weaker than at 30 m, with typical values of 20 cm/s or less.Horizontal temperature contrasts are much smaller than at 30 m and do not exceed 1.5 • .Also, the large cyclonic meander is clearly defined by a cold core, while the signal of the cold tongue advected by J1 near the surface is absent at this depth.
Figure 6 shows the potential temperature and velocity on 20 August at 300 m for year 2 (Fig. 6a) and year 3 (Fig. 6b).The dominant feature in both years is still the cold core, cyclonic eddy that was seen at 30 m, although now it appears as an independent eddy since the coastal jet is not present at this depth.The currents are weaker than at 30 m, with maximum values of 20 cm/s or less.Also, at 30 m the entire open sea part of the domain was covered with relatively cool water, except for the stream of warm water that was advected westward by the jet flowing along 33 • N.Here there is a clear separation of the cool water in the center of the eddy covering the southern half of the domain from the warmer water in the north, although maximum temperature differences are only around 1 • or less.
We conclude this section with Fig. 7, which shows snapshots of the salinity on 20 August at 100 m for year 2 (Fig. 7a) and year 3 (Fig. 7b).This depth was chosen since it is located near the bottom of the layer of inflowing, relatively fresh Atlantic Water (AW), which is most noticeable in summer (Hecht et al., 1988).The cold core, cyclonic eddy that dominated the potential temperature fields (Figs. 4 and 6) also appears here with a core of relatively high salinity.The westward flowing jet near 33 • N also forms a clear separation between relatively fresh water in the south and saline water in the north.The gradient across this front is as much as 0.1 psu.In both years the core of the eddy is well defined.In year 2 it is somewhat more saline, while in year 3 it is located further to the south (as shown in the potential temperature plots).The other noticeable difference between the two years is that in year 3 the core is a single isolated mass of saline water, while in year 2 the core of saline water is surrounded by alternating streaks of less saline and more saline, which appear to spiral into the center.The freshest water with salinity of about 38.7 psu, which flows into the domain from the southwest corner, advances along the coast, reaching almost as far as the northern boundary.This is most likely AW, which is known to appear off the coast of Israel (Hecht et al., 1988).We do not show the corresponding winter maps since in late February the water column is typically mixed to a depth of 150-200 m and, therefore, a clear signal of AW is absent.

Discussion
In the previous section we have presented some selected results from a simulation with our high-resolution model of the southeastern corner of the Levantine Basin, forced with cli-matological surface fluxes and lateral boundary conditions taken from an intermediate resolution model which was, in turn, nested in a coarse resolution full Mediterranean model.The model was stable throughout the three-year integration and spun up to a repeating annual cycle within the first year.Furthermore, the simulated values of temperature and salinity were found to be consistent with the observed climatological hydrography of this region (Hecht et al., 1988).
While it would be most instructive to fully assess the model's performance, a detailed quantitative analysis of the model's skill is beyond the scope of this study, primarily due to the lack of relevant verification data sets.Nevertheless, we can perform an additional qualitative check of the simulation in two ways.First, we can compare the results to the corresponding results of the intermediate model, which is not influenced by lateral boundary conditions in our region of interest.Second, we can compare the simulated velocities to observed direct current measurements that were conducted at various times and locations along the continental shelf and slope of Israel over a 10-year period.

Comparison with the intermediate model
In order to compare our results to the corresponding ALERMO results, we present Figs.8-12 which correspond to the snapshot fields shown above in Figs.3-7.Since the ALERMO results represent ten-day means rather than instantaneous fields, for a proper comparison in panel (a) of each of the figures we show the ten-day means from year 3 of our model, while panel (b) shows the corresponding ALERMO fields.Upon comparing our ten-day means to the corresponding instantaneous fields shown previously (panel (b) of Figs.3-7), we find that although some of the details are lost in the averaging process, nevertheless, the major features, both large and small, still appear.The overall appearance of the high-resolution and the intermediate model results are quite similar in all cases, although there are some noticeable and important differences.Bearing in mind that the 20 February fields (Figs.8a and 10a for potential temperature at 30 and 300 m, respectively) are 25.5 months into our integration, while the 20 August fields (Figs.9a and 11a for potential temperature at 30 and 300 m; Fig. 12a for salinity at 100 m) are 31.5 months into the simulation, the overall similarities indicate that the surface and lateral boundary forcing strongly control and constrain the circulation in this region.Nevertheless, the internal differences as noted above indicate that the model still has a certain amount of dynamic freedom.The fields along our two open boundaries appear to transition smoothly into the ALERMO fields with no indications of instabilities or spurious boundary effects.This confirms, a posteriori, the success of the nesting procedure.
While the similarities between the results of the highresolution and intermediate models are encouraging, it is also important to examine the differences between the two simulations, in order to fully understand the benefits of using the nested, high-resolution model.In the 30 m February potential temperature and velocity fields (Figs.8a and b), the two   main differences between the models are seen in the small warm core, anticyclonic eddy in the southeast (referred to as AE2 above) and in the westward flowing jet along the northern flank of the large, cyclonic meander in the center of the domain.Eddy AE2 appears in both models; however, in ALERMO it appears about 40 km further to the northeast and is followed by a small cyclonic eddy, which seems to have recently formed from a meander of the jet.The inability of the high-resolution model to reproduce this cyclonic eddy is most likely due to the suppression of small meander and eddy formation close to an open boundary, which is constrained by the specified lateral boundary conditions.The second, and perhaps more significant difference, is the westward flowing jet along 32.7 • N. In ALERMO it appears as a broad current, which forms the northern side of the large meander and which smoothly merges with the westward flow along the south side of the Shikmona gyre.In the highresolution model, however, the northeast corner of the large meander appears to be a dynamically active region where smaller meanders and pinched-off eddies are formed.One such cyclonic eddy can be seen at 33.7 • N, 34.5 • E in Fig. 8a.Also, as a result of this eddy formation, the westward flowing jet is stronger and more tightly defined and separated from the broad flow along the south side of the Shikmona gyre.
These clear and important differences in the interior of the domain of the high-resolution model are due to the finer grid spacing and the ability of the higher resolution simulation to better simulate smaller scale features.
In the 30 m potential temperature field in August, once again, the general picture is similar in both models, with the circulation dominated by the single large, cold core, cyclonic eddy in the center of the domain, although in the ALERMO case (Fig. 9b) the center of the eddy is located 50 km further to the northeast as compared to the high-resolution results (Fig. 9a).The major difference between the two models, however, is in the coastal current.In our model the strong jet is confined to the continental shelf zone, while in ALERMO it is much broader and extends out over the slope and into the open sea.In both models the coastal current begins to bifurcate around 33 • N but with very different results.In our model part of the current turns westward to form the warm jet separating the cold core eddy to the south and the colder water to the north, but most of the jet continues to follow the coast and flow northward out of the domain in the northeast corner.In ALERMO, the current separates from the coast and also turns westward, forming a warm and intense warm jet flowing along the northern edge of the eddy.Only a very small part of the current continues northward along the coast.This importance difference is also due to the higher resolution in our model and is probably closely linked to the representation of the bathymetry in the two different grids.At the deeper levels (300 m potential temperature and velocity in February, Figs.10a and b; and August, Figs.11a and b) and 100 m (salinity in August, Figs.12a and b), the results are quite similar in both models, although ALERMO has some difficulty capturing the northernmost branch of the slope current.
We conclude this subsection with Fig. 13, which shows a comparison between the free surface height (cm) in February for our model and ALERMO (Figs. 13a and b, respectively) and for August (Figs.13c and d).The free surface height represents the combined effects of the barotropic mode and the vertically integrated baroclinic circulation.It confirms the results shown in the previous figures in the sense that the most consistent and pronounced feature in the circulation is the large cyclonic meander in the center of the domain, which takes on the form of a well-defined eddy in summer.It also confirms the inability of ALERMO to spin off smaller cyclonic eddies at the northeast corner of the meander in winter (compare Figs. 13a and b), while in summer the jet along the northern edge of the eddy is more intense in ALERMO, leading to a larger and deeper cyclonic center (Figs.13c and d).Also, in the August ALERMO field (Fig. 13d), the seaward extension of the shelf-slope free surface gradient near 33.3 • N is associated with the separation of the current from the coast that does not occur in the highresolution model.Finally, except for the centers of anticyclonic eddies, the highest free surface heights appear along the southeast coast of our domain as a result of the wind set up associated with the strong and persistent northwesterly component of the wind in this region.

Comparison with direct current measurements
As we noted above, a quantitative verification of our simulations is beyond the scope of this project due to the lack of appropriate data sets.However, we do have a rather extensive set of direct current measurements that were conducted at various locations along the continental shelf and slope of Israel between 1987-1996 (Rosentraub et al., 2002).Due to the climatological nature of our simulation, it is not possible to conduct a direct comparison with the observations.However, we can compare some of the general properties of the simulated and observed currents to further assess our model's performance.In Fig. 14 we show the monthly mean alongshore (approximately northward) velocity from two current meters that were located on the inner shelf (dash-dot line) for several years.They were 18 m below the surface in water depth of 26 m.One was located in the south (31.7 • N) and one was located near the center (32.5 • N) of the coast of Israel.The observed mean currents are directed northward, except in September when it is close to zero.Typical mean velocities are 5-10 cm/s, and there is a clear bimodal seasonal signal, with the strongest mean northward flow in February and July and the weakest northward flow in May and September.At the closest corresponding grid points and depths, the model (solid line) correctly simulates mean northward flow throughout the year, with typical velocities also of 5-10 cm/s and a bimodal seasonal cycle but with a lag of about two months as compared to the observed currents.While we have no immediate explanation for this time lag, it has been suggested that the flux correction that was implicitly incorporated in the surface heat fluxes has a similar time lag (Pinardi, private communication).Another possible reason is the lack of spatial details of the local wind field in the climatological wind stress (Rosentraub, private communication).Both of these possibilities will be explored in future work.
In Fig. 15 we show the mean monthly alongshore (roughly northward) currents from two midshelf current meters (dashdot line) that were located 30-40 m below the surface in water depth of 120 m also off the southern and central coast of Israel.The model simulated currents from the closest corresponding grid points and depths are shown by the solid line.Here the mean observed currents are also directed northward throughout the year, with typical values of 5-10 cm/s and with a pronounced maximum of nearly 30 cm/s in June-July.The model successfully simulates mean northward flow throughout the year, with comparable velocities and a single pronounced maximum of 30 cm/s.However, here the peak occurs in September, so once again we see that the simulated seasonal cycle is off by about two months.
Finally, in Fig. 16 we show the currents from a single current meter that was located over the middle of the continental slope at 40-50 m, below the surface in water depth of 500 m off the central coast of Israel.Here the mean currents, also directed northward, are noticeably stronger than at the shallow and the midshelf stations.We should note that the observed southward flow in September is not representative, since it consisted of a record of only two weeks of measurements from a single year.Typical observed and simulated values are 10-20 cm/s with summer and winter maxima of 25-30 cm/s.However, once again, we notice that the simulated maxima lag the observations by about two months.

Summary and conclusions
In this work our main goal was to present the implementation of a high-resolution, doubly nested model for the southeastern corner of the Levantine Basin of the Mediterranean Sea.We were primarily concerned with the spin up and the long-term stability of the model, as well as the performance of the nesting scheme.Therefore, we focused on a multiyear, climatological simulation forced by perpetual year, monthly mean surface heat fluxes and wind stress, and by ten-day mean lateral boundary conditions taken from a similarly forced coarser grid model.In addition to spatial interpolation from the coarse grid, the lateral boundary conditions were subjected to an integral constraint in which the original, intermediate model, cross sectional area weighted total mass flux across each boundary is preserved by the high-resolution model.A comparison with the known hydrography and phenomenology of this region showed that the model is able to successfully reproduce many of the observed features in this region.Similarly, a comparison with the intermediate model simulation showed that the high-resolution model was capable of reasonably reproducing the larger (domain) scale circulation of the former with no adverse boundary effects.The requisite information was able to pass freely into and out of the high-resolution domain across the boundaries.At the same time our model was also capable of simulating some of the smaller scale dynamics not captured by the coarser resolution intermediate model, especially in the shelf region which is not adequately resolved by the former model due to its minimum depth cutoff of 50 m.Finally, a comparison with observed, long-term, direct current measurements conducted at several locations along the shelf and slope region proved to be quite encouraging, despite the limitations of the climatological forcing used in our simulation.The model correctly simulated the typical monthly mean observed current speeds and directions throughout the year, although the seasonal cycle in the speed lags the observed cycle by about two months.Since this model will eventually form part of a planned operational Mediterranean forecasting system, the next obvious step will be to repeat these experiments but using observed synoptic atmospheric forcing rather than climatological heat fluxes and wind stress.

Fig. 1 .
Fig. 1.Model domain showing the bathymetry and a portion of the grid.

Fig. 2 .
Fig. 2. Evolution of domain means of (a) potential temperature, (b) salinity, and (c) kinetic energy during the three years of integration.

Fig. 13 .
Fig. 13.Simulated free surface height: (a) 10-day average for 10-20 February of the third year, (b) 10-day average for 10-20 February from the intermediate model, (c) 10-day average for 10-20 August of the third year, and (d) 10day average for 10-20 August from the intermediate model.

Fig. 14 .
Fig. 14.Seasonal cycle of the year 3 simulated (solid line) and multiyear observed (dash-dot line) monthly mean currents at the shallow stations (current meter depth 18 m; water depth 26 m).