Residual circulation and freshwater transport in the Dutch Wadden Sea: a numerical modelling study

The Dutch Wadden Sea is a region of intertidal flats located between the chain of Wadden Islands and the Dutch mainland. We present numerical model results on the tidal prisms and residual flows through the tidal inlets and across one of the main watersheds. The model also provides insight into the pathways of fresh water originating from the two sluices at the Afsluitdijk (enclosure dike) through the use of passive tracers. All these results are obtained from threedimensional numerical simulations carried out with the General Estuarine Transport Model (GETM), at a horizontal resolution of 200 m and with terrain-following vertical coordinates with 30 layers. We concentrate on the years 2009–2010, for which we impose meteorological forcing, freshwater discharge, and boundary conditions for tidal forcing and storm surges. Results from the model show an excellent agreement with various observational data sets for sea surface height, temperature, salinity and transport through the Texel Inlet. The simulations show that although tides are responsible for a characteristic pattern of residual transport through the inlets, the wind imposes a large variability on its magnitude and can even invert its direction during strong southwesterly winds. Even though these events are sporadic, they play an important role in the flushing of the Dutch Wadden Sea, as they strongly diminish the flushing time of fresh water. In addition, wind can force a residual transport across the Terschelling watershed of the same order, if not larger, than through any of the main tidal inlets, despite the fact that its tidal prism is much smaller than any of those of the inlets. For the pathways of fresh water, the Terschelling watershed turns out to be more important than was previously thought, while the opposite holds for the Vlie Inlet.


Introduction
The Wadden Sea is located along the coast of the Netherlands, Germany, and Denmark between the European mainland and a chain of barrier islands, and it is connected to the North Sea by a series of tidal inlets found in between the islands.It is the largest system of interconnected intertidal flats in the world and was declared a World Heritage site by UNESCO in 2009.It possesses a large biodiversity, serves as feeding ground for migratory birds (Reise et al., 2010, and references therein), and also hosts economic activities such gas extraction and fishing.In terms of water movement, it presents a complex estuarine dynamic due to the combined effects of the bathymetry, tides, wind surges, and freshwater flows.
The Dutch Wadden Sea extends from the Texel Inlet (the westernmost inlet of the Wadden Sea) at the west up to the Ems-Dollard estuary, which marks the beginning of the German Wadden Sea to the east.Figure 1 shows the configuration of the Dutch Wadden Sea without including the Ems-Dollard Estuary since it does not form part of the present study.There are five barrier islands in the Dutch Wadden Sea with a watershed behind each of them.These watersheds divide the Wadden Sea into a series of tidal basins (Oost et al., 2012).The current paper addresses two basic questions, which have been of interest for several decades, about the hydrodynamics of the Dutch Wadden Sea.What is the residual circulation pattern in the Wadden Sea?What is the fate of the fresh water discharged into it?(In other words, what are the main gateways for exchange with the North Sea?) Two inlets have long been regarded as the principal gateways: the Texel Inlet (label A in Fig. 1) and Vlie Inlet (label C in Fig. 1).They have the largest tidal prisms of the Dutch Wadden Sea, estimated to be 1054 and 1078 × 10 6 m 3 , respectively.The values for the other inlets are much smaller: Eierlandse Gat, 207 (label B in Fig. 1); Borndiep, 478 (label D in Fig. 1); Pinkegat, 100 (label E in Fig. 1); Friesche Inlet, 200 (label E in Fig. 1).These empirical estimates were listed by Louters and Gerritsen (1994) and are all in 10 6 m 3 .
The residual circulation in the western Dutch Wadden Sea (WDWS) has been previously studied, both through numerical simulations and measurements.This region spans from the Texel Inlet to the Terschelling watershed (the watershed between the mainland and the island of Terschelling).Ridderinkhof (1988a) used a vertically averaged numerical shallow-water model with a 500 m horizontal resolution from the Texel Inlet to the Ameland watershed.The forcing consisted of tides alone, i.e. wind forcing and freshwater discharge were not included.He found a residual circulation consisting of an inflow through the Vlie Inlet (1145 m 3 s −1 ), and an outflow through the Texel Inlet (820 m 3 s −1 ), the Eierlandse Gat (250 m 3 s −1 ) and the Terschelling watershed (55 m 3 s −1 ).Using an analytical model, this circulation was explained mainly by differences in water-level amplitude between the Texel and Vlie inlets (Ridderinkhof, 1988b).However, it was later shown, by expanding the analytical model to include wind effects, that the residual circulation between the Texel and Vlie inlets is highly variable due mainly to local wind stresses (Buijsman and Ridderinkhof, 2007b).Between 2002 and2004, Rijkswaterstaat (part of the Dutch Ministry of Infrastructure and Environment) carried out three 13 h ADCP (Acoustic Doppler Current Profiler) campaigns around the Texel Inlet (Elias et al., 2006).They found an average outflow rate of about 2300 m 3 s −1 within a tidal period.In 2012, two 13 h transects were carried out at the Vlie Inlet (Gerkema et al., 2014).They measured net imports through the main channel with an average flow rate of 1300 m 3 s −1 and 3900 m 3 s −1 within the two tidal periods sampled; the difference seems to be associated with wind speed and direction.However, these values are unlikely to be representative of the long-term mean since they are taken during very specific conditions; such measurements are rarely taken during extreme weather events that might be of high importance to the residual flow.The only longterm semi-continuous transport measurements are obtained by the ADCP on board the ferry crossing the Texel Inlet (Buijsman and Ridderinkhof, 2007a, b;Nauw et al., 2014).Buijsman and Ridderinkhof (2007a) estimated with these measurements an average volume outflow rate through the Texel Inlet of 2910 m 3 s −1 for the period 1998-2002.The fate of the fresh water discharged into the WDWS has been a subject of study for even longer than the residual circulation (Postma, 1950;Zimmerman, 1976).This is probably due to the fact that the measurements of the freshwater distribution are easier to perform.On the other hand, no previous numerical results on this topic exist.The fresh water traversing the WDWS originates mainly from two sources (Fig. 1): the sluices at Den Oever and Kornwerderzand in the Afsluitdijk (the dike closing off the inner part of the basin, which was completed in 1932).On a decadal timescale, salinity in the WDWS shows only a slight decreasing trend (van Aken, 2008a), so for the period we are concerned with here (2009)(2010), the mass of fresh water entering from the sluices must be nearly identical to the mass leaving the basin through all inlets combined.Postma (1950) analysed several measurements of salinity at slack tides, including ones during which one sluice was closed.He concluded that all the fresh water from Den Oever leaves via the Texel Inlet, whereas the fresh water from Kornwerderzand is evenly split over the two main inlets, Texel and Vlie.Postma (1950) also emphasized that salinity cannot be used as a tracer indicative of water movements because of the strong mixing in the Wadden Sea.The latter aspect was explored further by Zimmerman (1976).On the basis of new data and a box model, he modified the earlier estimate by Postma (1950) for the fate of the discharge from Kornwerderzand, to 1/3 going towards the Texel Inlet and 2/3 into the Vlie basin (Zimmerman, 1976).Flushing times were estimated to be 15 tidal periods for Den Oever and 31 tidal periods for Kornwerderzand.As for the residual circulation, the residual transport of fresh water is highly variable.This was already pointed out by Zimmerman (1976) as he noticed an influence of the wind in the freshwater distribution and observed a clear decrease in the volume of fresh water inside the WDWS during a large storm surge.
With respect to the fate of the fresh water, two important aspects have still to be addressed.The first aspect is the quantification of the actual export of fresh water through the different tidal inlets.Although it was previously shown that the fresh water goes towards the Texel and Vlie inlets (Zimmerman, 1976), no estimates (from measurements or numerical models) of the net export exist.The model results to be presented here will shed light on the residual transport through the inlets and across the Terschelling watershed.In general, this exchange can be expected to show a lot of variability, depending on the freshwater discharge rate, wind direction and speed, water level surges, and on the phase in the spring-neap cycle.This variability also makes the application of the concept of flushing time less obvious, as will be discussed below.A second aspect concerns the unsuitability of salinity as a tracer, which is compounded by the fact that the fresh water from the two sluices is fundamentally indistinguishable.
In the current paper, we present answers to the questions at hand based on results from three-dimensional simulations of the Dutch Wadden Sea for the years 2009 and 2010 performed using the General Estuarine Transport Model (GETM).The model solves the primitive hydrostatic equations in three dimensions using a realistic bathymetry, freshwater discharge, tides and meteorological forcing.In the simulations, the fresh water is distinguished by tagging it with a different Eulerian passive tracer for each of the main sluices.These realistic simulations allow us to compute the residual water circulation and the residual flow of fresh water in a consistent and inherently mass-conservative way while obtaining a full three-dimensional synoptic view of the Dutch Wadden Sea.In other words, they allow us to revisit the results from the previous work described above, and to quantify the variability of the system by virtue of sufficiently long time series and a good temporal and spatial resolution.
The paper is organized as follows.Section 2 describes the setup: the domain, the bathymetry, and the forcing.For validation, a comparison with different observational data sets is carried out in Sect.3. The tidal prisms at every inlet are presented in Sect. 4. Results on the residual circulation and the fate of the fresh water from the sluices at Den Oever and Kornwerderzand are presented in Sects.5 and 6, respectively.Finally, the conclusions are outlined in Sect.7.

Methods
Fully three-dimensional numerical simulations are carried out using GETM for the years 2009-2010.The numerical model GETM is an open-source code that has been developed specifically for coastal areas.It is now commonly used, and hence, detailed descriptions are abundant in the literature (see e.g.Stips et al., 2004 andStanev et al., 2003).In addition, the documentation and the code are freely available at http://www.getm.eu.For these reasons, a complete description of the model will not be made here, but only the components and specific details relevant for the setup at hand are given in Appendix A.
The simulations include the most realistic possible (in terms of available data and computational resources) values for the bathymetry, atmospheric forcing, boundary conditions, and freshwater discharge.The simulations are initialized from rest on 1 November 2008.For salinity and temperature, climatological values for the month of November are used.Hence, two months are used to spin the model up in both the barotropic and baroclinic sense.This is sufficiently long given the typical flushing time of the basins of about one month (see Sect. 6.2).

Numerical domain and bathymetry
A rectangular domain rotated 17 • in the anticlockwise direction with respect to the east-west axis was defined with its corners located at (4.4693  2a).Within this domain, we defined an equidistant grid with a 200 m resolution using the Rijksdriehoek projection.This projection is the standard projection used by the Dutch Government, and hence it is convenient to use it here.
The Dutch coast is continuously monitored by the Dutch Government -in particular, by Rijkswaterstaat.As part of this monitoring programme, the bathymetry along the whole Dutch coast, from dry land until the 20 m isobath in the North Sea, is measured regularly.Every section of the Dutch Coast, including the Wadden Sea, is sampled at least once every seven years, with certain zones sampled more frequently due to their importance for navigation and other activities.These measurements are processed and made available in a 20 m grid (Wiegman et al., 2005).This data set is known as the vaklodingen, and the whole historical record is freely available (at http://opendap.deltares.nl).To construct the bathymetric map, we considered only the measurement closest in time (either in the past or in the future) to the years to be simulated: 2009-2010.Figure 2b shows the year during which the bathymetry at the different regions was measured.As can be seen, most of the data were measured around the simulated years with some small exceptions.
The process of constructing the bathymetry is as follows.First a 20 m grid with the same orientation as the model grid was constructed.An interpolation of the bathymetry from the vaklodingen was made onto this high-resolution grid.Then, the average depth of the area covered by 10 × 10 points (i.resolution bathymetric map compiled by TNO and the Dutch Navy was used. The bathymetric map was further smoothed for the stability of the model and to avoid numerical noise.This was done in two steps: (1) application of a weighted average filter and (2) using a heuristic linear approach (Dutour Sikirić et al., 2009) (see Appendix B for more details).

Lateral boundary conditions
At the open boundaries of the numerical domain, boundary conditions for the sea surface height (SSH), vertically integrated velocities and vertical profiles of salinity and temperature were applied.On the watershed at the eastern boundary (south of Rottumerplaat), a wall was placed since it is difficult to match the position of the channels in this region with the boundary conditions.
The Flather type boundary conditions are used for the barotropic mode.The vertically integrated velocities and sea surface elevation at the boundaries are taken from results of a two-dimensional model with data assimilation run by Rijkswaterstaat to predict the SSH along the Dutch coast (Plieger, 1999;Verlaan and Heemink, 1998).The extracted data imposed at the boundary have a temporal resolution of 10 min.
Salinity and temperature at the boundaries were obtained from simulations with a set of nested setups.The simulations for these setups were also performed with GETM.First, a two-dimensional North Atlantic setup with 4 nm resolution with only wind forcing was used to model surges.The results of this setup together with the Oregon State University Tidal Prediction Software (OSU-TPS), and climatology for salt and temperature were used for the boundary conditions of a three-dimensional North Sea setup with a 1 nm resolution and 30 vertical layers.In this setup, salinity, temperature, freshwater discharge and complete meteorological forcing are all taken into account.Finally, another threedimensional southern North Sea setup with a 600 m resolution and 42 vertical layers was used.The salinity and temperature data at a temporal resolution of two hours was extracted for the boundary conditions.In addition, temporal relaxation and a four-point sponge layer are used for the salinity and temperature boundary conditions.(More details are given in Appendix D.) We are aware of the fact that the usage of boundary conditions from two different models might cause slight inconsistencies.However, due to data assimilation, the Rijkswaterstaat tide/surge model shows a superior accuracy in the tidal and intratidal variability.
A complete validation and description of the nested setups, and in particular the 600 m southern North Sea, will not be presented in the current paper.The quality of the boundary conditions obtained is evaluated only through the quality of the results within the numerical domain considered in the current paper as shown in Fig. 2a.

Meteorological forcing and surface boundary conditions
Meteorological data were obtained from the operational forecast model of the German Weather Service (DWD).The available variables are wind speed and direction, air temperature, precipitation, cloudiness and dew point, which are discretized in a grid with a resolution of 1/16 • and have a temporal resolution of three hours.The surface stresses τ x s and τ y s are calculated from the wind velocity and are used to set the dynamic boundary condition at the surface: (See also Appendix A for definitions.) The solar radiation as a function of depth I (z) -needed in the transport equation for the potential temperature (Eq.A9) -is calculated according to the expression where I 0 the albedo-corrected radiation normal to the sea surface, a is a weighting parameter and η −1 1 and η −1 2 are the attenuation lengths (Paulson and Simpson, 1977).The parameters a, η 1 and η 2 depend on the amount of suspended material in the water.For the simulations presented in the current paper, a = 0.57, η 1 = 0.37 m −1 and η 2 = 3.54 m −1 .This is equivalent to a Jerlov coastal water type 3 (Jerlov, 1976) according to a non-linear least-squares fit (Stips, 2010).
The boundary conditions at the water surface for potential temperature and salinity are computed using the heat fluxes, precipitation, and evaporation following the work by Kondo (1975).

Freshwater discharge
Fresh water is discharged into the Dutch Wadden Sea from several sluices.Compared to this discharge, the contribution of precipitation minus evaporation is negligible (Zimmerman, 1976).The main sluices in terms of discharge are Den Oever and Kornwerderzand, located in the Afsluitdijk (enclosing dike).This dike separates the Wadden Sea from Lake IJssel (a freshwater lake).To include the freshwater discharge from the sluices in the simulations, the discharge rate for each one of them has to be specified.The location of the sluices implemented in the model are presented in Fig. 3, and the yearly average and maximum freshwater discharge for each of the sluices are presented in Table 1.
For the two main sluices, a record of the total discharge per tidal cycle, the times at which the gates are open and closed, and the water level inside and outside the sluices are available.To implement the discharge in the model, it is important to take into account that discharge only takes place when the gates are open and the water level at the Wadden Sea side is lower than the water level on the inside of the sluice.In other words, considering the average discharge per tidal period would probably not give the correct salinity distribution and variability.Using the available information and assuming that the Bernoulli equation holds for the discharge velocity, we reconstructed the discharge rate with a 10 min temporal resolution (see Appendix C).
For the sluices at Harlingen and Lake Lauwers, a record of discharge rate at a resolution of 15 min is kept by the authorities.In these cases, a simple linear interpolation is done to obtain the discharge rate at a resolution of 10 min.For the sluices at Helsdeur and Oostoever, values for daily average discharge are available.In these two cases, the discharge is spread over the periods when the SSH at the tidal gauge in the port of Den Helder (next to the Helsdeur sluice) is below 0.5 m.This is not the exact condition for discharge to occur, but it reproduces part of the variability.For all other sluices, the daily average discharge is interpolated to a time series with a resolution of 10 min.The modelled discharge rates for these small sluices are not as realistic as for the large sluices, but this is immaterial as they account (all together) for less than 0.5 % of all the fresh water discharged into the Dutch Wadden Sea.

Numerical schemes and computational details
The equations are solved on an Arakawa C-grid.The advection schemes for momentum are TVD-P2-PDM (third order) in the horizontal and TVD-SUPERBEE (second order) in the vertical.For salinity and temperature, the scheme is TVD-SUPERBEE for both horizontal and vertical advection.The same holds for the advection of turbulent kinetic energy and dissipation rate.GETM uses a micro time step for the external mode and a macro time step for the internal mode.The chosen micro time step is 4 s, and the macro time step is 40 s (i.e. a splitting factor of 10).Internal pressure is calculated using the approach by Shchepetkin and McWilliams (2003).In addition, a background horizontal momentum diffusion of 5 m 2 s −1 is used.In general, conservative model parameters (e.g.time step) were chosen to keep the model from crashing during extreme events (storms, large freshwater discharges).However, less conservative values can be used for most of the time.
GETM is an efficient code for parallel runs.The domain is subdivided into 137 subdomains.The simulations are then carried out in a computing cluster with eight nodes of which each has two processors (AMD Opteron 6272 2.6 GHz) with 16 cores each.This gives a total of 256 cores of which 138 are used in the simulations.With this configuration, 11 to 12 h are required to simulate one month.

Comparison with measurements
In order to estimate the accuracy of the simulations, their results are compared with several observational data sets.Although the Dutch Wadden Sea is monitored extensively, this is mainly driven by ecological reasons, storm surge safety and navigation.Hence, few of the data sets available are suitable for a comparison with results from the simulations.In addition, most of the measurements of physical variables are concentrated around the Texel Inlet, while little information www.ocean-sci.net/10/611/2014/Ocean Sci., 10, 611-632, 2014 on salinity and temperature is found elsewhere.In the current paper, the comparison is therefore restricted to: (1) SSH measurements at 14 tidal stations, (2) time series of salinity and temperature at one station with 30 min temporal resolution, and (3) gross transport through the Texel Inlet as measured by the ferry across the Texel Inlet.In all cases, the accuracy of the simulations are quantified using the coefficient of determination R 2 and the root-mean-square (rms) error rms obtained by comparing the measured data and the simulated data.

Sea surface elevation
Within the domain studied, there are 14 tidal gauges recording the instantaneous sea surface elevation every 10 min for the whole simulated period (2009)(2010).To make a comparison, we extract a time series of sea surface elevation at the grid point closest to the position of every tidal gauge.In addition to calculating R 2 and rms , we performed a harmonic analysis of both the simulated and the observed SSH time series using t_tide (Pawlowicz et al., 2002), and we compared the amplitude and phase of the different tidal constituents at all the tidal stations.
Figure 4 shows the position of the tidal stations, and the amplitude and phase lag (in minutes with respect to Greenwich) for the M2 tidal constituent, which is the dominant tidal constituent in the region.Station numbers from 1 to 4 correspond to stations in the North Sea; station numbers from 5 to 11 correspond to stations inside the WDWS; and numbers from 12 to 14 correspond to stations in the eastern Dutch Wadden Sea (EDWS).The tidal wave travels along the Dutch coast from Southwest to Northeast with a (measured) phase lag of 262 min in station 1 to 475 min at station 4. In other words, it takes about 3.5 h to travel from station 1 to station 4. In general, a good agreement is found between the measured and the simulated characteristics of the M2 component.
Figure 5 presents in more detail the error of the simulated characteristics of the M2 constituent.All stations have an error in the amplitude between 2 % to 12 % and an error in phase between −1 min and −20 min (0.2-3 % of the tidal period).A positive error in the amplitude means that the amplitude is underestimated by the simulations, while a negative phase error means that the tidal wave in the simulations is ahead of the measured wave.The largest errors are found for the EDWS stations.This is probably due to such stations being in a location (port/channel) that is not well resolved in the simulation with the current resolution, and to the wall that was placed as a boundary condition across the watershed at the eastern boundary of the domain.However, the error is still within an acceptable range.A similar analysis was performed for the five most dominant tidal constituents (M2, S2, N2, O1, M4) with similar results.For all five components the largest errors are found for the stations in the EDWS.However, in the WDWS, the error in phase remains below 15 min.In terms of error in phase, this is about the best possible result considering that the temporal resolution of the observed SSH time series is 10 min with the time-centring not being specified.The error in the amplitude of the S2, N2, and O1 components in the WDWS remains below 10 %.For the M4 component, the error in amplitude increases in certain stations (stations: 3, 7, 8, 9, 10, 11) to about 20 %.
The comparison discussed so far is based on a harmonic analysis, in which the tidal constituents are extracted from the full signal, which also includes both set-ups and set-downs due to wind.In shallow areas like the Wadden Sea, the propagation of each constituent is affected by these changes in the sea level.This renders the notion of tidal harmonic constants somewhat dubious.In this respect, it makes sense to also compare the full signals.
A correlation analysis between the measured and simulated SSH, and a linear regression was done for all 14 stations.As an example, we present the scatter plot of the observed vs. simulated SSH at station 8 (outside the Kornwerderzand sluice) in Fig. 6.The linear regression (red line) is very close to the ideal result (black line).The slight angle between the two lines means that the simulated SSH range is slightly smaller than the real range.This is consistent with the fact that the amplitude of the M2 component is slightly underestimated by the simulations.The spread around the line is due to the phase difference between the observed and simulated variables.In this case, R 2 = 0.99 and rms = 0.07 m, indicating an excellent agreement between simulated and measured SSH.The rms error remains below 0.12 m for all North Sea and WDWS stations, while it remains below 0.18 m in the EDWS.For all stations, R 2 > 0.96.

Salinity and temperature time series
Time series of temperature and salinity are compared against measurements at the NIOZ jetty, which is located in the Texel Inlet just off the southern coast of the island of Texel measurements are carried out continuously using an Aanderaa conductivity/temperature sensor 3211 about 10 m off the Texel sea dike and at a depth of between 1.5 and 2 m below mean sea level; the actual depth has varied in time within this range.The raw data, with a fixed temporal resolution of 12 s, are smoothed by considering 1 min medians.For the final available product, which is used for the comparison, the 1 min medians are stored only every 30 min.Salinity data are sometimes missing due to instrument fouling even though the instrument was serviced about once a month.For their comparison with the measured data, the numerical results of the closest grid cell were interpolated in the vertical to a fixed depth of 1.75 m below mean sea level.Figure 7a shows the time series of temperature at the NIOZ jetty station in 2009 and 2010 for both measurements and simulations.The seasonal cycle is well reproduced.Figure 7b presents the time series of the daily averaged salinity on the practical salinity scale (PSS).The daily average was computed to observe more clearly the trends over the year.In general, the time series of the simulated data follows well the trends of the measured data with fresher water during late winter/early spring and saltier water towards the end of summer.It can also be seen that the water in the autumn of 2009 is much saltier than the autumn in 2010.Figure 7c shows the time series of salinity at a 30 min resolution for a period of 20 days.It can be seen that the semi-diurnal variability in salinity is well reproduced by the model with some differences in the range of about unity, but without bias.
To quantify the differences between the measured and simulated values for both salinity and temperature, we performed a correlation of the time series.However, for the temperature time series, we compare the temperature anomaly which is obtained by removing the yearly cycle.This is done by fitting a sinusoidal function to the measured temperature time series, and then subtracting it to both the measured and simulated temperature time series.In this way we focus on the variability in weekly, daily and tidal time scales, which are much harder to reproduce than the seasonal variability.
Figure 8 shows the scatter plots of measured vs. simulated values of the temperature anomaly and salinity at the NIOZ jetty for the full years of 2009 and 2010.In the case of temperature, there is a very good agreement with R 2 = 0.85 and rms = 0.68 • C. The difference between the measured and simulated values rarely exceeds 1 • C.This indicates that the tidal and intratidal variability are also well reproduced.The largest deviations were found during the winter at the end of 2010.This might be due to the presence of very low temperature from the end of November, and consequently, sea ice, which is not taken into account in the model.For salinity, the spread is larger (R 2 = 0.78 and rms = 1.22), which is usually the case in numerical simulations of estuaries since salinity depends on the complex pathways of the fresh water, while temperature is forced more locally.However, the linear regression is close to the ideal case.

Transport through the Texel Inlet
Between the port of Den Helder and the island of Texel, across the Marsdiep channel, a ferry passes every 30 min from 06:00 to 22:00 LT.This ferry carries an ADCP and a CTD (conductivity, temperature, depth), which measure current velocities, salinity and temperature.This data set has previously been used to study the dynamics of the Marsdiep channel and was described by Buijsman and Ridderinkhof (2007a) and Nauw et al. (2014).We compare the water volume flow rate through the Texel Inlet as obtained from the simulations and the calculations based on the velocity measurements for the year 2009.A vertical velocity profile with a vertical resolution of 1 m is measured every 1.3 s.The velocity is decomposed into its components perpendicular and parallel to the ferry's trajectory.Then the perpendicular component is integrated over the vertical and along the ferry's trajectory.This approach is more direct than the one previously used by Buijsman and Ridderinkhof (2007a) and Nauw et al. (2014) since it avoids the projection into a regular grid, and hence involves fewer assumptions.However, it still assumes that the data obtained during one crossing are taken at the same time, while in reality, the journey takes about 15 to 20 min.
The comparison between the simulation results and the ferry measurements is presented in Fig. 9.The coefficient of determination is R 2 = 0.98.For the rms error, we obtained rms = 8.00 × 10 3 which is less than 10 % of the maximum volume flow rate.
Due to the irregularity (both temporal and spatial) of the ferry's trajectory, it is difficult to make a detailed one-to-one comparison for individual locations along the transect.This detailed comparison is therefore left for a later paper.The agreement between observations and simulations demonstrates that the model reproduces well the hydrodynamics of the modelled area in all its aspects (sea surface height, current velocities, temperature and salinity).Based on this, we are confident about describing, in the following sections, the exchange through the inlets between Dutch Wadden Sea and the North Sea using the numerical results.

Tidal prisms
As the tides rise and fall in the North Sea, large volumes of water are exchanged through the various inlets.To quantify this exchange, we compute the instantaneous volume flow rate at transects across the inlets and the Terschelling watershed (indicated by red lines in Fig. 10) by integrating the instantaneous velocity over the cross-sectional area.The tidal prism (the water volume going in or going out during a tidal cycle) at every inlet is calculated by integrating over time the absolute value of the flow rate from an instant with zero flow rate to the next.In addition, we define the average tidal prism as the average of the tidal prisms (both during flood and ebb) for every tidal period during the period 2009-2010.
The values for the average tidal prisms at all the transects -as obtained by integrating the velocity across the cross section of the transect -are shown in Fig. 10.The two major inlets in terms of their tidal prism are the Texel Inlet and the Vlie Inlet with about 10 9 m 3 each.In fact, each of them has a larger tidal prism than all the minor inlets combined.As expected, the tidal prism at the Terschelling watershed is much smaller than that at any of the inlets.The values for the inlets confirm the previous estimates listed by Louters and Gerritsen (1994)  first and foremost because they depend on the conditions during which the measurements are made, such as wind direction and spring-neap cycles.For example, Postma (1982) mentions 430 × 10 6 m 3 as the tidal prism for the Borndiep, which (coincidentally or not) is closer to our model value.
Given the accuracy of the model results in reproducing the records of the tidal gauges and the transport measurements by the TESO ferry (as demonstrated above), it is reasonable to regard the model results on tidal prisms as reliable, too.
In fact, the model results allow us to calculate the long-term average tidal prism including a large variety of conditions, something which obviously cannot be achieved by the usual short-term (13 h) measurements.Table 2 presents the average tidal prisms for the inlets and the Terschelling watershed.Here we distinguish between the flood and ebb tidal prisms, which are, in general, different.The difference is indicative of the average residual flow (also indicated in Table 2 and Fig. 11).For convenience, flood at the watershed is defined as having a flow from east to west; i.e. into the WDWS.It can then be seen that there is a net inflow at the Vlie Inlet and outflow at every other inlet.In the WDWS, the inflow at the Vlie is supplemented with an average inflow of fresh water through the sluices of about 20 × 10 6 m 3 per tidal cycle.Remarkably, the residual tidal prism for the watershed (−23 × 10 6 m 3 ) is of the same order of magnitude as those for the Texel and Vlie inlets (18 × 10 6 m 3 and 21 × 10 6 m 3 , respectively) even though the average tidal prism itself is two orders of magnitude smaller.
Table 2 also presents the standard deviation of the tidal prism at every transect which gives a measure of their variability.This variability is partly due to the spring-neap cycle and partly due to wind speed and direction.For all the inlets, the standard deviation is of the order of 20 % of the average tidal prism.Note that the estimates obtained through the measurements mentioned by Louters and Gerritsen (1994) and Postma (1982) fall well within the range of possible values obtained in the numerical simulations.It should also be noted that these empirical values were obtained several decades ago; in the meantime, the prism may have changed somewhat.For the Terschelling watershed, the standard deviation reaches about 200 % of the average tidal prism during flood.The variability is in all cases much larger than the residual tidal prisms.As we will discuss below, the wind (rather than exclusively tides) is a determining factor in the transport across the watershed.In fact, although the tidal prism is defined as the difference between the volumes at high and low water, its actual value at every tidal period is strongly influenced by non-tidal effects, notably wind forcing, which causes it to be highly variable.In conclusion, the average tidal prisms for ebb and flood give a good indication of the long-term net exchange of water volumes through the inlets, but in the short run, the variability overshadows these average residuals.For this reason, we study more in detail the temporal variability of the residual water transport in the next section.As a complementary view to tidal prisms, we consider the volume within the basin.In a basin with several inlets (as the one studied here), flood and ebb at different inlets occur at different times (see phases in Fig. 4), so the tidal prisms at the inlets cannot be simply added up to get the change in volume.Here, we consider the volume of the WDWS basin, which is delimited by the transects across the Texel Inlet, the Eierlandse Gat, the Vlie Inlet, and the Terschelling watershed.Figure 12 shows the time series of that volume for the years 2009-2010 at a 30 min resolution.This volume can be computed from the model output since it contains the water depth at every grid point.The average water volume in this basin for the 2 years of simulation is of 4730 × 10 6 m 3 , while the average maximum volume over the same period is 5724×10 6 m 3 , and the average minimum volume is 3700×10 6 m 3 ; i.e. 43 % of the mean water volume inside the basin leaves and enters every tidal cycle.The difference between the average maximum and minimum volumes is 2024 × 10 6 m 3 , which is smaller (by ∼ 6 %) than the sum of the tidal prisms as calculated from the volume flow rate through the transects limiting this region.This difference reflects the phase shift of the in-and outflow between the inlets.In addition, Fig. 12 shows that the maxima of the volume are highly variable (more so than the minima).The spring-neap cycle can be clearly observed.Besides, there are occasionally very large peaks, with the highest being 7341 × 10 6 m 3 or 1.65 times the average volume, on 4 October 2009.These extreme events can be ascribed to storm surges.

Variability of the residual circulation
For the calculation of residual transport it is, first of all, necessary to have a clear idea over what interval of time the residual is to be taken.At any one location (a tidal gauge, for example), this already poses a problem because the tidal period is not well defined: the time between low waters is, in general, unequal to the time between high waters, and both vary rather erratically over longer periods.This is true even if one were to consider a purely tidal signal (as in tidal predictions), disregarding the effects of wind.This is illustrated for data from tidal prediction at Vlieland Harbour (Gerkema et al., 2014).The underlying cause is that the tidal signal consists of tidal constituents that have incommensurable periods.A further complication is that subsequent low (or high) waters are unequal; in other words, the background state itself changes over the period.In the case at hand, the problem seems to be even more complicated since we have to deal not just with one location but with an entire area, involving differences in phases of the tide of up to 3 h.
In the previous section, where we computed local tidal prisms, the definition of the tidal period was based on the time between alternate slacks.Here, we will consider the multiple-inlet system as an entity, and the previous definition no longer works because of non-simultaneity of slacks.To define a tidal period valid for the whole region, we consider the time series of the water volume inside the WDWS shown in Fig. 12.We define the reference volume as the long-term average volume during the years 2009-2010.During a tidal cycle, the actual volume goes up and down, matching the reference volume at a certain moment during ebb and at a certain moment during flood.The tidal period is defined as the interval between two consecutive moments when the instantaneous volume matches the reference volume during flood.An advantage of this definition is that, for every period, the sum of all the residual inflows and outflows equals zero.(Strictly speaking, the calculation should be in terms of mass rather than volume, since volume is not conserved if a certain mass of fresher water enters and an equal mass of saltier water leaves, but this effect can be neglected here.)As can be seen in Fig. 12, for certain tidal periods the minimum is larger or the maximum is smaller than the average volume.This only means that the calculation of the residual transport will be done over a larger number of tidal periods (usually two or three at most) but still over an integer number of tidal periods according to our definition.
Figures 13a and 14a show the average residual volume flow rate for every tidal period as a function of time for the Texel Inlet, the Vlie Inlet, and the Terschelling watershed for the years 2009 and 2010, respectively.The Eierlandse Gat is omitted since, as already seen in Table 2 and Fig. 11, the residual flow there is much smaller.Most of the time, the residual circulation consists of an inflow through the Vlie Inlet and an outflow through the Texel Inlet, with a relatively small in-or outflow across the watershed.This is consistent with the differences between the average flood and ebb tidal prisms, and also with the results of Ridderinkhof (1988b).This occurs for 63 % of the tidal periods, but the situation can change completely due to wind effects.Figures 13b and  14b show the time series of wind speed and direction at the centre of the Texel Inlet as obtained from the forcing data.It can be seen that all periods with southwesterly to westerly winds in excess of 15 m s −1 result in an inverted residual circulation with an outflow rate through the Vlie Inlet and the watershed of at least 5 × 10 3 m 3 s −1 .These periods are www.ocean-sci.net/10/611/2014/Ocean Sci., 10, 611-632, 2014  Of great interest during the simulated years is a period when strong southwesterly winds blew continuously during about two weeks in November 2009.During this period, a very strong inverted residual circulation took place with a large import of water through the Texel Inlet, and an outflow through the Vlie and across the watershed.Although such events seem to be infrequent (another similar event, although smaller, occurred in December 2011), they could be of importance for the general state and evolution of the Wadden Sea since up to 50 % of the average volume of the basin gets refreshed within a period of about one week.On the other hand, strong easterly winds enhance the residual flow from the Vlie Inlet to the Texel Inlet.This is in agreement with the work by Li (2013) who investigated the effects of wind in an idealized three-inlet system.He found that: (1) wind can significantly influence the residual circulation in a multiple inlet with an inflow usually in the upwind inlet and an outflow in the downwind inlet; (2) the weakest residual flows are found in the middle inlet; (3) wind can produce high-magnitude net flows that tides cannot produce.However, it is clearly important for the WDWS to include the Terschelling watershed, i.e. include at least one more inlet, into such an idealized model.Since wind effects render the residual circulation highly variable, caution has to be exerted when defining the normal or typical conditions.For this, we present the histograms of the residual volume flow rate Q through the three inlets of the WDWS and the Terschelling watershed as shown in Fig. 15.The histograms are constructed by binning the residual flow rate for every tidal period into 100 bins between −8 × 10 3 and 8 × 10 3 m 3 s −1 .In addition, Table 3 presents the value for the median, mean, standard deviation and skewness for each of the inlets and the watershed.Particularly noticeable are the large dissimilarities between the mean and median values.
All four histograms present a clear asymmetry, meaning that the average residual flow over the simulated period is not indicative of typical conditions as it deviates from the peak in the histogram.Instead, the median values correspond well with the peak, giving a better indication of the typical conditions: an inflow through the Vlie of about 600-700 m 3 s −1 and at the sluices of about 450 m 3 s −1 ; and an outflow through the Texel Inlet of about 600-700 m 3 s −1 , and both through the Eierlandse Gat and the Terschelling watershed of about 100-200 m 3 s −1 .The order of magnitude and direction of these residual flow rates is in agreement with the results obtained by Ridderinkhof (1988b), particularly if we consider that the freshwater discharge was not taken into account in his study.Hence, we can assume that the median flow rates are representative of relatively calm winds and are mainly driven by tides.In addition, it can be noted that the typical values of the residual flow rates (including the freshwater discharge) add up to almost null.
As discussed above, much of the variability and the extreme events can be attributed to the varying winds that can enhance, weaken or even invert the tidally driven residual circulation depending on the wind speed and direction.The variability of the residual flow rate (which is represented by the standard deviation) is largest at the Texel Inlet, followed by the Vlie Inlet and the watershed, and finally, by the Eierlandse Gat.The fact that the variability is larger at the Texel and Vlie inlets can be attributed to the existence of a relatively good connection between the two inlets.On the other hand, the low variability at the Eierlandse Gat suggests that this inlet is rather isolated from the rest of the WDWS by the Texel and Vlieland watersheds.The large variability of the residual flow rate in the Terschelling watershed is also clearly associated with the wind.
The histograms of the residual flow for both the Texel Inlet and the Eierlandse Gat are positively skewed, while the histograms for the Vlie Inlet and the watershed are negatively skewed.This is a reflection of the predominant wind direction both in average and during most of the major storms, which results in more tidal periods where the typical tidally driven residual circulation with an inflow through the Vlie Inlet is weakened or reversed.Furthermore, note that the histogram for the residual flow at the watershed has a larger skewness and peakedness than those for the Texel Inlet and the Vlie Inlet.This shows that the transport through the watershed has a clear preferential direction, and explains the large net transport.However, even in the preferential direction, strong winds seem to be needed to drive a transport across this watershed.

Spatial distribution of fresh water
In order to study the freshwater distribution and transport, the fresh water from the two main sluices (Den Oever and Kornwerderzand) is tagged with an Eulerian passive tracer following a similar approach to Zhang et al. (2009) and Meier (2007).A different tracer is associated with each of the sluices, and the fresh water has a concentration equal to unity (c i = 1, where i refers to the sluice).In this way, we can distinguish the fresh water from the main sluices from the fresh water from other sources (the Rhine, other sluices, and precipitation) in an unambiguous way.
For the advection of the passive tracers, the Framework for Aquatic Biological Modeling (FABM) passive tracer module (see http://fabm.sourceforge.net/)coupled to GETM is used.GETM solves for the advection of the passive tracers using the same methods as for salinity and temperature.The initial concentration in the domain is set to zero, and the model is spun up during two months.We consider this time to be sufficient since the typical transport timescales in the Wadden  Sea are shorter than a month.The flushing time of the fresh water will be discussed in more detail in Sect.6.2.
In describing the fate of fresh water, we could mean two different things: how the water gets saltier along the way, or where the fresh water goes in the sense of Lagrangian tracers (e.g.drifters, buoys).The latter is purely advective, the former is diffusive as well.A tracer concentration, defined with a certain value at the sluice, but initially zero elsewhere, behaves identically to salinity if the diffusion coefficients are the same and if that sluice acts as the only source as in the work by Zhang et al. (2009).When dealing with multiple freshwater sources as in the Wadden Sea or Baltic Sea (Meier, 2007), a couple of points must be clarified.
1. Since there are multiple sources, the way the fresh water from one sluice gets saltier is influenced by the other sources.For example, if fresh water from Den Oever meets fresh water from Kornwerderzand there is no diffusion since the salinity is identical.If on the other hand, the Den Oever fresh water would instead meet North Sea water, there is a gradient in salinity.
2. Since the advection-diffusion equation is linear in salinity or tracer concentration, the working of the equation is the same for these quantities, even if they differ by a constant factor or by an added constant.
From the above two points, it follows that the behaviour of the tracer released in Den Oever would be identical to that of fresh water if the Den Oever sluices were the only source of fresh water.We can follow the tracer concentration, calculate the amount in a volume V : and calculate the transport through the transects delimiting the WDWS: with j denoting the transect of interest, u n the velocity component normal to the transect, and A j the cross-sectional area of that transect.
Figure 16a shows the average distribution of the concentration in the uppermost layer of the Den Oever tracer and Fig. 16b of the Kornwerderzand tracer for April 2009.In general, the Den Oever tracer remains in the southern part of the WDWS close to the Texel Inlet, while the Kornwerderzand tracer is distributed all the way from the Texel Inlet to the Friesche Inlet.It is clear that for this distribution to occur, the fresh water from Kornwerderzand has to be transported through the Terschelling watershed.Figure 16c shows the average salinity distribution in the uppermost layer.By comparing the salinity distribution to the tracer distribution, it can be observed that the regions of low salinity are a superposition of the regions of high tracer concentration.
High concentrations of the tracers on the tidal flats or other shallow areas are not equivalent to large volumes of fresh water.In other words, a map of concentration in such a complex bathymetry as the one of the Wadden Sea, gives little information on the position and destination of the fresh ter. Figure 17 shows the average distribution of the volumes of fresh water for the month of April 2009.While the concentrations are higher on certain tidal flats, the largest volumes of fresh water are located in the channels.The fresh water from Den Oever is located in the channels leading to the Texel Inlet, and the fresh water from Kornwerderzand is distributed along the channels in the whole of the WDWS and even around the Borndiep.This is in line with the findings of Postma (1950) and Zimmerman (1976), who found that the fresh water from Kornwerderzand gets distributed in almost equal parts between the Texel and Vlie basins, whereas the water from Den Oever travels mostly towards the Texel Inlet.We now examine where it actually leaves the WDWS.
Figure 18a shows the time series with a 30 min resolution of the instantaneous flow rate of the Kornwerderzand tracer through the Texel Inlet, the Vlie Inlet and across the Terschelling watershed.As can be seen, the magnitude of the instantaneous transport is higher at the inlets than at the watershed.The tracer changes its destination inlet as a function of time depending on the conditions (e.g.wind direction and freshwater discharge rate).Once the tracer reaches one of the inlets, it does not exit the Wadden Sea permanently.It goes in and out through the inlets with the tide, and a clear outgoing residual trend cannot be deduced from the instantaneous flow rate.To compute the residual flow of the tracers through the inlets, we use the same approach to the one used in the previous section to calculate the residual water transport.Figure 18b and c show the residual transport through the Texel Inlet, the Vlie Inlet and across the Terschelling watershed of the tracers associated with the two major sluices.The tracer from Den Oever (Fig. 18b) exits mainly through the Texel Inlet.A few exceptions occur during strong southwesterly winds when there is an outflow through the Vlie Inlet and across the watershed.These periods correspond with the periods when the residual circulation in the WDWS is inverted such as was the case in November 2009 (see Fig. 13).During the period 2009-2010, the Den Oever tracer exited the WDWS 81 % through the Texel Inlet, 12 % through the watershed, 3 % through the Vlie Inlet and 4 % through the Eierlandse Gat.
For the water out of Kornwerderzand (Fig. 18c), there is a more balanced distribution of the residual outflow through the possible outlets.Remarkable, however, is the importance of the residual outflow across the Terschelling watershed, which is usually disregarded.In fact, 40 % of the Kornwerderzand tracer exits across the watershed freshening the water in the contiguous basin, which corresponds to the Borndiep Inlet.About an equal amount of the tracer from Kornwerderzand (43 %) exits through the Texel Inlet.Surprisingly, the Vlie is not an important outlet for this fresh water with only 10 % of the total outflow, which is comparable to the 7 % of the outflow through the Eierlandse Gat.The time series of the outflow of the tracer through the Terschelling watershed exhibits large peaks.A visual comparison between the residual outflow of the tracers and the wind  shows that the large peaks are related to strong wind events.
For both tracers, the total averaged outflow across the transects was compared with the inflow at the sluices.The computed outflow is within 5 % of the inflow, giving confidence in the results.Note that in principle there is no reason for them to be identical since there is a lag between the inflow and the outflow.

Flushing frequency of fresh water
Transport timescales, such as the residence time and the flushing time, are used to estimate the renewal and ventilation of estuaries and other coastal areas, and are useful to explain several processes of importance in biology and ecology (see e.g.Abdelrhman, 2005;Monsen et al., 2002, and references therein).Note that actual names and definitions vary greatly within the literature (Abdelrhman, 2005).The flushing time T f is commonly defined as the ratio of the mass of a constituent in a water body M i to its rate of renewal (inflow) F i , i.e.
where the subscript i is associated with origin of the specific constituent; see e.g. the work by Bolin and Rodhe (1973) where the flushing time is actually referred to as turn-over time -and by Zimmerman (1976).However, Eq. ( 6) was derived for a steady-state condition, which is not the case for the Wadden Sea for the timescales relevant here.In addition, we are interested in studying the variability of the flushing time of fresh water and its dependence on the external forcing.Some of the limitations of considering Eq. ( 6) for highly variable systems have been previously discussed by Monsen et al. (2002).
For the steady-state case, the inflow and the outflow are identical, and so, the flushing time (Eq.6) can be defined as the ratio of the mass of the constituent in the water body to the outflow.For the unsteady case, this definition means that the flushing time is the time it would take to empty the available mass at the current outward volume flow rate.For the case studied here, this still presents a difficulty since the outflow (and the inflow) are null during some periods, resulting in infinite flushing times.We find it useful to define the tidally averaged flushing frequency (the inverse of the flushing time) for every tidal period T as where a sum of the residual volume flow rate Q i over all the inlets and the watershed (sum over j ) is considered, and the volume of the specific constituent is used instead of the mass.
Figure 19 shows the tidally averaged flushing frequency i as a function of time for the tracers associated with Den Oever (DO) and Kornwerderzand (KDZ).Negative frequencies mean that there is a net import of fresh water from the North Sea or the adjacent basin to the East.This can occur, for example, after a large flushing event when a large portion of the fresh water is located outside the WDWS, especially if the wind rotates close to 180 • as during the passage of cold fronts.
A typical value for the flushing frequency of the fresh water is given by the median frequencies: DO = 0.0363 days −1 (equivalent to a flushing time 1/ DO = 27.5 days) and KDZ = 0.0259 days −1 (equivalent to a flushing time 1/ KDZ = 38.6 days).These estimates for the flushing times are 2.5-3 times longer than the ones found by Zimmerman (1976) of about 8 days for the fresh water from Den Oever and 16 days for the fresh water from Kornwerderzand.Nonetheless, the values obtained by Zimmerman (1976) are still well within the possible range.The difference can be due to the conditions during the period that the measurements were taken since the flushing time is highly variable.In addition, in Zimmerman's box model, there is a sensitivity of the calculation of the freshwater volume to the choice of a reference salinity, which may explain discrepancies with our results.
There are important flushing episodes which are again associated with strong winds (Fig. 19).The arrows in Fig. 19 correspond to the same strong wind events signalled in Figs.13-14.During the years 2009-2010, the most important flushing episode occurred in November 2009 when the flushing frequency reached a value larger than 1 day −1 .This event was characterized by strong persistent southwesterly winds (also discussed in Sect.5) that forced a strong persistent outward flow that drastically reduced the volume of fresh water inside the WDWS.In fact, the volume of fresh water inside the basin decreased one order of magnitude in 10 days from 550 × 10 6 m 3 on 18 November 2009 to 65 × 10 6 m 3 on 27 November 2009.The end of this flushing period corresponds to some of the highest salinities registered at the NIOZ jetty, as shown in Fig. 7b.

Conclusions
In this paper, we have shown model results for the hydrodynamics of the Dutch Wadden Sea, with a focus on the western Dutch Wadden Sea, for the years 2009-2010.For this purpose, we used the GETM/GOTM model at a 200 m resolution in the horizontal and with 30 vertical layers.Boundary conditions impose tides, wind surges, freshwater discharges, wind stress, precipitation, evaporation and surface heat fluxes.
Comparisons with different kinds of observations show an excellent agreement.These include a comparison with the observed SSH at 14 tidal gauges, with transport flow rates through the Texel Inlet as derived from ADCP measurements beneath the TESO ferry, and with measurements of salinity and temperature at the NIOZ jetty.The agreement between these data sets and the simulations makes it reasonable to assume that the model results are close to reality also at places where no measurements are available to compare with.Consequently, we have created as a by-product a solid 2-year hindcast database of the hydrodynamic conditions of the Dutch Wadden Sea for further use by sediment transport, and especially, ecological studies that can benefit from an upgrade to an overall spatial resolution of 200 m and a temporal resolution of 30 min.
Qualitatively, conclusions from earlier studies are mostly confirmed by our results.For example, we also find that most of the fresh water discharged at the sluice at Den Oever leaves via Texel Inlet, in agreement with Postma (1950) and Zimmerman (1976).We find, on average, a residual flow from the Vlie Inlet to the Texel Inlet, confirming Ridderinkhof (1988b), and average values of the tidal prisms of the various inlets are close to earlier (empirical-based) estimates listed by Postma (1982) and Louters and Gerritsen (1994).
In addition, the present study sheds light on the key aspect of variability due to wind, which had previously been unapproachable for a quantitative analysis, either because measurements lack spatial coverage and are mostly biassed to mild weather, or because models lacked wind forcing.In fact, this is the first three-dimensional modelling study of the Dutch Wadden Sea that includes, besides tides, realistic wind forcing and freshwater discharge.
This has allowed us to calculate not only the tidal prisms for each inlet but also their variability (standard deviation), which turns out to be typically about 20 % of the average tidal prism.For the Terschelling watershed, the variability even exceeds the average tidal prism.In all cases, the variability in the tidal prism is much larger than the residual prism, the long-term average of the difference between ebb and flood prisms (Table 2).
We also examined the time series that show the episodes responsible for this variability.We find that, for each inlet, the probability distribution of the residual flow rate spreads around the typical value which is approximated better by the median than by the average.In fact, the median and mean are very dissimilar (Table 3) due to a clear skewness of the distributions dictated by the wind's direction.In fact, the threeinlet system of the WDWS is reminiscent of the idealized system studied by Li (2013) in which the outer inlets presented the largest response to wind and the middle inlet had the weakest.However, the large residual flow across the Terschelling watershed during strong southwesterly winds is a crucial component, whose role in the overall circulation in the Dutch Wadden Sea is much larger than was previously assumed.
To calculate the residual transport, we propose an alternative to the more common (tidal) harmonic analysis on the time series of transport flow rate or velocity.There, the residual of the harmonic analysis is assumed to be the residual flow.However, this method fails here because set-up by winds affects the tidal propagation, and hence, the total signal is not just the linear sum of independent components.This is particularly true at the Terschelling watershed, since due to its shallowness, it can be dry for several hours a day.This clearly renders the time series of the transport rate highly non-linear.For this reason, we adopted an alternative way by calculating the residuals directly from time integration of the volume flow rates.The definition of the time interval, i.e. the tidal period, is based on the volume inside the basin, and hence, it is valid for the whole basin.
With regard to the discharge of fresh water from the sluice at Kornwerderzand, it had never been established previously where it actually leaves the Wadden Sea, although it was presumed that the Vlie Inlet was a major exit route along with the Texel Inlet (Postma, 1950).It now becomes clear that the Vlie Inlet, in fact, plays only a minor role in the net exchange and that, instead, the Terschelling watershed together with the Texel Inlet are predominant.
Most of the time, the hydrodynamic state of the Dutch Wadden Sea, with regard to the residual circulation and freshwater flow patterns, is close to the typical, i.e. median state.However, at other times it is upset or reversed with a skewed tendency due to the wind direction.As a result, the yearly mean values of residual volume flow rates or flushing times are not necessarily representative.The variability here has an episodic character, and this is a key feature of the net transports and flushing times.This is the main outcome of this paper, along with the quantitative substantiation of it that the model provides.

Figure 1 .
Figure 1.Map of the Dutch Wadden Sea up to the beginning of the Ems-Dollard estuary.The names of the islands and inlets are indicated.The position and name of the two major sluices are also indicated; more details on these and other sluices are given in Sect.2.4.The grey line indicates the −5 m isobath.
Figure 2. (a) Numerical domain and bathymetry.(b) Dates of the measurements used to compose the bathymetry up to the 20 m isobath in the North Sea.

Figure 3 .Figure 4 .
Figure 3. Location of sluices discharging fresh water into the Dutch Wadden Sea.The sluice numbers correspond to those inTable 1.

Figure 5 .
Figure 5. Amplitude error (in percentage of the measured amplitude) and phase error (in minutes) of the simulated M2 tidal constituent when compared with the measurements.The dashed lines divide North Sea, WDWS and EDWS stations.

Figure 6 .
Figure 6.Scatter plot of observed vs. simulated SSH (in metres) at station 8 (Kornwerderzand).The red line represents the linear regression given by a coef + b coef × SSH observation .The centre black line represents the ideal linear regression (simulation = observation), while the side dashed black lines are the ideal regression ±0.3 m.The coefficient of determination R 2 and the rms error rms are also indicated.

Figure 7 .
Figure 7.Comparison of temperature and salinity measurements at the NIOZ jetty station against results from the numerical simulations.(a) Time series of temperature with a 30 min resolution for the years 2009-2010.(b) Time series daily averaged salinity for the years 2009-2010.(c) Time series of salinity with a 30 min resolution for the days from 25 March 2009 to 15 April 2009.Gaps in the salinity time series are due to malfunctioning of the sensor or biofouling.The numerical results were interpolated to a fixed depth of 1.75 m below mean sea level.

Figure 8 .
Figure 8. Scatter plots of measured vs. simulated (a) temperature (every 30 min), (b) salinity (every 30 min).The red lines correspond to a linear regression.The centre black line represents the ideal linear regression (simulation = measurement), while the side black lines are the ideal regression ±1.For every plot the coefficient of determination (R 2 ) and the rms error min are indicated.

Figure 9 .
Figure 9.Time series of volume flow rate across the Texel Inlet as obtained from the TESO ferry measurements (black dots) and from the numerical simulations (red line) for 16 days in August 2009.

Figure 10 .
Figure 10.Average tidal prisms (10 6 m 3 ) of each inlet and the Terschelling watershed.The red lines mark the position of the transects used to compute the tidal prisms and transport.The numbers represent the tidal prism in 10 6 m 3 .

Figure 11 .
Figure11.Average residual tidal prisms (flood − ebb) (10 6 m 3 ) of each inlet and the Terschelling watershed.The average input of fresh water per tidal period through the Den Oever and Kornwerderzand sluices is also presented.The numbers (in 10 6 m 3 ) and the length of blue arrows represent the magnitude of the residual tidal prism.The direction of the arrow represents the direction of the average residual tidal prism.The red lines mark the position of the transects used to compute the tidal prisms and transport.

Figure 12 .
Figure 12.Time series of the volume inside the WDWS (black line).The grey line represents the mean volume over the period 2009-2010.The dashed grey lines represent the average maximum and minimum volumes over the period 2009-2010.

Figure 13 .
Figure 13.(a) Residual volume flow rate as a function of time for the Texel Inlet, the Vlie Inlet, and the watershed in 2009.Positive flow rate corresponds to an inflow, while negative flow rate corresponds to an outflow.(b) Wind speed and direction next to the Texel Inlet as a function of time.The wind velocity is obtained from the data used for the atmospheric forcing.The black arrows point to strong southwesterly to westerly wind events which correspond to periods when the residual circulation is inverted with respect to the typical case.

Figure 15 .
Figure 15.Histograms of the residual flow rate through the different inlets of the WDWS and the Terschelling watershed.The solid red line represents the median value, and the dashed lines represent the median value plus/minus one standard deviation.The green line represents mean value.

Figure 16 .
Figure 16.Average concentration in the uppermost layer of the passive tracer from (a) Den Oever and (b) Kornwerderzand for the month of April 2009.(c) Average salinity in the upper layer for the same month.

Figure 17 .
Figure 17.Average volume (10 3 m 3 ) per horizontal grid cell of the fresh water from (a) Den Oever and (b) Kornwerderzand for the month of April 2009.
Figure 18.(a) Time series of the instantaneous volume flow rate of the tracer corresponding to water from Kornwerderzand at the Texel Inlet, the Vlie Inlet, and the Terschelling watershed.The temporal resolution is 30 min.(b) Tidally averaged volume flow rate of fresh water from Den Oever through the Texel Inlet, the Vlie Inlet and across the Terschelling watershed.(c) Tidally averaged volume flow rate of fresh water from Kornwerderzand through the Texel Inlet, the Vlie Inlet and across the Terschelling watershed.

Figure 19 .
Figure 19.Flushing frequency (days −1 ) of the fresh water from the sluices of Den Oever and Kornwerderzand. i

Table 1 .
Average and maximum freshwater discharge in m 3 s −1 at the different sluices for the years 2009-2010.The sluice numbers correspond to those in Fig.3.

Table 2 .
Average, standard deviation (SD), and residual of the flood and ebb tidal prisms for all inlets and the Terschelling watershed.For convenience, flood at the watershed is defined as having a flow from east to west.

Table 3 .
Median, standard deviation and skewness of the residual flow rate at the inlets of the WDWS and the Terschelling watershed.