Tidal dynamics of the Terminos Lagoon, Mexico: observations and 3D numerical modelling

The tidal circulation patterns in the Terminos Lagoon were studied based on the analysis of 1 year of measurements and numerical simulations using a baroclinic 3D hydrodynamic model, the MARS3D. A gauging network was installed consisting of six self-recording pressure–temperature sensors, a tide gauge station and two current profilers, with pressure and temperature sensors moored in the main lagoon inlets. Model simulations were validated against current and sea level observations and were used to analyse the circulation patterns caused by the tidal forcing. The numerical model was forced with eight harmonic components, four diurnal (K1, O1, P1, Q1) and four semi-diurnal (M2, S2, N2, K2), extracted from the TPX0.7 database. The tidal patterns in the study area vary from mixed, mainly diurnal in the two main inlets of the lagoon, to diurnal in its interior. The tidal residual circulation inside the lagoon is dominated by a cyclonic gyre. The results indicate a net flux from the southwest Ciudad del Carmen inlet (CdC) towards the northeast Puerto Real inlet (PtR) along the southern side of the lagoon and the opposite in the northern side. The results indicate two areas of strong currents in the vicinity of the inlets and weak currents inside the lagoon. The area of strong currents in the vicinity of the CdC inlet is larger than that observed in the PtR inlet. Nevertheless, the current analysis indicates that the highest current speeds, which can reach a magnitude of 1.9 m s−1, occurred in PtR. A further analysis of the tide distortion in the inlets revealed that both passages are ebb dominated.


Introduction
The Terminos Lagoon (TL) is the largest coastal lagoon in Mexico  and is located in one of the largest Mexican natural reservations. It has an area of 7,050 km 2 and is surrounded by one of the biggest wetlands of the region. The TL is located at the south Gulf of México (Fig. 1). Kjerfve and Magill (1989) classify it as a restricted lagoon in their comparative oceanography of coastal lagoons. It has a maximum length of 75 km oriented along the shore and a width of 35 km, in the cross shelf direction. The mean lagoon depth is 3.5 m, and it is separated from the Gulf of México by two islands: Carmen Island, which is 38 km long and 2.5 km wide, and Aguada Island, which is 40 km long and 900 m wide. Both of these islands are Holocene calcareoussand barrier islands (Gutierrez-Estrada and Castro-del-Río 1988;Yáñez 1963). The lagoon is connected with the Gulf of Mexico (GoM) in three points. The first point is the Ciudad del Carmen (CdC) inlet, which is located in the southwest, 3.4 km wide, and has a maximum depth of 19 m, with an average depth of 4.5 m. The second point is the Puerto Real (PtR) inlet, which is 3.2 km wide with a maximum depth of 12 m and an average depth of 3.5 m. The third point is the Sabancuy inlet, which is located 22.6 km away from the main lagoon and connected to it by a narrow and shallow channel. In this inlet, the observed flows are significantly smaller than those observed in the other two inlets. Therefore, this inlet is neglected in this study. The inlets are deeper than the rest of the lagoon. Although the TL has been previously studied, France and Mexico launched the bi-national project Joint Environmental Study of Terminos Lagoon (JEST) in 2009, with the aim of performing a multidisciplinary study looking toward the main physical, chemical and biological components and environmental processes. One of the main objectives of this project is the full understanding of the lagoon hydrodynamics and its interaction with the GoM.
The TL climate can be characterised by three distinct climatic seasons (Yáñez -Arancibiia A and Day 1982): the rainy season (Jun-Sep) when the trade winds weaken and the stable atmosphere allows strong convection processes that generate rainfall; the Nortes season (Oct-Mar), when highpressure fronts transit over the GoM towards the Isthmus of Tehuantepec, influencing the TL region with wind speeds higher than 12 m s −1 and the dry season (Mar-May), during which the sea-land breezes blow with a strong diurnal signal that in April-May reaches its major mean speed of 4.0 m s −1 . The mean direction alternates southeasterly to northwesterly (Taylor -Espinosa 2009;Gille et al. 2005). In the months of spring when the maximal temperatures occur, the wind speeds decay and the rivers' discharges are at their minimum (David and Kjerfve 1998); for these months, the TL circulation is strongly influenced by the astronomical tide. During this season, the lagoon reports the saltiest values of the year, and because of the diurnal wind signal, the water column homogenises with salty seawater (David and Kjerfve 1998;Robadue et al. 2004).
Through the different seasons of the year, the circulation varies due to the magnitude of the multiple driving forces, including river discharges, wind, shelf circulation and tides. The astronomical tides are particularly important in the TL circulation, which makes it important to understand the circulation conditions without considering river discharges and wind forcing. Knowledge of the circulation generated only by astronomical tides is vital to the understanding of the physical processes involved in the occurrence of strong currents at the inlets. The knowledge of the residential time during tidal circulation (i.e. to evaluate the risk of prolonged periods of polluted water inside the TL in case of an oil spill) is also crucial.
The shelf circulation over the Eastern Campeche Bank (ECB) of the GoM is characterised by a southwesterly current that flows, almost year-round, parallel to the coast. The current only reverses from September to November, when the current over the southwestern GoM shelf reverses because of the wind seasonality (Dubranna et al. 2011;Zavala-Hidalgo et al. 2003). There is evidence of coastally trapped waves (CTP) travelling from the north and passing along the southeast coast of the Campeche Bay, with a period of 6.1 to 10.7 days. The CTP are most likely generated in the northern Tamaulipas-Veracruz shelf by wind (Dubranna et al. 2011).
Tidal propagation in the GoM displays, for the main diurnal and semi-diurnal components, an amphidromic point north of the Yucatan Peninsula, and all components have an anticlockwise propagation (Kantha 2005;Egbert and Erofeeva 2002). The diurnal components dominate the GoM; however, due to the extension of the continental shelf in the west part of Florida, the west side of the Mississippi delta, and in the southeast of the Campeche Bay (SECB), the semi-diurnal components are amplified (Clarke and Battisti 1980), generating a semi-diurnal behaviour in the first two places mentioned and a mixed, mainly diurnal, behaviour in the SECB. In the SECB, most components are amplified northeastwards.
The literature dedicated to the circulation inside the TL is not abundant. It is based upon scant measurements of currents, water level, temperature and salinity (Mancilla-Peraza and Vargas-Flores 1980;Candela 1983;Jensen et al. 1989;David and Kjerfve 1998;Espinal-González 2002;Kjerfve et al. 1988). Those studies are based on the measurements of periods smaller than 3 months and involve limited spatial coverage. Semi-permanent southwestward circulation has been observed; the water masses preferentially enter into the lagoon via the PtR inlet and leave it through the CdC inlet (Mancilla-Peraza and Vargas-Flores 1980;Jensen et al. 1989;David and Kjerfve 1998). Nevertheless, during the winter season, during the occurrence of northers, the lagoon responds to the external circulation with increased water level and modification of the mean circulation (David and Kjerfve 1998;Candela 1983;Espinal-González 2002). At the end of the rainy season, when river discharges increase, the lagoon circulation changes because of the freshwater input. Model results indicate freshwater export through the CdC inlet, but some models also indicate export through the PtR inlet, reversing the southwestward current (Dressler 1982;Jensen et al. 1989;Kjerfve et al. 1988;David and Kjerfve 1998;Espinal-González 2002).
The first formal study of tidal circulation inside and outside the TL was conducted by Mancilla-Peraza and Vargas-Flores (1980); this study reports a flow entering from PtR and going out at CdC. The authors found a peak at the time frequency of the O 1 harmonic component using the current spectral analysis. The report of Grivel-Piña et al. (1982) addresses tide measurements inside the lagoon. The authors suggest the presence of four tidal harmonic constituents (K 1 , O 1 , M 2 , S 2 ), obtained at five different points inside the lagoon. Candela (1983) reports that, during neap tide periods, the semi-diurnal tide components' amplitude reduction decreases and the phase propagation suffers a smaller delay inside the TL. The diurnal harmonic component suffers a time delay of 27 min, arriving first at the CdC inlet. Also, his results show an increase in the amplitude of the harmonic components and a delay in the phases inside the TL. David and Kjerfve (1998) concludes that the TL presents a predominantly diurnal mixed tide signal in both inlets, and they also report a harmonic analysis in three points, presenting three diurnal harmonic components (Q 1 , O 1 , K 1 ), three semidiurnal components (N 2 , M 2 , S 2 ) and M 4 . They also conclude that the three major diurnal and semi-diurnal constituents progress from the northeast to the southwest in the GoM with high waters first occurring at PtR then at CdC. Their measurements show an amplitude decay inside the TL and a phase delay. According to their research, the tide inside the lagoon is diurnal. They also present a current analysis at the two main inlets, in both cases, at a depth 3.5 m below the surface, presenting the current harmonic components in the form of ellipses with the current intensity of the major and minor axes, as well as the major axis orientation and the component phase. The analysis was performed for the same harmonic components previously mentioned.
There are four studies about the numerical modelling of tides in TL; the first one was done by Graham et al. (1981), who uses a 2D hydrodynamic numerical model, and the tidal boundary conditions were estimated from an analysis of the National Ocean Survey (NOS) Tide Tables (1979). The data are tidal range and phase. The phase delay and amplitude differences between CdC and PtR are solved, increasing 171 s to the phase in CdC and increasing 0.60 m the amplitude. Graham's investigation concludes that, during light wind conditions, the TL circulation could be strongly determined by tidal forcing. In 1982, Dressler constructed a 2D model, which was forced using the four harmonic components reported by Grivel-Piña et al. (1982). The modelling used an artificial bathymetry, closing and opening the inlets artificially. The harmonic components were individually propagated. Kjerfve et al. (1988) used the four harmonic components reported by Grivel-Piña et al. (1982). Their 2D model included river and wind forcing. The researchers concluded that the tides are mixed, but primarily diurnal, in the inlets and the tides at CdC lead those in the PtR inlet by 1 h.
Finally, Espinal-González (2002) used the harmonic component reported in David and Kjerfve (1998) and reported a tidal circulation starting from PtR and later in CdC. The study focused on the stability of the TL inlets, considering waves, currents and sediment properties.
From the literature, we can conclude that the tidal forcing in the generation of currents is an important factor, particularly in the two principal inlets. It is possible to conclude that, in the TL, the tides are mixed (mainly diurnal) at the inlets and diurnal inside. Nevertheless, consensus regarding the directivity of tide propagation in the inlet has not been reached. Some researchers advocate a southwest to northeast direction, while others report opposite findings. Additionally, the amplification or reduction of harmonic components within the lagoon is not clear, and evaluation of circulation within the TL has relied predominantly on artificial bathymetry.
Using data from CO-OPS of NOAA and from the Servicio Mareografico de Mexico (Mexican Hydrographic Service), from several tide gauge stations along the GoM coast, stretching from Louisiana to Campeche, it is possible to observe amplitude amplification and a southeast phase delay in several harmonic components. This information, compared with the TPX0.7 database, presents similar features for the GoM.
The aim of this article is to clarify the tidal circulation in the TL and to determine the dynamics of the circulation features generated by tides at the inlets and inside the lagoon, considering a double inlet system. Measured currents, sea level observations and a numerical simulation model (MARS3D) are used to analyse the response to tidal forcing. The study is based on the analysis of 1 year of measurements taken with the gauging network described in Sections 2 and 3. A high-resolution numerical simulation of the TL was developed as part of this research. Its calibration and validation with the measurements are described in Section 4. In Section 5, cotidal charts and some features of amplitude decay and phase delay at the inlets are described. In Section 6, we explain some of the main tidal circulation patterns using current ellipses over the entire domain. In Section 7, tidal currents are discussed. In Section 8, Lagrangian residual velocity and the residual current patterns are analysed. Section 9 presents a discussion section, and finally, a summary and the conclusions are presented in Section 10.

Measurements
Six pressure-temperature self-recording sensors (TP) were strategically placed in the lagoon. The network was deployed from Jan 27, 2010, to Jan 28, 2011. Current profilers (CP) were installed bimonthly. One was deployed in the CdC inlet and performed measurements for 6 months, and another profiler measured in the PtR inlet for 8 months. All the instruments were programmed to take a sample every 10 min. In the CdC inlet, tide gauge station (TG) data were collected every 2 min. Measurements also included meteorological variables, including atmospheric pressure, air temperature, relative humidity, precipitation, wind speed and wind direction. The temperature-pressure sensors were deployed near bottom at an average depth of 1.1 m, except for TP-01, which was deployed at 2.5 m for safety reasons, due to vandalism in that region. CP-07 was deployed at 16.7 m and CP-00 was deployed at 7.6 m. The CPs were configured to compute currents in layers of 50 cm. Details regarding the geographical position, measurement periods and sampling rates are given in Table 1. Information about the weather stations is summarised in Table 2. Two of the weather stations were operated by the Mexican Weather Service, and the last station, a Davis 6162 Vantage Pro 2, was installed for this study by the JEST project. The atmospheric pressure was measured to remove the anomaly of the inverse barometer contribution to the water level.

Bathymetric information
The bathymetric chart serial numbers SM 842, SM 842-201, SM 842-301 and SM 84001 were digitalised. To update the information in the two main inlets, bathymetric measurements were collected using an echo sounder. From all the collected data, the total flooded area of the TL was calculated to be 1,936 km 2 and the mean water volume was calculated to be 4.65 km 3 . The area estimated is slightly smaller than the 2,000 km 2 value that is generally used (David and Kjerfve 1998). For the offshore bathymetric information, the GEBCO_08 Grid was used, considering all the data below 200 m. Above that depth, we used information provided by the Mexican Navy. One important bathymetric feature is that the study area is located in an area of the Gulf of Mexico where the continental shelf is approximately 130 km wide. Another important bathymetric feature is the existence of submerged deltas in both inlets. In the case of the PtR inlet, the delta is located inside the lagoon, and in the CdC inlet, the delta is located outside the TL (see Fig. 1).

Tidal analysis (sea surface elevation and currents)
For the harmonic analysis of sea level and current variability, we used the SHOM (Service Hydrographique et Océanographique de la Marine) package (Simon 1974;Bessero 1979). The algorithm computes the harmonic components from the surface variability as well as the harmonic components of the observed current variability. The harmonic components are presented as current ellipses with the magnitude of the major and minor axes expressed in centimetres per second. The major axis orientation is shown in azimuth degrees with respect to magnetic north, as well as the phases, with respect to the Greenwich meridian. The weather stations were equipped with wind speed-direction sensor (WSPD, WD), atmospheric pressure (Pa), relative humidity (RH), precipitation (RAIN) and air temperature (Ta) H is amplitude in cm, g are degrees refered to Greenwich 3 Data

Analysis of the pressure sensor data
All water pressure data were corrected for atmospheric pressure variations and density variations to remove these contributions to the sea level. For the corrected time series, a harmonic component analysis was performed using the SHOM package (Simon 1974;Bessero 1979). The resulting amplitude and phase for the eight more energetic components are shown in Tables 3 and 4. The Dietrich (1944) classification of the tides based on the ratio of the amplitudes (K 1 +O 1 )/(M 2 +S 2 ) provides a measure of the relative importance of the diurnal and semi-diurnal components. The classification is as follows: less than 1.5 mixed, mainly semi-diurnal; 1.5-3.0 mixed, mainly diurnal; and greater than 3, diurnal. The values corresponding to the tide type are shown in the last column of Table 3. The tide range presented in Table 4 was obtained by adding the amplitudes of all the resolved harmonic components obtained with the SHOM package. In this period, the maximal spring tide should occur with the sum of all the components.
The three sensors in the inlets and the two closest to the inlets (CP-08, CP-07, TG-08, TP-01, TP-06) indicated a mixed mainly diurnal tide, and the sensors on the interior of the TL indicated a diurnal tide. The amplitude of K 1 and O 1 are two times larger than M 2 and four times larger than the rest of the considered components. The amplitudes of the diurnal components P 1 and Q 1 are larger than the S 2 amplitude. Inside the lagoon, the tide amplitude of most of the components decreases compared with its values outside the lagoon. The amplitude of K 2 is larger in the northeast region of the lagoon. The phase increases from the inlets into the inner lagoon. The shallow water tidal components inside the TL (Table 5) present a small amplitude in general, although some of them report an amplification inside the lagoon. Due to the small amplitude of those components, we decided to avoid a further analysis of them.
For the remainder of this study, only eight of the resulting diurnal and semi-diurnal components were considered. These components explain more than the 80 % of the observed variability. All components that were not considered had an amplitude smaller than 0.8 cm and were therefore considered insignificant.

Analysis of current data
To analyse the time series of the measured current profiles, the following methodology was applied: For every layer of the current profilers (30 for CdC and 11 for PtR), the current components U and V were individually postprocessed with the SHOM package. The results are shown in Tables 6 and 7.
The currents in the PtR inlet are more intense than in the CdC inlet. In the PtR, the semi-major axis of component K 1 is in the range of 55 to 64 cm s −1 with a mean of 60.45 cm s −1 ; for O 1 , the mean is 55 cm s −1 , and for M 2 , the mean is 51.30 cm s −1 . The velocities of the other components are smaller. In the PtR inlet, K 1 and O 1 are 3.5 times larger than P 1 , Q 1 , K 2 and N 2 .
For the CdC inlet, the largest amplitude corresponds to O 1 , with a mean speed of 46.18 cm s −1 , whereas for K 1 , the value is 35 cm s −1 , and M 2 is of the same order. The rest of the components have smaller amplitudes. For the CdC inlet, the relations between K 1 , O 1 , and M 2 , and the rest of the components are 2.5 times smaller. The maximum variance currents have an orientation of 112°and 125°for the CdC and PtR inlets, respectively (Fig. 2).

Tidal asymmetry analysis
The currents have also been analysed to determine if the inlets are flood-dominant (having shorter-duration and highervelocity floods) or ebb-dominant systems (having shorterand higher-velocity ebbs) (Friedrichs and Aubrey 1988;Dronkers 1986). This tidal distortion phenomenon is one of the main forcing factors in inlet stability and determines if a system is importing or exporting sediments to the ocean (Dronkers 1986).
In the PtR inlet, the flood mean speed is 0.568 m s −1 with a standard deviation of σ=±0.318 m s −1 (southeast) and the ebb mean speed is 0.833 m s −1 with a standard deviation of σ= ±0.51 m s −1 (northwest) (Fig. 3). These effects shift the lagoon filling time to an average of 14.37 h (σ=±4.53 h), and the ebb mean time is 8.50 h (σ=±2.87 h). Therefore, the inlet is ebb dominated.
The analysis of the tidal asymmetry reveals a lagoon system with ebb-dominant tidal distortion. This condition is consistent with a stable erosion-sedimentation condition of  Histogram of the current speeds the inlet systems because the exportation of sediments deposited in the channels (Dronkers 1986) should balance the river contribution. To study the variability of the channel inlets, a more accurate study should be performed that takes into consideration a rising sea level condition.

Numerical modelling of harmonic components
MARS3D from IFREMER was used to study the circulation of the lagoon (Lazure and Dumas 2008;Lazure and Jegou 1998). The governing equations of the hydrodynamic model express the conservation of mass, momentum, and salinity and assume the hydrostatic and Boussinesq approximations (Lazure and Dumas 2008). The equations are solved in the σ-coordinate system using finite differences, and the resolution is based on the mode-splitting technique (Mellor and Blumberg 1985). One characteristic of the model is that the barotropic mode (free surface wave propagation) is semiimplicitly solved using an alternate direction implicit scheme. This allows a coupling with the baroclinic mode (internal motion) using identical time discretisation. Grid cells emerging at low tide have the ability to dry and wet in a massconservative way (Leendertse 1970). Spatial discretisation was performed using a staggered Arakawa 'C' grid (Arakawa and Lamb 1977), resulting in a finite volume type formulation. The model allows nesting using the Agrif2 tools (Lazure and Dumas 2008). This capability was used to refine the grid inside the lagoon.
As a result of the shallow nature of the TL and the significant depth of the inlets, as well as the strength of the friction effects, they are better simulated using a 3D model.

Boundary conditions
Boundary conditions were obtained from a global barotropic inverse tide model (TPX0.7) and are described by Egbert and Erofeeva (2002). Eight constituents (M 2 , K 1 , S 2 , K 2 , Q 1 , O 1 , N 2 , P 1 ) were used to force the sea surface elevation at the boundaries of the larger domain. The amplitude and phase of the TPX0.7 database were interpolated to the resolution of the largest domain. The 3D model was forced by the free surface elevation, and a zero normal gradient condition was used for the velocity along the open boundary. The currents were considered to be equal to zero at the beginning of the simulation. The spin-up time for tidal dynamics is very short due to their high velocity.

Model calibration using sea level data
Using the model results, the amplitude and phase of the main harmonic components were computed using the SHOM package (Tables 8 and 9). The amplitudes calculated are similar to the values extracted from measurements. At sites TP-04, TP-05 and TP-06, the model slightly overestimates the amplitude for the first larger three components. For the S 2 constituent, the model underestimates the amplitude inside the lagoon and overestimates the values at the inlets. For the smaller components, the model and measurements have very good agreement. For K 2 , the model reproduces the observed amplitude increase on the northeast part of the lagoon.
In general, the model phases present good agreement with the measurements, except for S 2 , for which the model overestimates the value with a difference larger than 100°.
To validate the numerical simulation, a comparison between measurements vs. model was conducted, using the nearest model grid point to the corresponding gauging station. The methodology used was as follows: & A correlation was conducted between the time series of water level observations and the time series reconstructed with the harmonic components obtained with the SHOM package for each site. & A correlation was conducted between the model extracted time series and the measured time series. & A correlation was conducted between the synthetic time series obtained with the resulting harmonic components from observations calculated using the SHOM package and the numerical model time series.
The resulting determination coefficients (r 2 ) are shown in Table 10. From the results in Table 10, it is possible to conclude that the time series from the harmonic analysis has a good correlation with the observations. These results indicate that the astronomic tide explains an important percentage of the observed surface variability. The remaining unexplained variability likely represents wind forcing, river discharges and the circulation in the adjacent ocean. The average determination coefficient r 2 between the observations and the model is r 2 =0.78. Although the correlations between the observations and model data are high, they are lower than those between the observations and the series based on the harmonic analysis data. The greater correlation coefficients were obtained between the time series based on the harmonic analysis and the numerical model with an average determination coefficient r 2 =0.93. These results validate the performance of the numerical simulation. The low correlation at point TP-10 can be explained by the length of the measurement period (only 62 days). The analysis requires a longer time period for a more accurate harmonic decomposition.

Calibration using measured currents
As was performed for the sea level variations, the model results were compared to measurements collected using current profilers. In the model grid points closest to the measurement points and along the water column, the currents calculated by the model were extracted. A harmonic analysis was conducted to calculate the ellipses of the current for the different components of the tide (see Tables 11 and 12). The ellipses calculated by the model have characteristics comparable to those obtained from the measurements. For site CP-00 on the surface layer, the semi-major axis calculated by the model for components O 1 , K 1 and M 2 is slightly greater than that from measurements. In the lower part of the water column, the values are closer to the measurements. The main directions of the major axes are very similar, and the differences are approximately 2-3°maximum. The phases at the origin are also comparable. The ellipses of the other components have good agreement between the model and measurement values. For CP-07, the model and measurements displayed similar results. In this case, the model slightly underestimates the value of the major axis of O 1 , M 2 and S 2 . It overestimates the value of K 1 .
To better assess the agreement between the model and measurements, the U and V components of model currents were compared with the measured currents (OBS/MODEL) with data reconstructed from the tidal analysis of measurements (harmonic analysis (HA)/MODEL). The HA dataset is considered as a signal in which the current variability not associated to the tides has been eliminated. Finally, the reconstituted signal was compared with the measured currents (OBS/HA). The results of these comparisons are shown in Tables 13 and 14. The determination coefficient r 2 is slightly smaller near the surface and near the bottom in the case of the current meter at the CdC inlet. The comparison between the reconstructed signal and the model results is particularly good, with r 2 values larger than 0.9 in most cases. The bias for all cases is smaller than 2 cm s −1 , and the RMSE is smaller than 20 cm s −1 for all comparisons. The index of agreement (IOA) (Willmott 1981) was calculated at the points presented in Tables 11 and  12 The comparison between the observations and the model outputs of sea level and velocities indicates good agreement between the two data sets. This result indicates that the model represents the tidal sea level variations and currents well.

Cotidal charts
To study the propagation of the tides in the lagoon, the amplitudes and the phases of eight tidal components were calculated in the model domain, and charts were generated (Fig. 5). Analysing the cotidal charts, it is possible to observe that, in the two inlets, the amplitudes of all tidal components present a similar behaviour. From the ocean to the lagoon, the amplitude decreases very rapidly and the phases delay significantly. To quantify the amplitude and phase gradients from the ocean to the lagoon, the amplitude and phase along a section perpendicular to the inlets were obtained. The locations of these sections are shown in Fig. 6. For K 1 and O 1 , the amplitude decrease is approximately 9 cm in the CdC inlet (see Table 15 and Fig. 6) and slightly less in the PtR inlet, but the phase delay is almost 10°larger in the PtR inlet than in the CdC inlet. For M 2 , the decrease is approximately 9 cm. For the other components with smaller Fig. 5 Cotidal charts of the main eight components, with amplitude in centimetres and phase in degrees relative to the Greenwich meridian amplitude, the gradients are also smaller. Globally, the reduction is approximately 55 % for the diurnal components and approximately 63 % for the semi-diurnal components (see Table 15). Both inlets have large amplitude decreases and a Fig. 5 (continued) phase delay; however, gradients in the PtR are larger. The phases of all eight components grow slowly inside the lagoon from west to east, with a maximum in the eastern side of the lagoon. The phase varies from approximately 50°to 60°. The difference in the gradient behaviour in the inlets could be related to the geomorphologic differences between the inlets.
The PtR channel inlet is shorter and has a delta inside the lagoon. In contrast, the CdC channel inlet is larger and its mean depth is deeper (see Fig. 2).
The amplitude of the different harmonic components outside the lagoon increases between the CdC inlet and the PtR inlet. The phase follows the same patterns. This could be Fig. 6 Change of amplitude and phase delay through the TL inlets for the K 1 component explained by the general propagation of the tide inside the GoM, which propagates from west to east in the southern GoM (Winant 2007). The amplification is generated by the friction effects due to the wide continental shelf in the southeast part of the GoM, as was explained by Clarke and Battisti (1980). The phenomenon is shown in the TPX0.7 database (Egbert and Erofeeva 2002).

Circulation patterns
Current ellipses were calculated in the entire model domain for all layers. Considering the current ellipses' orientation and eccentricity, the different tidal components have similar patterns but with different intensities. Considering the magnitude of the semi-major axis, K 1 , O 1 and M 2 are the more important components. P 1 has current components three times smaller. Q 1 , K 2 and S 2 are four to five times smaller, and N 2 is six to eight times smaller. The surface and near-bottom current ellipses of the component K 1 are shown in Figs. 7 and 8. The main characteristics of the current ellipses are that, in the two inlets, the semi-major axes are larger and, inside the lagoon, they are smaller. The K 1 surface speeds at the CdC inlet reach 51 cm s −1 , which is approximately 5 to 6 cm s −1 in the centre of the lagoon and approximately 1 cm s −1 in the eastern part of the lagoon. For all the components near the bottom, the ellipses are smaller than at the surface, with values of the bottom layer being 10 to 15 % smaller than at the surface. In the PtR inlet, the maximum speed reached is 80.5 cm s −1 .
The area with strong currents is larger in the CdC inlet than in the PtR inlet (Fig. 9a, b). This area extends largely on the two sides of the CdC inlet, reaching over 20 km offshore of the lagoon. In the PtR inlet, the area with high currents is smaller and is mainly located within the lagoon.
In the case of CdC, greater speeds are observed in a 26-kmlong cross section, along the channel. This extension is smaller in the PtR inlet, with a maximum length of 12 km. Another characteristic is that, in the CdC inlet, the faster speeds are outside the lagoon (see Fig. 9), whereas in the PtR inlet, the faster speeds occur inside the lagoon (see Fig. 9). In other areas, the lagoon has weaker currents with a west to east orientation. In the southwest part of the lagoon, in the Atasta lagoons subsystem, the current speed increases because of the narrow and shallow passes between sub-lagoons. These current patterns are very similar at the bottom, with lower speeds, as is possible to observe in Figs. 7 and 8.

Maximum current speeds
A contour plot was generated using the maximum current speed modelled for a 28-day tidal period. The map shows how the area of the current speeds faster than 0.2 m s −1 is more extensive in the CdC inlet than those obtained in the PtR inlet and also regions where speeds are slower than 0.5 cm s −1 . In the middle of the lagoon, the maximum speeds have values around 2 cm s −1 . Based only on tidal currents and their low speeds, the region in the middle of the lagoon would generate a sediment deposition area, but this phenomenon is not observed in the bathymetry. The mentioned feature suggests that, inside the TL, the tidal circulation is important but not the major forcing.
Eight points are shown in Fig. 10. The sea level and current information at these positions was extracted to analyse the variability patterns between them. To obtain the current patterns, the values at opposite points were compared. The points P-01 and P-05 are located in the inlet channels. Therefore, the magnitude of the currents and sea level variation were expected to be of the same order of magnitude. The same procedure was performed between northeast and southwest regions. Figure 11 shows the results of the different comparisons. The U and V velocity components between points P-03 and P-06 show a peculiarity: at point P-03, the negative values of the zonal component are smaller than the positive ones. In the case of point P-06, the negative values of the zonal component Fig. 8 Current variance ellipses of the near-bottom layer. This current field corresponds to the K 1 component Fig. 9 Variance ellipse of surface currents related to the K 1 harmonic component at the Cd. del Carmen inlet and at the Pto. Real inlet are larger than the positive ones. A similar comparison between points P-07 and P-08 reveals that the positive meridional component in P-08 is larger than the negative one, whereas at point P-07, the opposite occurs. This imbalance in the currents causes the existence of a residual flow. The residual flow is discussed further in the following section.

Lagrangian residual velocities
Previous studies have estimated residual tidal transports in the TL. One technique is to follow particles during a tidal cycle. After a tidal cycle, the particles return back near their starting point, but due to friction, the end position is not the same as the initial position. By measuring this difference and dividing by the time, the tidal excursion is obtained as a Lagrangian residual velocity (Jouon et al. 2006;Garreau 1993;Salomon et al. 1988;Andrews and McIntyre 1978). In general, this residual velocity is calculated for each component separately, allowing the following of the trajectories of the particle for a relatively short period proportional to the period of the component.
To obtain better knowledge of the distances travelled by particles, we followed them on a monthly tidal cycle of 28 days (Fig. 12). Inside the lagoon, the long-term transport is relatively slow and usually simple. In regions near the inlets, the trajectories are more complex and the residual larger (Fig. 11). There is always a periodic motion due to the ebb and flood currents, and the trajectories dramatically change the shape between inside and outside of the lagoon, partially due to the geomorphology of the lagoon and partially because the particles can be completely expelled from the lagoon and embedded in areas where they are trapped. In those cases, it is impossible to define a barycentre or a velocity representative of the transport. These results indicate that it is possible to calculate the residual velocities inside the lagoon but not near the inlets. Figure 12 shows the 28-day integrated trajectories in grey. The deployment points are marked in red dots, and the blue lines represent the averaged position every 6-day period. The black stars are the end points. The integrated paths show the residual circulation and the water input in the west inlet (CdC) and the output in the PtR inlet.
To study the effects of the tidal cycle over the residual velocity inside the lagoon, the following methodology was implemented: a 28-day period (a lunar month conformed by two springs and two neap tides) was considered. In each case, to evaluate the flood and ebb flow conditions, Lagrangian particles were released in 148 different positions, starting at the maximal high water level of the spring or neap tide (Fig. 13). These particles are advected in the following 168 h, using the modelled current fields with a time step of 180 s. The same process is repeated Fig. 10 Maximum current speeds after a tidal cycle of 28 days. Speeds are in metres per second every 3 h the following 24 h. In this manner, the differences in the initial conditions during ebb and flood are considered. The speed and direction of every particle is calculated, using the initial and the final position of each Lagrangian trajectory, as  Fig. 10 for location) well as its barycentre trajectory position. Figure 14 shows the results from the first spring tide cycle with its eight residual velocities. The residual velocities calculated were, as expected, very low, reflecting currents of few millimetres per second inside the lagoon. South of the inlets, the velocities are slightly larger than 1 cm s −1 . Inside the lagoon, the general transport is northeastward in the southern part of the lagoon and southwestward in the north of the lagoon. In the central lagoon area, but slightly eastward, the Lagrangian residual circulation displays a large cyclonic gyre.

Discussion
From the analysis of the observations, it is possible to conclude that, outside of the TL, the eight diurnal and semidiurnal components propagate from southwest to northeast. This propagation direction is different than the direction proposed by David and Kjerfve (1998). It is worth to mention that their results are based on a short period of measurements and our analysis is based on more than a year of measurements and includes a better coverage of the lagoon. Under these conditions, we can consider the model results more reliable. As shown in Fig. 6 and Table 15, the amplitudes of the semi-diurnal components are significantly reduced at the inlets, although diurnal components are reduced as well. The difference in amplitude reduction allows the internal TL tidal circulation to be mainly diurnal, which reinforces the already dominant diurnal components in the adjacent ocean.
The reduction of the amplitudes is greater in PtR, which could lead to a smaller diurnal-semi-diurnal interaction, noticeable in a reduction of amplitude of the shallow water tidal components on the northeast region and amplification in the southwest region. The harmonic analysis shows shallow water tidal components are relatively small inside the lagoon; however, some of them amplify on the southwest part. This phenomenon could be related to the fact that the PtR inlet has a semi-submerged sand delta inside the TL and the water flows through the delta channels, increasing the dissipation caused by bottom friction. The tidal asymmetry at PtR makes it an ebb-dominated inlet, but this does not explain the delta inside the lagoon. This issue should be considered in future studies.
Current measurements show that the astronomical tides explain more than 80 % of the variability observed in the inlets' currents. The large variability of the currents associated to the tides and the extension of the areas where the currents are considerably strong, obtained with the model simulation, confirms the importance of a detailed knowledge of the influence of the tidal circulation in the TL.  The numerical modelling shows a complex residual circulation. There is a clear cyclonic vortex inside the lagoon, where the net flow shows water entering through CdC and exiting through PtR. Considering only this circulation, it should be a strong sedimentation area in the middle of the TL because of the low current speeds. This phenomenon is not observed in the bathymetry, which means that the circulation inside the TL is not primarily driven by astronomical tides, but by local wind, river discharges and ocean circulation

Summary and conclusions
Observations of sea surface elevation in different sites of the TL and current profiles in the main inlets were used for tuning a 3D hydrodynamic simulation using MARS3D. The observations and the numerical model were used to understand the influence of the astronomical tide on the circulation in the lagoon and its interaction with the GoM. The principal findings of this investigation are the following: 1. The TL presents a mostly diurnal mixed tidal behaviour outside and inside the lagoon inlets. Because of the inlets' cross-section reduction and shoaling effects, the harmonic components suffer an amplitude decrease and a phase delay at the inlets themselves; however, the dimensions and depth of the lagoon are important in the tidal propagation. For the diurnal components, the averaged amplitude reduction is approximately 45 % relative to the open sea amplitude and there is an average 3.0-h phase delay. The semi-diurnal components suffer an averaged amplitude reduction of 63 % and a 1.8-h phase delay. Inside of the lagoon, the tides are diurnal. 2. Propagation of all tidal components suffers a smooth decay and phase delay with a convergence between Isla del Cayo (north) and Balchacaj lagoon (south). 3. In both inlets, there are ebb-dominated regimes, and this phenomenon is stronger in the PtR inlet. The larger measured speeds are in the PtR inlet during ebb events for spring tides. A study about the import-export characteristics of both inlets should be conducted. In addition, it is necessary to consider the effects of a rising sea level over the coastal lagoon stability in future studies (Zavala-Hidalgo et al. 2011). 4. The hydrodynamic model MARS3D was successfully validated to reproduce the effects of the astronomical tide in the lagoon. As was explained, the numerical model was nested, so the boundary conditions are far away from the  lagoon region; this characteristic allows us to model in a realistic way the areas surrounding the inlets. From the results, is possible to observe how the current fields are influenced, considering that in those regions, in few kilometres, the tide components suffer an important amplitude reduction and phase delay. Based on this fact, considering that in previous studies the boundary conditions were established at the inlets, the described phenomenon is not correctly considered on them. 5. In the inlets, the astronomical tide explains 80 % of the observed variability in both the sea level and currents. Inside the TL, the tide explains more than 70 % of the sea level variability. 6. Tidal modelling reveals a residual current flow from the southwest inlet (CdC) to the northeast inlet (PtR) and a cyclonic circulation inside the lagoon with a vortex core in the northeast.