Calculating residual flows through a multiple-inlet system: the conundrum of the tidal period

The concept of residual, i.e., tidally-averaged, flows through a multiple inlet system is reappraised. The evaluation of the residual through-flow depends on the time interval over which is integrated, in other words, on how one defines the tidal period. It is demonstrated that this definition is ambiguous and that different definitions (based on, e.g., high waters, slack tides, etc.) yield very different results for the residual, also in terms of their long-term statistical properties (median and standard deviation). A basin-wide applicable method of defining the tidal period, in terms of enclosed water volume, is analyzed. We compare the different methods on the basis of high-resolution model results for the Western Dutch Wadden Sea. The multitude of tidal constituents together with wind variability creates broad distributions for the residuals, with standard deviations much larger than the mean or median residual flows.


Introduction
Multi-inlet estuarine and lagoon systems are a common feature along the coast and possess high economical and ecological importance. Examples of this type of systems are the Venice Lagoon (Italy), The Wadden Sea (The Netherlands, Germany, and Denmark), and several found along the coast of the USA such as the Indian River lagoon, Florida; Palmico Sound, North Carolina; and the whole of Georgia's coast. The water exchange between these coastal systems and the adjacent shelf seas is the main way to refresh the system, expel pollutants, and import nutrients and sediment.
This exchange is largely due to the tidal cycles of ebb and flood. However, if the flow is integrated in time over a tidal period, one finds that there is also a residual flow. This residual flow can be driven by the interaction of tides and the topography (Zimmerman 1978), density gradients (see, e.g., the review work by MacCready and Geyer (2010), Geyer and MacCready (2014)), wind (Smith 1990;Li 2013), and the inflow of fresh water into the system. In the current paper, we are not interested in these processes but in how to compute the magnitude of the resulting net residual flow through a multiple-inlet system. In such a system, one may find that, on average, flood dominates in some inlets, while ebb dominates in others. This may be due for example to differences in tidal phase, amplitude, mean sea level, and inlet dimensions (van de Kreeke and Cotter 1974;van de Kreeke and Chiu 1980). Other factors such as the fresh water discharge and wind modify or can even dominate over the tidally-induced through-flow (see, e.g., Smith 1990; Buijsman and Ridderinkhof 2007). The total residual through-flow is then obtained by integrating the residual flow velocities over the cross-section of the inlets connecting the basin with the sea, or in an equivalent way, by integrating the instantaneous volume flow rate through the inlets over a tidal cycle. This quantity can be used to characterize the system by computing, for example, typical transport time scales (see, e.g., Bolin and Rodhe 1973;Monsen et al. 2002). In addition, the residual through-flow can play an important role in the stability of the inlets (Salles et al. 2005).
The concept of the residual through-flow seems straightforward, but it is hard to measure. Its magnitude is usually much smaller that the maximum instantaneous flow rate, and moreover, one needs to cover the entire cross-sections of all the inlets over a year or longer (to take into account the variability due to wind). In numerical modeling, the necessary spatial and temporal coverage is more easily attained, an early example being the study by (Ridderinkhof 1988) of the Western Dutch Wadden Sea (WDWS). He included the dominant tidal constituent M 2 (and its internally generated higher harmonics), but no wind effects. Likewise, the theoretical model by Loder (1980) for the residual circulation around George's Bank concerns the dominant tidal constituent alone. The residual flow then refers to the flow integrated over one tidal period defined by that single dominant constituent, i.e., every tidal period has a constant duration of 12 h 25 m 14 s, the period of the M 2 constituent. As a result, the residual flow can be determined without ambiguity.
In reality, however, it is far from obvious what one should take as the duration of the tidal period. (The underlying causes of this ambiguity will be discussed in Section 2.) This problem seems to have received little attention in the literature, or perhaps it has not been perceived as a problem at all. One reason may be that the focus in tidal analysis is often on the (main) individual tidal constituents, for which the periods are well-defined. Indeed, the harmonic method developed by Kelvin exploits this fact (see, e.g., Pugh 1987;Cartwright 1998) making it possible to predict high and low waters precisely by adding up the different constituents after their amplitudes and phases have been determined empirically for the location in question. The period between subsequent high (or low) waters is then simply an outcome of this method. Another reason for neglecting this problem may be that the main interest was in computing a representative quantity such as the yearly average residual flow through the inlets. For such quantities, the definition of the tidal period is not as relevant since one integrates over a much longer period. However, the tidal period is relevant when the focus lies on variability at shorter time scales, e.g., due to wind.
Hence, the central point of the present paper is to examine how the mean, median, and the standard deviation of tidally averaged residual flow depend on how one defines the tidal period, i.e., on the limits of integration. We discuss the advantages and disadvantages of several possible definitions (the time between alternate mean sea-level crossings, the time between alternate slacks, a constant value (M 2 ), and the time between alternate mean-volume crossings).
As a case study, we consider the WDWS (delimited by the red transects Fig. 1): this is the western part of the Wadden Sea, a coastal lagoon which extends along the coast of The Netherlands, Germany, and Denmark. The WDWS is connected to the North Sea by three inlets: the Texel Inlet (also known as the Marsdiep), the Eierlandse Gat, and the Vlie Inlet. To the east, it is connected to the eastern Dutch Wadden Sea through the watershed (i.e., tidal divide) between the island of Terschelling and the mainland. We study the problem using a realistic numerical model for the Dutch Wadden Sea, covering the years 2009-2011; a brief description of the model is given in Section 3. Results from this model were already used to study the residual flow of water and fresh water through the WDWS by Duran-Matute et al. (2014). It was observed that the residual flow is highly variable, mainly due to wind. Statistical measures (mean, median, standard deviation, and skewness) were used to characterize the magnitude of the through-flow. To define the tidal period, they selected phases of rising tides and picked the moments when the water volume matches the reference volume, the long-term mean. The time between consecutive matches was then taken to be the tidal period. The instantaneous volume flow rate was then averaged over each tidal period for the years 2009-2010. The questions that we tackle now is whether, how, and why the choice of the definition of the tidal period affects the resulting statistical measures for the tidally averaged through-flow. Furthermore, we examine the sensitivity of the volume method to the choice of the boundary of the tidal basin(s) and explore the applicability of this method for observational data sets (which are spatially less comprehensive than model data).
In this paper, we focus on the residual transport of water itself, but the relevance of the problem at hand-how to define the tidal period?-extends directly to sea-level residuals as well as residual transports of sediment, nutrients, pollutants, etc. in multiple-inlet systems.

The tidal period-what is it?
In back-barrier systems like the Wadden Sea, large amounts of water are flushed back and forth through each inlet with ebb and flood. What interests us here is the net result; this means that we have to integrate the (instantaneous) transport over a tidal cycle. However, the duration of that cycle appears to be ambiguous. Three distinct causes underlie this problem.
The first cause derives from the nature of the tidal constituents. Each of them has a well-defined mean period, which originates from the astronomical motions. In general, however, the periods of different constituents are incommensurable, their fractions being irrational. This fact, too, has its origin in astronomy, with its incommensurable periods of the mean solar day, synodic month, tropical year, lunar nodal cycle, etc. This is also the reason why it has proved to be impossible to construct a properly working calendar based on those periods (Richards 1998). For the tides, the implication is that the sum of the tidal constituents shows a distinct lack of periodicity. Thus, the tidal period becomes an elusive concept. For example, the time between high waters will in general be different from the time between low waters, and both vary in an erratic way. This is even true if one considers purely tidal prediction signals that involve no wind effects (an example is shown by Gerkema et al. (2014)).
Second, the effects of wind blur the definition even more. For example, the time between consecutive low waters can be prolonged if the second low water is delayed by a surge due to a storm.
As an example, we show in Fig. 2 the duration between consecutive low-waters and between consecutive highwaters at a tidal gauge at Vlieland Harbour inside the WDWS at the eastern tip of Vlieland. Clearly, the tidal periods based on consecutive low and consecutive high waters can differ considerably, sometimes more than an hour. For both definitions, the spread is large, with values lying anywhere between 11 and 14 hours. In the course of time, the periods vary in a complex way, although one can discern a diurnal inequality as well a variation with the spring-neap cycle (periods are generally shorter during spring tides).
Collecting such results for 3 years of data (2009-2011), we can make histograms of the duration of the tidal period based on consecutive low or consecutive high waters. They are shown in Fig. 3a,b. We have added a third possible definition of the tidal period: we consider consecutive transitions through the (local) mean sea level during phases of rising tides. (Sometimes no such transitions occur when the set-up due to wind keeps the water level elevated above the mean-sea level for some days; in such cases, the values for   Fig. 3, the conclusion must be that the different ways of defining the tidal period can result in quite different histograms of those periods. As an aside, we draw attention to an intriguing coincidence hidden in the tidal-gauge records. If we calculate the long-term average for each of the definitions used in Fig. 3, involving 3 years of data, we obtain 12 h 25 min 13 s (average of tidal periods based on consecutive low waters), 12 h 25 min 15 s (based on consecutive high waters), and 12 h 25 min 13 s (based on alternate transitions through mean sea level). This is the M 2 period to within 1 s! At first sight, it is puzzling why this specific period should come out of a complex time series that involves many other tidal constituents as well as wind effects. However, the underlying cause is purely mathematical. This can be seen from a simple exercise: if one adds two sinusoidal waves with unequal amplitudes and with periods whose fraction is irrational (e.g., those from the M 2 and S 2 tidal constituents), then the period between consecutive minima ('low waters') varies, but its long-term average precisely matches the period belonging to the larger of the two sines. It is only for a narrow regime when amplitudes are very nearly equal that the long-term average can lie somewhere in between the two fundamental periods. It is for this reason that the M 2 period (being the strongest constituent at Vlieland Harbour) can be deduced very accurately from the tidalgauge data-without any recourse to astronomical considerations. We are not aware of any earlier mention of this remarkable fact.
From the above, it should already be clear that there is no unique definition of the tidal period. The question of which of the various options is the most suitable depends on the specifics of the problem and on the location under study. In our case, we focus on the question: what is the best definition to compute the residual flow in a multiple-inlet system?

Towards a basin-wide definition of the tidal period
The two causes underlying the difficulty of defining the tidal period (i.e., the presence of incommensurable tidal constituents and the effects of wind), as sketched in the previous section, are supplemented by a third one. Definitions of the tidal period based on local conditions at any one inlet do not represent the conditions in the multiple-inlet system as a whole since, for example, moments of low and high water vary from one inlet to the other. In the WDWS, the phase difference between the two major inlets, Texel and Vlie, amounts to 1.5 h. This lack of simultaneity needs to be taken into account if one wants to calculate the residual flow through the system as a whole. To address this last point, we first need to inquire into the nature of 'periodicity' in the present context.
The very usage of the word 'period' presumes an underlying periodicity, i.e., a return to the same state. But this begs the question as to what constitutes the state. For example, going from one low water to the next, the system shows periodicity only in the sense that the water level returns to a minimum. In other respects, the system may have changed; in particular, the water level in general changes from one low water to the next (e.g., due to diurnal inequalities and wind), and so does the volume of water in the basin. Similarly, the volume generally changes between alternate slacks. Over a tidal period thus defined, a change in freshwater content or sediment load will in part be simply due to the fact that the water volume inside the basin has changed over that period.
This suggests a way to define the tidal period in a global sense by taking the water volume as a reference, as proposed and used by Duran-Matute et al. (2014), and already described in the introduction. This definition reflects a periodicity in terms of volume: over one period, we return to the same volume. In the current paper, we compare the duration of the tidal periods and the magnitude of the tidally averaged residuals obtained with this definition to those obtained with other more classical definitions of the tidal period (e.g., time between alternate slack tides or mean sea-level crossings).

Model description
To study the effect of the definition of the tidal period for the calculation of the residual flow through a multipleinlet system, we will use results from numerical simulations of the hydrodynamics of the Dutch Wadden Sea already described in detail by Duran-Matute et al. (2014). These three-dimensional numerical simulations were performed using the General Estuarine Transport Model (GETM) with a horizontal resolution of 200 m and terrain following σ -coordinates with 30 layers in the vertical. Realistic meteorological forcing was provided by the German Weather Service (DWD; http://www.dwd.de) operational model with a spatial resolution of 1/16 • and a temporal resolution of 3 h and reconstructed times series of fresh-water discharges rates at 12 sluices were imposed. Bathymetric and freshwater discharges rates of the main sluices were provided by Rijkswaterstaat, part of the Dutch Ministry of Infrastructure and Environment. Fresh-water discharge data for smaller sluices were provided by the local water board responsible for the specific sluice.
For the current discussion, the boundary conditions at the open ocean boundaries are of particular importance. We have used Flather type boundary conditions that include realistic tidal elevations and storm surges, with their corresponding depth-integrated horizontal velocities. These data were obtained from a larger-scale operational model used by Rijkswaterstaat to predict the sea surface height along the Dutch coast. This model has a Kalman filter routine that is used to assimilate water level observations giving a realistic and accurate sea-surface height evolution (see, e.g., Plieger 1999). The computed sea level was compared with observations of sea-level height at 14 different tidal gauges in the numerical domain. For all stations, the coefficient of determination R 2 > 0.96 and the root-mean-square (rms) error was below 0.18 m. For the stations inside the WDWS, the rms error was below 0.12 m. A harmonic analysis was performed to compare every tidal component individually. For the dominant M2 tidal component, the relative error in the amplitude was below 10 % and the phase error below 15 min (2 %) for the stations inside the WDWS.
In addition, the results of the model were validated against several other observational datasets. They show a good agreement with salinity and temperature at one fixed station, and the instantaneous transport through the Texel Inlet as measured by the ferry crossing that inlet. For more details about the comparison with observations, we refer to the paper by Duran-Matute et al. (2014).
From the simulations, we extracted the current velocity data from transects across the first four inlets (Texel Inlet, Eierlandse Gat, Vlie Inlet, and Borndiep) and two watersheds (Terschelling watershed and Ameland watershed) as shown in Fig. 1. These transects, together with the mainland, define two basins: the WDWS (also known as the Marsdiep-Vlie basin and delimited by red lines in Fig. 1) and the Borndiep basin (delimited by purple lines in Fig. 1) at the East. The instantaneous water volume as a function of time for both of these basins and the instantaneous water transport across each transect have been computed with the model. As a convention, we consider the transport into the WDWS as positive and the transport outwards as negative.

Statistics on the tidal period
Several definitions of the tidal period can be used to calculate the residual transport through a tidal inlet. We first examine the effect that such a choice has on the statistics of the tidal period itself. For this, we consider the model data for the years 2009-2011 at the Texel inlet, one of the two  Figure 4 illustrates four of the definitions of the tidal period used in the current paper. These definitions are described in detail below.

Alternate mean sea-level crossings
As a first definition, we consider the tidal period to be the time between alternate moments when the sea-level, averaged over the transect, crosses the mean sea-level. Crossings occur during phases of rising (illustrated in red in Fig. 4) and falling tides. Hence, there are two possible definitions, depending on which phase is selected. The corresponding histograms are shown in Fig. 5a. It may happen during storm surges that the minimum level at low water stays above the mean level (or the opposite at extremely low levels: when the maximum level at high water remains below the mean). In those exceptional cases, the formally calculated period based on mean-level crossings includes several (maximum three) tidal periods. We then apply a correction by dividing the calculated time interval by the whole number of semi-days involved.
As expected, there is a large variability of the tidal period; moreover, the two options for the phases (crossings during rising tides or falling tides) yield different results. The histogram for the tidal period defined by mean sea-level crossings during falling tides is much wider (ranging from 9.5 to 15.5 h) than the one for rising tides (ranging from 10.5 to 15 h) as shown in Fig. 5a. This can also be seen in the values of the standard deviation shown in Table 1. The difference between the two histograms is due to the diurnal inequality, i.e., the fact that two consecutive tidal periods are unequal. A different phase of the diurnal tidal components is encompassed in consecutive tidal periods depending on the choice of the definition of the tidal period. This is further discussed in Section 5.5 and explained in Appendix A.

Alternate slack tides
The time between alternate slack tides as defined by zero instantaneous transport (averaged over the transect) is probably the most natural definition to calculate the tidally average flows. As with the previous definition, there are two distinct possibilities with regard to the phase: taking the slack tides from flood to ebb or from ebb to flood (illustrated in blue in Fig. 4). The choice between these possibilities has little effect on the shape of the histogram of the tidal period ( Fig. 5b), on the range of the tidal periods (from 11 to 14 h) or on its statistical measures (see Table 1).
The two definitions so far presented are local and have to be calculated for every transect individually. Remarkably, the histograms in Fig. 5a,b are quite different even though they deal with same tidal inlet and with the same years. This demonstrates how strongly the statistical properties of the tidal period depend on how it is defined. For example, both options based on alternate slack tides result in sharper histograms than those using alternate mean sea-level crossings.

Fig. 5
Histograms of the duration of the tidal period at the Texel Inlet. The tidal period is defined as the time between a mean sea-level crossings during rising tides (red/full) and during falling tides (blue/empty), b slack tides from ebb to flood (red/full) and from flood to ebb (blue/empty), and c mean-volume crossings during rising tides (red/full) and falling tides (blue/empty) a c b

Constant 12.42-h period
The crudest definition of the tidal period is a fixed 12.42h period since it completely disregards the variability in the system. (The probability distribution is a delta function, and for this reason, no histogram is presented.) This definition can be applied to the whole basin at once, as opposed to the previous two local definitions, so the calculation of the residual flows at the different inlets is now synchronized. This is a clear advantage in a multiple-inlet system. For this definition, the starting point is arbitrary. However, it is usually taken to be a slack tide. With this, however, the definition becomes local again, since slacks do not occur simultaneously in different inlets. This sacrifice is necessary to synchronize the residuals at the different inlets. As a starting point for the first tidal period, we consider here two options: the first slack tide from ebb to flood (illustrated in Fig. 4) and the first slack tide from flood to ebb, at the Texel Inlet.

Alternate mean-volume crossings
Finally, we define the tidal period as the time between two consecutive matches of the water volume inside the WDWS with the long-term average volume (see the area enclosed by red lines at inlets A, B, C, and watershed 1 in Fig. 1). Like the other definitions, this definition of the tidal period can be specified in two ways: mean-volume crossings during rising tides (flood) or during falling tides (ebb). The option of mean-volume crossings during rising tides, used also by Duran-Matute et al. (2014), is illustrated in green in Fig. 4.
Using alternate mean-volume crossings has the advantage of automatically producing a closed balance: all influxes and outflows must add up to zero for every tidal period. (Strictly speaking, it is mass, not volume, that is the conserved quantity, but the error is negligible here.) In addition, this definition is practical to compute the residual transport of freshwater or sediment, since their values cannot be distorted by a change in water volume inside the basin.
However, it has a similar shortcoming to the definition based on mean sea-level crossings: during storm surges, the minimum volume at low water can be higher than the mean volume, or in the opposite case, the maximum volume at high water can be lower than the mean. These are exceptional situations happening an average of 13 times per year, for which the method will include several tidal periods. We then apply a correction by dividing the calculated time interval by the integer number of semidiurnal cycles involved.
The two options of the definition based on mean-volume crossings produce nearly identical results (Fig. 5). Both show a large variability of the tidal period with a range from 10 to 15 h. Their corresponding histograms and the values of the standard deviation are similar to the ones obtained using the mean sea-level crossings during falling tides.

The effect of the diurnal inequality
The choice of the definition clearly affects the duration of the tidal period as reflected in their histograms (Fig. 5). The most noticeable differences lie in the width of the histogram (i.e., the standard deviation of the tidal period). This means that the choice can translate into a larger variability of the tidal period. This is due in a large degree to the diurnal inequality, i.e., the relative amplitude of the diurnal components and their phase lag with respect to the main semidiurnal component (see Appendix A for more details).
We have performed a harmonic analysis using T TIDE (Pawlowicz et al. 2002) on the different tidal signals considered to define the tidal period. The sum of O1 and K1 components have an amplitude of ∼28 % relative to that of the M 2 , for both the sea-level and the volume inside the WDWS. In contrast, this ratio is only 14 % for the transport. This already partly explains why the tidal period defined a b Fig. 6 Time series of the residual flow rate through the Texel Inlet for every tidal period during January 2009. a The residual transport is computed using the definitions for the tidal period based on slack tides from ebb to flood, or based on instants when the mean sea-level is passed during rising tides, or a constant 12.42-h period starting at a slack from ebb to flood, or matches of the volume inside the WDWS with its long-term average during rising tides. b Same, but now starting at the opposite phase: using the definitions for the tidal period based on slack tides from flood to ebb, or based on instants when the mean sea-level is passed during falling tides, or a constant 12.42-h period starting at a slack from flood to ebb, or matches of the volume inside the WDWS with its long-term average during falling tides using slack tides is less affected by the diurnal inequality, and hence, its distribution has a smaller standard deviation.
Actually, in terms of the amplitudes of the O1 and K1 components relative to that of the M 2 , the time series of the volume and that of the sea level are not very different (∼ 16 % for the O1 and ∼ 12 % for the K1). It is the phase of the K1 component that is responsible for the distinction between the histograms corresponding to tidal periods defined between crossings during rising tides and falling tides.

Statistics on the residual flows
We now address the question of how the definition of the tidal period affects the computation of the residual flow. For this, we computed the residual tidally averaged volume flow rate through the Texel Inlet by taking the average of the instantaneous volume flow rate during each tidal cycle for the years 2009-2011 using all the definitions of the tidal period described in the previous section. As an example, Fig. 6 shows the resulting time series for January 2009. While the residual flows based on 'local' definitions of the tidal period show an erratic behavior, as does the one based on a fixed 12.42-h period, the volume-based method produces a much smoother evolution. Although we only show the data for January 2009, this behavior is found to be consistent throughout the whole period studied.
There is a little difference between the tidally averaged volume flow rate computed using the definitions based on alternate slacks and constant 12.42 h. This suggests that, as usually thought, the error made by taking a constant 12.42 h starting at a slack is indeed small when compared to taking the actual instants of slack tides. However, the flow rates computed using both of these definitions show a large variability compared to the other ones. This is particularly true if the definition is based on slack tides from flood to ebb. Figure 7 shows the probability distributions of the tidally averaged flow through the Texel Inlet for every tidal period-defined using all of the options outlined beforefor the years 2009-2011. For the definitions based on alternate zero-crossings of the sea level, the distributions do not depend much on whether such crossings are taken during rising or during falling tides. For the definitions based on alternate slack tides, or on a constant 12.42-h period starting at a slack, the distributions highly depend on the phase of the starting point: from ebb to flood or from flood to ebb. The large variability observed in Fig. 6b, where the slacks from flood to ebb are selected, results in a much wider probability distribution in Fig. 7b,c. In fact, for these cases, the probability distribution loses its shape completely. For the definitions based on volume matches (Fig. 7d), the distributions are very similar and show no sensitivity to the choice of the starting point. Table 2 shows some statistical measures of the probability distributions of the tidally average flow through the Texel Inlet. All mean and median values are negative and hence The residual flow is calculated using the tidal period based on : a on alternate mean sea level crossings during falling tides (blue/empty) and rising tides (red/full); b alternate slack tides from flood to ebb (blue/empty) and from ebb to flood (red/full); c constant 12.42-h periods starting at the first slack tide from flood to ebb (blue/empty) and the first slack from ebb to flood (red/full); d alternate mean-volume crossings during falling tides (blue/empty) and rising tides (red/full) The small variations are due to the difference in the start of the first tidal period and end of the last tidal period. On the other hand, the median depends significantly-as much as a factor four-on the definition used for the tidal period. At the same time, all mean and median values are dwarfed by the corresponding standard deviation, which is one order of magnitude larger. This is a reflection of the great variability in the system. In practical terms, this means that measurements made on any one day are extremely unlikely to yield a residual flow anywhere near the typical (i.e., median) value. The volume-based definition of the tidal period, with its relatively low standard deviations, offers in this respect the best option. This is also reflected in the smoother time series of the tidally averaged flow rate computed using the volume based definition in Fig. 6. All distributions are positively skewed, with a longer tail and more extreme values on the positive (inward) direction. Duran-Matute et al. (2014) attributed this to the predominantly south-westerly winds in the region.
In Fig. 7b, c, the distributions using flood-to-ebb phases for the definition of the tidal period are particularly broad, resulting in a large standard deviation. This is in large part due to the diurnal inequality, as indicated in Fig. 6b. For these cases, we can see that the tidally average flow consists largely of up-and-down jumps from one semidiurnal period to the next. Interestingly, it is for these definitions that the tidal period itself has the least variation, while the flow rate varies the most. This suggests that the variation on the tidal period based on alternate sea-level or mean-volume crossings absorbs the diurnal inequality to provide a smoother time series of the volume flow rate through the inlets. In turn, this translates into a smaller dependence of the volume flow rate on the phase of the starting point.
The effect of the diurnal inequality, and hence, the variability in the computed residual transport would be greatly reduced if one considers two tidal periods instead of one. Further calculations indeed show that the time series of the tidally averaged flow computed using a constant 24.84 h definition is closer to the one obtained using the volumebased definition. A similar improvement occurs for the definitions using alternate slack tides or zero-crossings of sea level. However, lengthening the period over which the residual flow is considered comes at the cost of lowering the temporal resolution.
So far, we have considered the model data of three consecutive years together, as in Table 2. But how representative are these values (the median, for example) of each individual year, or of any other year? In Table 3, we present the results for the individual years. The mean residual flow varies strongly between the years, even in sign. The median is more robust and persistently negative (i.e., outflow), but still varies considerably. The implication of this All values are in m 3 s −1 , with the exception of the skewness which is dimensionless result is that one cannot in a meaningful way speak of a typical yearly-mean residual flow since it varies so much between individual years, but one may broadly indicate a typical range for median values. This variability must be due to the difference between the wind patterns in the respective years since the tidal constituents hardly vary from year to year and the averaged fresh water discharged varies only in the order of 100 m 3 /s between different years. Of the years studied, 2010 was the calmest in terms of wind and storms. For this year, the mean and median are fairly close, and the skewness is smaller. These values are indeed close to the ones found by Ridderinkhof (1988) in his model study that involved tides but no wind. This confirms a relatively weak contribution of the wind for that particular year, which is otherwise the main cause of variability.

Discussion
The upshot from the previous section is that the volumebased definition of the tidal period, while having itself a large standard deviation, actually renders the standard deviation of the residual flow relatively small compared with alternative definitions. This can be seen as a distinct advantage. The other advantage, of course, is that it is a globally defined quantity that applies equally and uniformly to the whole basin under consideration.
Nonetheless, two potential downsides must be considered. The first is whether the outcome for the residual flow depends on the choice of the basin (and basin size). This question becomes especially acute at watersheds, in our case the eastern boundary of the WDWS (see Fig. 1). Surely, the outcome for the residual flows should not depend on which basin, of the two at either side of the watershed, one chooses. This is verified in Section 7.1.
Second, the volume-based definition works with results from numerical simulations, where data of the water levels are available at all the (grid) positions in the basin and at every time step, hence the volume, too. But in the context of field measurements, such extensive information is lacking. However, if a certain number of tidal gauges happens to be available in the multiple-inlet system, one may construct a time-series of the volume by interpolating between the stations. We will here test that method and will check how close the tidal periods thus obtained are to those from the model. All values are given in m 3 s −1 , except for the skewness which is dimensionless. Here, the tidal period is defined using zero mean-volume crossings during rising tides

Dependence on demarcation of the basin and the problems at a watershed
The watershed between two basins presents a number of challenges and particularities. For instance, the complex bathymetry, composed of many channels, gullies and extensive tidal flats, makes it very difficult to measure the transport across the watershed. The watershed also presents a dilemma for the volume-based definition: it is the boundary between two different basins, and hence, either one could be chosen to define the tidal period. The question to be addressed here is whether that choice would affect the calculated residual transports at the watershed. Figure 8 shows the instantaneous transport through the Terschelling watershed during January 2009 as obtained from the numerical simulations. It is much unlike a sinusoidal signal. This is caused by nonlinear effects, the clearest of which is the effect of shallowness: whenever high waters are lower than the height of the tidal flats, the transport is confined to the deepest channels, but when high waters are higher, the transport takes place over the whole cross-section. Around low waters, there may be even occasions without any transport (some instances can be seen in Fig. 8). Due to these irregularities, it is hard to define the tidal period based on local slack tides or mean sea-level crossings.
In Fig. 9, we show the histograms of the tidal period and of the residual flow rate, comparing in each case two options: using either the western basin (delimited by red lines in Fig. 1) or the eastern basin (delimited by purple lines in Fig. 1) to define the volume-based tidal period. In general, the differences are slight, which demonstrates the robustness of the method.

Empirical application of volume-based definition
In previous sections, we discussed several ways of defining the tidal period and argued that the volume-based method seems the most appropriate one in a multiple-inlet system. However, there is one practical drawback to it: while the other methods depend on local data (a tidal gauge or current velocity data), the volume-based method requires information on the entire basin for a long stretch of time. In a numerical model, the volume can be readily calculated at any time step, but in the field such data is, of course, impossible to collect.
However, an estimate of the volume may be made indirectly if a sufficient number of tidal gauges and an accurate bathymetric map are available. The idea is to interpolate between water levels simultaneously recorded at the tidalgauge stations. In the Dutch Wadden Sea, these records have a sampling rate of once every 10 min. At any given position, an interpolation may not be expected to give a very reliable representation of the local water level, but for bulk quantities like the volume, one may plausibly expect that inaccuracies mostly cancel out. We will here test how good the method works.
We split up the WDWS into five triangles (Fig. 10), which together form a basin that covers most, but not all, of the basin defined in Fig. 1. So, the volume calculated here will be somewhat smaller than the one from the model, but the purpose really is to identify the moments of meanvolume crossing, and from that, the tidal period and its long-term distribution. From crossings of the mean volume at rising tides, we can define the tidal period as the time between successive crossings, as explained before. The tidal periods, calculated with identically applied methods, based on the model and on tidal gauges can now be compared: we find that the mean absolute difference is 7.1 min, which falls well within the sampling rate of the tidal-gauge records; the linear correlation coefficient is 0.99.
Finally, we calculate the respective distributions, see Fig. 11. It turns out that they are quite similar. In other words, in areas with a good coverage of tidal gauges, such as in the WDWS, the volume-method can be applied even if no model results are available. This greatly enlarges the scope of applicability of the this definition of the tidal period.

Integration vs. subtraction or filtering of tides
Two methods commonly used to extract the subtidal motions out of both sea-level height and currents records are (1) obtaining the tidal constituents by harmonic analysis and subtracting them from the full signal (see, e.g., Horsburgh and Wilson 2007;Guyondet and Koutitonsky 2008) and (2) filtering the tides from the full signal (see, e.g., Chant 2001;Buijsman and Ridderinkhof 2007). The data that results from applying either of these methods is a time series of the subtidal signal with the same temporal resolution as the original time series. This contrasts with the result obtained by integrating over every tidal period-i.e., one value per tidal cycle. Hence, the choice of integration over subtraction or filtering of tides depends on the overall aim of the analysis. For example, filtering out the tides from the sea surface height signal is very useful for the understanding of the causes of peak sea surface elevations to plan coastal protection (Brown et al. 2012). Integration over a tidal period is less useful for this, as the relevant time scales can be shorter than the tidal period.
Both methods-subtraction and filtering of the tideshave been studied, developed, and compared extensively (Jay and Flinchem 1999;Brown et al. 2012). Hence, most of their shortcomings are well known. For example, computing the tides using harmonic analysis to later subtract them from the total signal is particularly challenging in shallow areas like the Wadden Sea. Harmonic analysis rests on the assumption that the tide is a superposition of the tidal harmonics with constant phase and amplitude. However, strong nonlinear interactions that modify temporarily the phase and amplitude of the individual tidal components can occur in shelf seas during storm surges (Horsburgh and Wilson 2007), due to stratification , and due to the interaction with the bathymetry (Friedrichs and Aubrey 1988). Moreover, as mentioned by Brown et al. (2012), residuals are typically small compared to the full signal, and hence, they are very sensitive to inaccuracies resulting from the harmonic analysis. This was borne out by a recent comparison between a direct calculation of the residual current from model data in a tidal inlet (using the 'volume method') analysed in the current paper) and the residual calculated by harmonic analysis from the same data; the discrepancy turned out to be extremely large (Sassi et al. 2015).
For extreme cases of nonlinear interactions and irregular signals, both harmonic analysis and filtering break down. It is clear that it is practically impossible to get meaningful results from applying either of these methods to a signal as irregular as the one for the volume flow rate through the Terschelling watershed (Fig. 8). For this reason, it is practical to use a basin wide definition to have a direct way of calculating the residual through-flow at least once per tidal cycle.

Conclusions
We have studied the residual flow over the tidal periods within 2009-2011, comparing different possible definitions of the tidal period. One definition is based on the volume of water inside the basin and thus holds basin-wide: the same tidal period applies synchronously to all points in the basin. This stands in contrast to other more common definitions that are based on local dynamics, notably when the tidal period is defined as the time between alternate slacks.
The basin-wide validity is not the only advantage of the 'volume-method'. We also found that it is robust with respect to the choice of the starting point. The method uses the moments when the actual volume equals its long-term mean. There are two kinds: the moment can coincide with rising water levels or with falling ones. Either can be used for the definition of the tidal period, but our results indicate that neither the distribution of the tidal periods nor the distribution of the residual flows is significantly influenced by those choices. In this sense, the method is robust. This too stands in contrast with using slacks, where the distributions differ a great deal depending on whether one selects slack from ebb to flood or from flood to ebb to define the tidal period. Between the different ways of defining the tidal period, the resulting long-term median residual flow differs considerably (Table 2), but the mean is nearly the same for all.
It is, however, difficult to attach meaning to a yearly mean value of the residual flow through an inlet beyond its validity for that specific year. As summarized in Table 3, the yearly mean residual flow varies strongly between the years, even in sign. This is probably due to changing wind conditions between years since wind can have a large effect in the residual circulation in multiple inlet systems (Smith 1990;Li 2013;Duran-Matute et al. 2014). This means that decadal time-series would be needed to speak in a statistically meaningful way of a long-term mean residual flow. But in the course of such time-scales, the morphology of the basin itself may have changed, rendering the notion again doubtful. Similar issues are presumably at stake in the transport of suspended sediment.
Physically more meaningful, it seems, is to consider shorter periods in which certain conditions prevail (e.g., days with easterly winds) and to study the dynamics and residual flows through the system in that setting. Then, the question of how to define the tidal period becomes acute, and this paper offers ways to deal with them.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A: Explanation on the effect of the diurnal component
The effect of the diurnal component on the distribution of the tidal period can be already seen by considering a simplified tidal signal f (signifying sea-surface height or transport) with one semidiurnal and one diurnal component: where the frequency of the semidiurnal component is 2ω, and its amplitude has been normalized; φ is the phase lag between the semidiurnal and the diurnal components, and C is the amplitude of the diurnal component. We are interested in consecutive zero-crossings, i.e., in the solution of the equation f (t) = 0. This equation has, in general, no analytical solution except for the cases when φ = nπ/2 with n = 0, 1, .... These cases already provide an important insight. For φ = 0, the zero-crossings occur for ωt = 0, π − arccos(C/2), π, π + arccos(C/2), 2π etc.
(1) This means that if the tidal period is defined as the time between alternate zero-crossings, starting at t = 0 (phase of rising tides), then they are all equal to π/ω. On the other hand, if we start at t = [π − arccos(C/2)]/ω (phase of falling tides), then two different tidal periods occur: one alternately finds a shorter period 2 arccos(C/2)/ω and a longer period [2π −2 arccos(C/2)]/ω. Notice that these differences disappear when the diurnal component is omitted: C = 0 renders all periods equal to π/ω, irrespective of the starting phase.
In this simple example, the probability distributions are given by one or two delta functions. In the former case, there is one delta function at the semidiurnal frequency and the standard deviation is hence equal to zero. In the latter case, there are two delta functions, but the average of the overall signal remains equal to the semidiurnal frequency. On the other hand, the standard deviation is equal to |[2 arccos(C/2) − π]/ω|. This shows that the diurnal component can affect the shape of the distribution and can create a distinction, depending on the phase at which one starts to define the tidal cycle.
For general values of φ, not only zero-crossings during falling tides can produce alternate tidal periods, but also those during rising tides. Because of the periodicity of the signal, there can be only two different tidal periods involved in each. Each distribution thus has just two peaks. A special case occurs when the phase difference is equal to π/4 and the distributions for the tidal periods defined by crossing during both rising or falling tides are equal. In this case, the diurnal inequality increases the standard deviation of the distribution of the tidal period, but does not produce a difference between the distributions obtained by considering crossing during rising or falling tides.
To get a continuous distribution of tidal periods, yet another element has to be introduced. In the previous example, the ratio of semidiurnal to diurnal frequencies was chosen to be exactly 2, for simplicity. If we relax this (unnatural) assumption and introduce irrational fractions (e.g., by taking the semidiurnal constituent to be M 2 , but the diurnal constituent K 1 ), then the resulting distribution gets smeared out. Because of the presence of a diurnal component, one still finds different histograms depending on the phase at which one starts to define the tidal cycle. Thus, with even these very limited ingredients, some key features of the histograms presented in previous sections can already be elucidated.
This can be seen by performing harmonic analysis in a real tidal signal and plotting the histograms obtained with the superposition of the different tidal constituents. In this case, we consider the sea-surface elevation at the Texel inlet. The harmonic analysis was performed using T TIDE (Pawlowicz et al. 2002), and the seven most important tidal constituents with their respective amplitude and phase are given in Table 4. Figure 12 shows the histograms of the tidal period of the 3 years simulated considering: (a) the seven components excluding the diurnal components, (b) the seven components excluding the O1 component, (c) all seven most important components. The histogram of the tidal period without diurnal constituents shows a distinct peak at the semidiurnal tidal period with a small standard deviation (8 min for tidal periods between consecutive crossings during rising tides, and 7 min for tidal periods between consecutive crossings during ebb). When the first tidal period (h) all seven most important tidal constituents. In filled red: tidal periods between consecutive crossings during flood. In empty blue: tidal periods between consecutive periods during ebb diurnal component (O1) is introduced, the standard deviation increases about a factor four to 32 and 26 min, and the distribution obtains a clear binormal shape for both definitions. The O1 tidal components is 34 • out of phase with the M2 tidal constituent, and hence, the results of both definitions-crossings during flood or during ebbare similarly affected. In contrast, the K1 constituent is out of phase with the M2 by 180 • , and its inclusion creates a large difference between the resulting histograms of the two definitions.