## Abstract

Persistent reports exist that the tides in coastal basins are often accompanied by regular or irregular oscillations of periods ranging from minutes up to several hours. A conceptual model relating the two is proposed here. It employs an almost-enclosed basin, connected to an open (tidal) sea by a narrow strait. Such a basin is a Helmholtz resonator, which is dominated by the “pumping” mode. Its response is governed by an ordinary differential equation that is forced by the tide, damped by friction and wave radiation, and whose restoring term is nonlinear due to the sloping bottom. When forced resonantly by a single frequency tide, due to this nonlinearity, the basin may exhibit multiple equilibria. Its response can either be amplified or choked, depending on the precise initial conditions. The presence of a second forcing term may, on a slow timescale, kick the system irregularly from the amplified into the choked regime, yielding a chaotic response. This may happen when either two nearby frequencies, for example, a combination of semidiurnal lunar and solar tides, are near resonance (and the frequency difference provides the beat), or when a small-amplitude, resonant perturbation is modulated by a large-amplitude, low-frequency tide. The aforementioned observations of irregular tides are discussed in the light of analytical and numerical results obtained with this model for these two regimes.

## 1. Introduction

“Chaotic tides” sounds like an oxymoron. As the etymology of the word “tide,” being time, season, or period according to Webster's (1913) *Revised Unabridged Dictionary* implies, it has, from ancient times on (Pugh 1987), often been held to be synonymous with periodicity. Indeed, this regularity made tidal heights one of the earliest known and best predictable phenomena in physical oceanography, with relative accuracies of prediction often exceeding 90%. The predictive capacities of the corresponding tidal *currents* is, admittedly, less impressive, but is still often in excess of 50%. Less known, however, is the fact that spurious, but persistent, reports on the *irregularity* of the tides exist for over a century (Honda et al. 1908; Krümmel 1911). It is this irregularity in the tide that we here aim to identify as being possibly due to the chaotic response of certain coastal basins to the tide at the open sea.

Tidal predictions have traditionally been made for individual ports on the basis of previously observed tidal elevation records at those locations. Amplitudes and phases of a number of precisely known “tidal frequencies,” stemming from the gravitational potential determined by celestial mechanics, are estimated, which can then be used to sum the corresponding harmonic series and thus to predict the future occurrence of tides at that location; see, for example, the review by Godin (1991).

Although the gravitational (or tidal) potential carries only about 5 relevant, independent frequency components, which, through an expansion of the tidal potential in harmonic terms leads to perhaps 11 principal tidal components (Platzman 1971), in the shallow, coastal areas, local nonlinear effects lead to leakage of energy to higher harmonics and sum and difference (combination) frequencies. In this way tidal energy fills in entire spectral bands surrounding the principal components (Pugh 1987, p. 188). As there are, in principle, no difficulties in also resolving the amplitudes and phases of these combination frequencies, running harmonic analysis with 150 frequency components or more has become routine procedure. This empirical *local* analysis, however, does not automatically guarantee that the amplitudes and phases of these independently determined frequencies also show a *spatially* coherent and stationary picture. Quite the contrary.

A spatially coherent description of the tides can be obtained by solving the governing *dynamical* equations, the Laplace tidal equations (see, e.g., Cartwright 1977). In a realistic ocean domain this is achieved by numerical methods (Schwiderski 1980). Outside continental margins, tidal predictions for these principal components, particularly when constrained by observations from deep-sea tidal stations, are in fair agreement with observations obtained at different seaports and deep sea tidal stations (global root-mean-square difference between modeled and observed elevations of about 10 cm). With the advent of satellite altimetry providing tidal observations, this agreement has in recent years actually been dramatically improved (global root-mean-square difference of about 3 cm); see Andersen et al. (1995).

Except for these directly driven, principal components, together perhaps with their first few harmonics (Davies 1986; Lynch and Werner 1991), spatial coherency is cumbersome and agreement with locally observed tidal components poor, particularly in the shallow, coastal areas. Waves at combination frequencies are generated locally, get amplified, or damped, to a degree depending on the particular resonance properties of that locality, and do not persist as free waves so that they decorrelate quickly. Treating the nonlinearly generated compound waves locally as linearly independent (as harmonic, or Fourier analysis wants it) is therefore somewhat contorted but still defendable when the local nonlinear transfer of energy is *frozen* so that harmonics appear as the result of a stationary process. However, when analyzing different, independent datasets from the same location, variations in tidal amplitudes and phases often still occur, to the extent that some of the minor components appear, in fact, *unresolvable* (Gutiérrez et al. 1981; Godin 1991). Although such variations in locally estimated tidal amplitudes and phases are usually attributed to nonstationary effects due to wind (van Ette and Schoemaker 1966; Gutiérrez et al. 1981), this need not necessarily be its only source. Nonstationarity may also be due to nonlinearity of the hydrodynamic system itself (Pugh 1987), which may not only change the tidal elevation profile, by giving rise to *superharmonics* (overtides, which stay fixed in time), but may also modulate its amplitude, by giving rise to *subharmonics.* Finally, strong nonlinearity may provoke a cascade of such subharmonic bifurcations, giving rise to *chaotic behavior.*

Is there any observational basis for chaotic behavior of the tides, other than the aforementioned unresolvability, or noisiness, of tidal “constants”? In an analysis of the tides in Venice Lagoon, at the head of the Adriatic Sea, where the tides seem to pick up because of near-resonancy of the basin, Vittori (1992) observed that consecutive tidal maxima are highly irregular. She argued this to be indicative of low-dimensional chaos. Whether the low-order dynamics to which this is due is either inherited from the dynamics of the local wind fields or of a genuinely oceanographic nature is not clear. Similar changes in consecutive maxima also appear in long-term, tidal elevation records in the Wadden Sea basins (J. T. F. Zimmerman 2000, personal communication), which are again close to resonance (see Maas 1997, hereafter referred to as M). Frison et al. (1999) find evidence of nonlinear tides in U.S. coastal estuaries from the nontrivial shape of the attractor obtained in a diagram where observed elevation at some moment is plotted against that observed at some previous instant.

Surprisingly, direct observations of irregular oscillations, accompanying the tides, have been reported for over a century (see Krümmel 1911, 157–185; Defant 1961, p. 187). The suggestion that these “secondary oscillations” were actually related to the “primary,” tidal oscillation was firmly put forward in an impressive study by Honda et al. (1908), who had intensively studied the seas around Japan. This study eventually led to a complete classification of bays in terms of “periodic,” “quasiperiodic,” etc. (Honda et al. 1908; Nakano, 1932). In the absence of a dynamical framework, such as recently offered by the field of nonlinear dynamics, this classification was, however, not further interpreted, and these observations seem, if not completely forgotten, neglected in contemporary tidal literature. One reason for this might be that these secondary oscillations generally seem inconspicuous, as they consist of high-frequency, low-amplitude waves superimposed on a low-frequency, high-amplitude (primary) tide. Typical periods range from some minutes to several hours, and amplitudes do not exceed a few percent of the tidal range. However, an interesting twist is given to this interpretation by recent observations of such irregular tidal elevations in a Norwegian fjord (Golmen et al. 1994). By making simultaneous observations of the currents in the strait connecting the fjord to the sea they found that these irregular small-amplitude elevations were accompanied by irregular *O*(1) variations of the tidal current (see Fig. 1). As velocities are determined by elevation *differences* over the strait, by inference the *velocity* field can be amplified when, at any moment, the difference in tidal elevation also amounts to just a few percent of the elevation itself, thereby becoming of comparable magnitude to the difference between the (small amplitude) high-frequency elevation that is resonantly excited within the basin but is practically absent at sea. Such *O*(1) velocity variations are not only important for nautical reasons, but are clearly equally relevant for the flushing of the fjord: the exchange of water, sediments, and dissolved gases or nutrients with the connected sea. It is believed (LeBlond 1991) that this property will transcend the global relevance that both dissipation as well as resonance of tides in the coastal environment may have in setting the boundary conditions for the global tide (Garrett 1975), although Munk and Wunsch (1998) recently made the interesting suggestion that tides might be playing a key role in ocean circulation.

In order to illustrate these issues we will consider the resonant response of a short, deep basin having a sloping bottom, to the tide in the adjacent open sea to which it is connected by a narrow strait—a nonlinear Helmholtz resonator (see M). This choice of geometry differs from the mostly wide bays with their quarter-wave resonancies, encountered by Honda et al. (1908), but it ensures that the tide within the basin takes the simplest form possible, as it is governed by the pumping, or Helmholtz mode. This mode is the lowest and generally also the most important frequency in the basin's spectrum, and is characterized by a spatially uniform response. It can thus be represented by a single state variable (the excess volume of water, related to the free surface elevation), whose evolution is therefore described simply by an *ordinary* differential equation. Although nonlinearity can also be present in the frictional damping term (Zimmerman 1992), it is the nonlinearity in the restoring term that gives rise to multiple equilibria and chaos here. The latter nonlinearity is due to the sloping bottom (Green 1992) and can be understood as follows. The current through the strait is driven simply by the elevation difference over the strait. However, the elevation change within the basin, affected by the transport through the strait, clearly depends on the surface area that the incoming water has to cover, which is nonuniform with depth when the walls of the basin are not vertical. The time needed to produce a particular elevation change therefore depends on the preexisting sea level; see Fig. 2. In M its free, forced, and forced-and-damped response was discussed in the case where the tidal forcing is “pure,” that is, contains one frequency component only. In inviscid circumstances, a single frequency forcing can provoke a chaotic response in a small parameter range. However, when damping is added, the chaos fails to persist. Since the inclusion of damping terms is mandatory in any realistic physical context, this “Hamiltonian chaos” appears to have no physical relevance. It is one of the goals of the present paper to determine under what conditions tides, within basins of the kind considered here, can become chaotic, even in the presence of damping. Multiple frequency forcing seems to be a prerequisite for this. Therefore, in light of the observations discussed above, we will consider forcing with either two nearby frequencies, as is the case for a combination of lunar and solar tides, a “mixed” type of tide (see Defant 1961), or with a single low-frequency tide, accompanied by a small-amplitude, high-frequency, resonant perturbation, as is the case in a fjord.

In section 2 the nonlinear tidal Helmholtz resonator will be introduced. It represents the response of a bay to tidal variations at the open sea, to which it is connected by a narrow channel. Its governing ordinary differential equation will be derived. This section also introduces the Poincaré phase plane, which gives a more comprehensive way to follow its evolution than that obtained by direct numerical integration. The tidal response of the bay typically shows a gradual modulation of amplitude and phase, which can be captured by averaging the original governing equation (section 3). In section 4 this modulation equation is applied to the case that forcing is at two nearby frequencies. When these are close to resonance, the response of the basin may, under certain conditions, be chaotic. The case that a single, low-frequency tide is accompanied by a small-amplitude, high-frequency resonant perturbation is very similar to this, and is discussed in section 5. In section 6 the results will be summarized.

## 2. Nonlinear Helmholtz resonator

An almost-enclosed, short basin, connected to a tidal sea by a narrow strait (modeled here as a *pipe*), is generally governed by the pumping (Helmholtz) mode, in which the water level within the basin rises and sinks in unison. Elevation is scaled with (maximum) depth *H,* volume with *A*_{0}*H* (where *A*_{0} is the surface area of the basin at rest), while time *t* is scaled with the Helmholtz frequency of the basin *σ*_{H},

with *B* and *L* denoting the strait's width and length, respectively, and *g* the acceleration of gravity. In M it is shown that the evolution of the (nondimensional) amount of water in excess of that present in the basin in the absence of tides—the *excess* volume, *υ*—is (neglecting quadratic damping terms) governed by

Here *ζ* and *ζ*_{e} denote the surface elevation within the basin and in the exterior sea, respectively. The latter is a prescribed function of time. At any time their difference provides the pressure difference that drives the flow *u* along the connecting strait: *du*/*dt* = *ζ*_{e} − *ζ.* This flow, multiplied by the (unit) cross-sectional area of the strait, equals the rate at which the volume within the basin changes, *u* = *dυ*/*dt.* From these two equations the inviscid form (*c* = 0) of (2) follows. The basin hypsometry (horizontal area *A* as a function of vertical coordinate *z*) determines the excess volume within the basin. With respect to mean sea level *z* = 0, this is defined as

The inverse of this relation, needed in (2), yields the nonlinear restoring term *ζ*(*υ*), which can be made explicit for some simple basin shapes. For a basin with vertical sidewalls, *A* = 1, for instance, it follows from (3) that the restoring term is linear, *ζ* = *υ.* However, when the basin area increases linearly with depth, *A* = 1 + *z* (the deepest point of the basin being at *z* = −1), this restoring term is given by

This shape of the basin is the only one addressed here and its hypsometric behavior is thought to be characteristic of a typical basin with shoaling sides (see M). Other basin shapes and the effect of this shape on the dynamics of the basin are addressed elsewhere (Doelman et al. 2001, submitted to *Phys. D,* hereafter DKM). A nonempty basin (and analyticity of the solution) requires *υ* > −1/2. In the right hand of (2) we find forcing by the prescribed external tide, *ζ*_{e}(*t*), and damping, respectively. Although this linear damping formally stands for radiative damping only, we may also consider it to incorporate a linearized version of bottom friction and form drag (see M). For further analysis we treat (2) as a dynamical system by writing for the basin's excess volume *x* = *υ,* and for the strait's flow rate *y* = *dυ*/*dt,* and obtain, with (4),

As usual, a dot indicates a time derivative. Although the emphasis will here be on the response to quasiperiodic (two frequency) forcing, some results obtained without forcing and with “single frequency” forcing, pertinent to the discussion to follow, will be recapitulated first.

### a. Poincaré plane, free response, and periodic forcing

A numerical example of an integration of (5) is given in Fig. 3a for the case without forcing. It presents the evolution over two tidal periods of volume and current, (*x*(*t*), *y*(*t*)), as a function of time, increasing upward. Projections of this helical curve at the back (*x*–*t*) and side (*y*–*t*) planes give the time evolution of volume and current proper, such as they are commonly presented (except that time now points upward). The third projection, on the bottom *x*–*y* plane, suppresses time, which now acts as a parameter *along* such orbits, and is called the phase plane. In this strongly viscous example time can still be discerned in the gradual decay of the oscillations. However, when damping is also absent, solution curves will “circle” around a common center point, which is the only fixed point present in this case. Each solution curve is egg shaped and is described by elliptic functions in terms of which the inviscid, unforced version of (5) can be solved (M). Different curves represent different initial conditions (volume and current). The phase plane is limited by an outer boundary (the limiting orbit) because the argument of the square root in the restoring term of (5) needs to be nonnegative. It can be found in parametric form as

for *t* ∈ (−3, 3), which is periodic with period 6. Owing to the nonlinearity of the oscillator, the period of such a free, finite-amplitude oscillation is a function of its amplitude. For the present geometry this period spans just a narrow range, decreasing from 2*π,* at the center, down to 6 at the outer boundary. The frequency range of other geometries is considered in DKM.

A comprehensive way of monitoring gradual changes of the tidal response due to forcing, damping, or nonlinearity is obtained by determining subsequent intersections of the helical curve with a sequence of periodically placed, horizontal planes at *t*/*T* = 0, 1, 2, … , with *T* indicating the forcing period (so-called Poincaré planes of section). These subsequent intersections define a map, called the Poincaré map, 𝒫. A strictly periodic signal (of period *T*) thus yields the same intersection. Compiling each of these subsequent intersections on one and the same horizontal (phase) plane, a stroboscopic impression of the slow evolution of the nearly periodic behavior is obtained. An example of such a Poincaré plane is presented by the dots in Fig. 3b. Here, the spiral is the projection of the continuous helical curve of Fig. 3a on the bottom plane. The phase at which the forcing is sampled can still be varied but is in this study taken to correspond to the moment that the external tide is at low water. From here on we will mostly use this phase-plane representation of the evolution, accepting the suppression of time information. This is illustrated in Fig. 4 with some stroboscopically sampled solutions of (5). When forcing is absent, damping brings any initial motion to rest, in this Poincaré plane represented by a sequence of points approaching the origin (Fig. 4a). When damping is absent and a single forcing component,

is near resonance (*ω* ≈ 1), the amplitude of the basin tide will grow because of the nonlinearity of the restoring term. However, doing so, it changes its period so that the resonator is detuning itself, and the tide at sea will actually damp the basin response. Once it is small again, this process may start afresh, and the whole cycle takes a time span, which one calls the modulation period, that is large compared to the tidal period. This slow growth and decay thus reveals itself in the Poincaré plane as a single closed loop that is traversed in one modulation period. Other confocal loops represent other amplitudes of oscillation. However, when the forcing period is close to resonance, it is clear from the three curves of Fig. 4b that there may be different types of responses, which will be discussed further in the next section. With forcing and relatively strong damping, a regular oscillation of definite amplitude and phase eventually results, so stroboscopic intersections will approach one single spot, away from the orgin (Fig. 4c). When damping is weak and forcing is near resonance, depending on the precise initial state, the system may lock into *different* states (Fig. 4d).

To summarize, in this Poincaré plane, a strictly periodic oscillation will be represented by a single (fixed) point, a quasiperiodic oscillation by multiple points lying on a closed curve (limit cycle), and an aperiodic state (encountered in later sections) by an irregularly located set of points (strange attractor).

### b. Qualitative effect of quasiperiodic forcing

The multiple equilibria near the Helmholtz resonance correspond to an amplified and a choked tidal response (see the two centers in Fig. 4d, at large and small volumes |*x*|, respectively). These states are stable and attracting when damping is present (Fig. 4d). Without damping, however, these states would retain their original distances to the fixed points (as the curves in Fig. 4b), and each state can be characterized by its energy level and period (with which it is traversed). The qualitatively different regimes, recognizable in Fig. 4b, are separated from each other by particular special orbits, the separatrices. One outer separatrix separates the outermost orbit from the other two, and one inner separatrix separates the banana-shaped orbit and small circular orbit. The period needed to traverse these separatrices approaches infinity, as the intersection of the two separatrices acts like a saddle. A state can creep up to this relative maximum in potential energy by exchanging some of its “kinetic” energy, allowing it to move forward along its orbit, in favor of “potential” energy.

When friction is introduced, it pulls the state down into one of the equilibria, except for the odd point sitting right at the saddle and the two orbits that limit exactly on the saddle point in forward time. Note that these orbits merge with the aforementioned separatrices as the friction decreases to 0. Apart from these exceptions, all initial states thus belong to one of two domains of attraction (amplification or choking). When one starts to weakly and slowly “shake” the system (i.e., the saddle and the separatrices) by introducing a large-period perturbation, a new &ldquo=uilibrium” may arise. This happens because it is clearly quite a sensitive matter near the separatrices to which equilibrium their state will recede. Therefore, when, due to shaking, the separatrices move a little, while the state is slowly receding to one equilibrium, it may find itself a little later in the attraction domain of the other. Making again its way to the other equilibrium, something similar may happen, and, in effect, the neighborhoods of the separatrices may trap the state.

This new “steady” state might turn out, however, to be chaotic. This is because, close to the separatrices, the period of the inviscid orbits (the modulation period of the Helmholtz oscillator) varies greatly, and thus, when the state of the system varies due to the weak forcing, neighboring states fastly diverge from each other.

The way to verify that such a new dynamical equilibrium state turns into existence under the addition of perturbative damping and extra forcing is by performing an energy budget study. For this, one picks one of the original inviscid orbits, characterized by a particular period, and one then evaluates, along that orbit, the net increase or decrease of energy over one period due to both the work done by the perturbative forcing as well as the energy dissipation by the perturbative damping term. Only when these two are in balance (and there is thus no net energy increase or decrease) will the point on average stay near the original, unperturbed orbit. The single orbits considered in the present study are the previously introduced separatrices (homoclinic orbits), along which the period turns to infinity. The function determining the net energy change on the homoclinic orbits is called the Melnikov function. This function is usually given the alternative interpretation of measuring the distance of the perturbed orbits leaving and approaching the saddle. A (traditional) search for its zeros—implying intersections of the two perturbed orbits so that points both belong to the sets of points approaching and leaving the (perturbed) saddle, and are therefore “stuck”—is thus equivalent to the requirement that there is no net change in energy content along its trajectory. Mathematically, all we can show for the moment is that certain points (that belong to this invariant set) get trapped near the unperturbed separatrices, and that they traverse that region in an unpredictable, chaotic way (section 4a). However, we lack the ability to rigorously show that this set is actually *attracting* (Wiggins 1990): Suggestions that it is follow from numerical analysis (section 4b).

Similar such phenomena appear at other resonances, like near *ω* ≈ 1/2 + ɛ, and in particular near *ω* ≈ 2 ± ɛ (see M).

## 3. Modulation equations

Much of the behavior obtained by numerically integrating the exact equation can, to a remarkable degree, also be found analytically in the solutions of the modulation equations that will be considered now. The advantage of using modulation equations is that fixed points, homoclinic orbits, Melnikov functions, etc. can be obtained explicitly, which is very helpful in addressing the main issue, namely, whether (and under what conditions) chaos can occur in the nonlinear Helmholtz resonator in the presence of damping. The stroboscopic intersections in Fig. 4 trace out continuous curves, whose locations can be obtained approximately when averaging the governing equations over the forcing period. Doing so, one obtains modulation equations that govern the amplitude dynamics on a slow timescale, which reveal where and under what conditions the dynamical response can become “complex.” When the excess-volume fluctuation *υ*(*t*) is of small amplitude (and when also the forcing and frictional parameters are small), (2) can be approximated by rescaling *υ* = ɛ*V, **c* = ɛ^{2}*C, **f* = ɛ^{3}*F* (with ɛ ≪ 1). With sinusoidal forcing (6), this leads to

This equation is solved by a multiple-scale perturbation expansion and has an amplified response at rational values of *ω* (Nayfeh and Mook 1979, p. 196). By requiring the absence of secular terms in the equations governing the subsequent orders of the perturbation expansion near the primary resonance (which, for consistency, requires *ω* = 1 + ɛ^{2}*σ*) the solution is, up to second order, given by

The response is thus frequency locked to the forcing and exhibits a “drift” term (the constant term), providing a net displacement of the mean state. Amplitude *R* and phase Φ vary on a slow timescale and are governed by

where a dot indicates differentiating with respect to the slow time variable *T* = ɛ^{2}*t.* Equating the right-hand sides to zero determines the steady states and allows one to obtain the frequency response curve for these equations *σ*(*R*; *F, **C*); see Fig. 5. However, by rescaling *t* → *t*/*σ, **R* → *σ*^{1/2}*R, **C* → *σC,* and *F* → *σ*^{3/2}*F,* the detuning frequency *σ* can be scaled out and can, therefore, without loss of generality, be set equal to one. We may use this property to argue that the actual detuning *defines* the small parameter ɛ ≡ (*ω* − 1)^{1/2}, which would otherwise be left arbitrary. This is consistent with the two rescalings, employed previously. In Cartesian coordinates, *X* = *R* cosΦ and *Y* = *R* sinΦ, these modulation equations read

where *R*^{2} = *X*^{2} + *Y*^{2}.

For the single-frequency forcing one can assume, without loss of generality, *θ* = 0 (which amounts to a rotation of our coordinate system). Its explicit presence is retained here however for later use, when, in the double-frequency forcing case, both *F* and *θ* will be varying on the slow timescale. Note that these modulation equations are the same as those found for the Duffing equation [in polar coordinates by the method of multiple scales, see Nayfeh and Mook (1979), and in both polar and Cartesian coordinates by the averaging method, see Guckenheimer and Holmes (1983)]. The extension to two forcing frequencies, discussed later, should thus have relevance to the quasiperiodically forced Duffing equation (Wiggins 1990; Yagasaki 1990, 1993), and extends it to the “double resonant” case of two nearby frequencies, both approximately equal to the natural frequency.

### a. Comparison of exact and modulated systems

It is useful to understand the system of modulation equations because this *local* approximation has a strong correspondence with the Poincaré plot of the exact equation. To appreciate this, compare some of its numerically obtained solution curves with Poincaré plots of the exact equation (Fig. 6). Not only is the topology of fixed points and separatrices preserved but, indeed, so is much of their location (when giving proper notice to the rescaling factor ɛ, which is present in between these two systems). The modulation equations, however, have no restrictions on the region attainable in *X*–*Y* phase space, while the exact equations do. It is probably this restricted access of the original phase plane that makes a (formally) local solution perform so well in approximating its “global” dynamics. We will find that all of the behavior found near resonance in the original phase space is mimicked in the phase plane of the modulation equations (including, as we will see, the occurrence of chaos), which therefore forms a useful substitute.

### b. Hamiltonian description

In the next section, particular attention will be given to the inviscid (*C* = 0) limit of (11)–(12) since these equations then turn into Hamiltonian form:

with (for *θ* = 0) Hamiltonian

Here a constant has arbitrarily been included to simplify the expression. Multiple equilibria, in this case, exist for |*F*| < 8/3. Because of the symmetry of the modulation equations (11)–(12), **X** → −**X**, when *F* → −*F,* it is sufficient to consider *F* to be single signed only. From here on we take this to be the interval: −8/3 < *F* < 0. Equilibria are, of course, more readily identified from (11)–(12). With *θ* ≡ 0, they have *Y* = 0, while *X*(*X*^{2}/12 − 1) = *F*/2. With auxiliary variable *α* ∈ (−*π*/2, *π*/2), defined through *F* = −8 sin(*α*)/3, these steady states are in this case given by

and *k* ∈ {−1, 0, 1}, which implies *X*_{−1} < *X*_{0} < *X*_{1}.

### c. Homoclinic orbits

The motivation for studying the inviscid case (*C* = 0) is related to the fact that the modulation equations can then be solved explicitly. This offers the means to calculate Melnikov functions explicitly, once, at a later stage, weak (slow timescale) forcing and damping are added to the modulation equations. Zeros of the Melnikov functions determine the presence of invariant chaotic sets and suggest the presence of chaotic solutions. The conditions under which zeros of these functions appear will map out regions in parameter space where one may find “chaotic tides.”

In the inviscid case (*C* = 0) explicit solutions of (11)–(12) are obtained by introducing

which measures the radial “distance” to the “special” radius *R*_{s} = 23. Note that it obeys the equation

where we have used *θ* ≡ 0. Hence, by multiplying (11) by *F*/6 one obtains, upon integration,

where *K* is an integration constant related to the energy level. By using *Y* = ±(*R*^{2} − *X*^{2})^{1/2} in (19) and then eliminating *X* from (20), and *R*^{2} from (18), one obtains

where explicit expressions for *S*_{1,2,3}, in terms of *F* and *K,* are found by identification of both expressions (see appendix A). Plots of the graphic *P*_{4}(*S*) appearing in (21) are given in Fig. 7a. These correspond, for any given energy level, to the orbits in the *X*–*Y* phase plane, given in Fig. 7b. Of particular interest is the case that the central zeros of the quartic coalesce, and the corresponding orbits represent the homoclinic orbits *γ*_{in,out}. These homoclinic orbits correspond to the separatrices described in section 2b. Equation (21) can be integrated and its solutions can be expressed in terms of elliptic functions. However, we will have no need here for the general expressions and will only describe the homoclinic orbits coming from the saddle point (*X*_{1}, 0), which will be used in section 5a to explicitly calculate Melnikov functions.

The solution of (21) for the outer homoclinic orbit, *γ*_{out}, for which *S*_{2} < *S* < *S*_{1} and *K* = −4*S*_{1} − 3*S*^{2}_{1}, is obtained in appendix A as

where *S*^{′}_{2,3} = *S*_{2,3} − *S*_{1}. This shows that *S* → *S*_{1} as *T* → ±∞ and *S* = *S*_{2} for *T* = 0. Similarly an expression for the inner separatrix can be obtained by simply interchanging *S*_{2} and *S*_{3}, which has *S* = *S*_{3} for *T* = 0 and which also approaches the saddle point *S* → *S*_{1} as *T* → ±∞. With these expressions for *S*(*T*), the trajectories *γ*(*T*) = (*X*(*T*), *Y*(*T*)) of the inner and outer homoclinic orbits follow from (20) and *Y*(*T*) = ±(*R*^{2} − *X*^{2})^{1/2}; see Fig. 7b.

## 4. Forcing at two nearly resonant frequencies

The tidal forcing term *ζ*_{e}(*t*) is of course not a perfect cosine. Here we consider a somewhat more realistic double-frequency forcing:

The first component of this external tidal forcing might represent the semidiurnal lunar *M*_{2} tide (dimensionally, *ω*_{1} = 1.4056 × 10^{−4} s^{−1}); the second component the semidiurnal solar *S*_{2} tide (*ω*_{2} = 1.4544 × 10^{−4} s^{−1}). Of course this means that |*ω*_{1} − *ω*_{2}|/*ω*_{1} = 3.47 × 10^{−2} ≪ 1. As a consequence we cannot apply the ideas for the analysis of weakly, quasiperiodically forced oscillators, developed in, for instance, Beigie et al. (1991). However, we will see that the closeness of the two frequencies *ω*_{1} and *ω*_{2} can provide an additional periodic forcing term in the modulation equations, on the slow timescale. This periodic forcing term will break up the homoclinic orbits of the modulation equations (in the integrable limit) and thus might lead to (transversely) intersecting stable and unstable manifolds for the Poincaré map associated to this new forcing term.

Thus, taking into account an additional solar, and therefore almost resonant, tidal forcing term gives a mechanism for the construction of chaotic solutions (see also Yagasaki 1990, 1993).

In section 3 we defined the small parameter ɛ as a detuning. Analogously, we set here (now again employing nondimensional frequencies),

The above-described mechanism can be constructed if we assume that *ω*_{1} − *ω*_{2} = *O*(ɛ^{2}); that is, if

and ΔΩ = *O*(1). The forcing amplitude and the phase now vary on the slow timescale *T* = ɛ^{2}*t*:

where the combined amplitude and phase are defined by

The derivation of the modulation equations is now identical to the derivation given for the single-frequency case in section 3 (now one averages over 2*π*/*ω*_{1} instead of 2*π*/*ω*) and the modulation equations are again of the form (9)–(10), or, in Cartesian coordinates, (11)–(12), but now with slowly modulating forcing *F*(*T*) = *f*(*T*)/ɛ^{3} and phase *θ*(*T*), Eqs. (27) and (28) replacing their constant counterparts (see Nayfeh and Mook 1979, section 4.4). The right-hand sides of the modulation equations (11)–(12) now depend on the slow time explicitly, and, in principle, this offers the opportunity for the occurrence of aperiodic solutions.

When the two components are of equal magnitude, *f*_{2} = *f*_{1}, the phase is simply *θ* = (ΔΩ *T* + *θ*_{1} + *θ*_{2})/2, which redefines the effective frequency to be equal to the *average* frequency. The amplitude gradually changes sign over a long timescale *f* = 2*f*_{1} cos[(ΔΩ *T* + Δ*θ*)/2]. The block version of this, *f* = 4*f*_{1}(1/2 − Θ{cos[(ΔΩ *T* + Δ*θ*)/2]}, is piecewise integrable. Here Θ(*x*) = 0 (1) for *x* < 0 (>0) denotes the Heaviside function. It has a saddle point on the *X* axis, whose position alternates with respect to the center of the *X*–*Y* phase plane with subsequent phases of the forcing. It is thus reminiscent of Aref's (1984) blinking vortex with its ensuing chaos. This suggests similar results for a sinusoidal forcing, which is confirmed by numerical analysis (see section 4b), although this cannot be substantiated analytically.

When the second component is but a small perturbation with respect to the first one *f*_{2}/*f*_{1} ≡ *δ* ≪ 1, then *f* ≈ *f*_{1}[1 + *δ* cos(ΔΩ *T* + Δ*θ*)] and *θ* ≈ *θ*_{1} + *δ* sin(ΔΩ *T* + Δ*θ*). Note that this occurs quite naturally in the context of *M*_{2} and *S*_{2} tides. In contrast to the previous case, this limit is more amendable to further analysis and will be addressed now. Without loss of generality we set *θ*_{1} = 0, fixing the origin of the “fast” time *t,* and Δ*θ* = −*π*/2, fixing the origin of the “slow” time *T,* whence the amplitude *f* = *f*_{1}[1 + *δ* sin(ΔΩ *T*)] and phase *θ* = −*δ* cos(ΔΩ *T*), so that, with *f* = ɛ^{3}*F,* for small values of *δ,* apart from a factor 2, the forcing terms in (11) and (12) are approximated as *F* sin*θ* ≈ −*δF* cos(ΔΩ *T*) and *F* cos*θ* ≈ *F*[1 + *δ* sin(ΔΩ *T*)], respectively.

### a. Perturbative second forcing component and Melnikov analysis

With a weak perturbative second forcing component (*δ* ≪ 1) and weak damping (*C* → *δC*) the averaged Eqs. (11) and (12) become, after an appropriate shift in the origin of the slow and fast time,

with *H* as in (15). These equations can be written symbolically as

Here a dot denotes differentiation with respect to the slow time *T.* Note that it is assumed implicitly that 0 < ɛ ≪ *δ.* One has to include higher-order terms in ɛ and perform a much more detailed perturbation analysis when *δ* = *O*(ɛ) [see Guckenheimer and Holmes (1983); one has to be careful with applying Melnikovs method in the averaged system in that case; see also DKM for a detailed discussion]. System (31) has a hyperbolic periodic orbit, with period 2*π*/ΔΩ, that is *O*(*δ*) close to the saddle point *P*_{0} = (*X*_{1}, 0) of the unperturbed (*δ* = 0) problem (see Fig. 7b), or, equivalently, the stroboscopic or Poincaré map 𝒫^{T0} associated with (31) has a fixed point *P*^{T0}_{δ} of saddle type *O*(*δ*) close to *P*_{0} (see also Guckenheimer and Holmes 1983, chapter 4, for more details). Time *T*_{0} defines which section is taken to construct this map. The stable and unstable manifolds of *P*^{T0}_{δ} are also *O*(*δ*) close to those of *P*_{0} in the unperturbed system. Thus, the perturbation will in general break open both homoclinic loops of the unperturbed system.

Neglecting for the moment the additional forcing term [last term in Eq. (30)], the effect of the damping is to draw trajectories inward across separatrices. Multiplying (30) with ∇*H,* the rate of change of the Hamiltonian is then obtained, on the separatrix (where *K* = −3*S*^{2}_{1} − 4*S*_{1}), as *dH*/*dt* = −3*C*(*S* − *S*_{1})(3(*S* + *S*_{1}) + 4)/2. Since on *γ*_{out}, *S* ≥ *S*_{1}, it has *dH*/*dt* ≤ 0, while on *γ*_{in}, the inner separatrix, *dH*/*dt* ≥ 0, because *S*_{3} ≤ *S* ≤ *S*_{1} and 3(*S*_{3} + *S*_{1}) + 4 = 4 − 3[−8*S*_{1}(1 + *S*_{1})]^{1/2} ≥ 0, for −2/3 ≤ *S*_{1} ≤ 0. Since outside *γ*_{out} the Hamiltonian has relatively large (*K*) values (see Fig. 7b) and since its variation is strictly negative (*dH*/*dt* < 0), this implies that all states outside it are drawn inward. Since the Hamiltonian is again relatively large inside *γ*_{in}, any trajectories that reach up to that separatrix are further drawn inward, to the choked state. All other states end up in the banana-shaped area in between *γ*_{in} and *γ*_{out}, approaching the amplified state (with relatively low values of the Hamiltonian *K*). With the additional forcing term, multiplication of (30) with ∇*H* gives the general expression for the rate of change of the energy (i.e., that of the perturbed Hamiltonian, *dH*/*dt*). If we evaluate that at the separatrices, we specify *X*(*T*) and *Y*(*T*) to satisfy the parametric relation corresponding to the particular separatrix (section 3c). The net change of the Hamiltonian, following this trajectory over one period (which is infinite on the separatrices), then amounts to the integral given below. Its vanishing would imply that the net change in energy is zero, and hence forcing and damping are in equilibrium and the state does not drift, representing a new, dynamic equilibrium.

A geometric interpretation of the Melnikov method is that it measures the splitting distance between the stable and unstable manifolds of *P*^{T0}_{δ} that used to be part of a homoclinic loop to *P*_{0} in the unperturbed case. In this case there are two Melnikov functions: one that measures the splitting distance in the perturbed, outer loop, *M*_{out}(*T*_{0}), and one for the inner loop, *M*_{in}(*T*_{0}). In both situations this distance is for instance measured at the intersections of the stable/unstable manifolds with the *Y* = 0 axis; Fig. 7b. Following, once again Guckenheimer and Holmes (1983, chapter 4) we see that *M*_{out,in}(*T*_{0}) are, at leading order, given by

where *γ*_{0}(*T*) = (*X*_{0}(*T*), *Y*_{0}(*T*)) is either the outer or the inner homoclinic orbit of the unperturbed system (see section 3c); **g** = (*g*_{1}, *g*_{2}) and **h** = (*h*_{1}, *h*_{2}) are defined by (30) and (31).

A (nondegenerate) zero of, for example, *M*_{out}(*T*_{0}) corresponds to a transversal intersection of the “outer” stable and unstable manifolds of *P*^{T0}_{δ} for the map 𝒫^{T0}. It is a standard procedure to associate a Smale horseshoe map to this situation (Guckenheimer and Holmes 1983). Thus, it follows that system (31) has chaotic solutions when either *M*_{out}(*T*_{0}) or *M*_{in}(*T*_{0}) has zeros. Note, however, that the existence of chaotic solutions does not automatically imply the existence of a chaotic *attractor*; see section 4b.

In this paper we only consider the “classical way” to construct chaotic solutions, namely due to intersections of either the inner stable and unstable homoclinic orbits, or the outer ones. Other possibilities, like the intersection of the *outer* unstable manifold of *P*^{T0}_{δ} with its *inner* stable manifold, are discussed elsewhere (DKM).

For the outer homoclinic orbit, Eqs. (30), (31), and the definition of the unperturbed homoclinic orbit *γ*_{0}(*T*) = (*X*_{0}(*T*), *Y*_{0}(*T*)), the Melnikov function (32) is calculated in appendix B. It is obtained as

where auxiliary variables *ψ* and *k* are related to forcing *F* and frequency difference ΔΩ by *F* = −83 cos^{2}*ψ*/(1 + 2 cos^{2}*ψ*)^{3/2} and *k* = −ΔΩ(3 + tan^{2}*ψ*)/2 tan^{2}*ψ,* where *F* ∈ (−8/3, 0) for *ψ* ∈ (*π*/2, *π*). The derivation of the equivalent of *M*_{out}(*T*_{0}) for the *inner* homoclinic orbit, *M*_{in}(*T*_{0}), is completely similar (see appendix B). There is a critical value of the damping *C* = *C*_{out} = *C*_{out}(ΔΩ, *ψ*) = *C*_{out}(ΔΩ, *F*) for which the Melnikov function (33) vanishes, thus implying that there will be transverse homoclinic intersections for the *outer* orbit and chaotic solutions, when

depicted in parameter space in Fig. 8a. A similar critical surface *C*_{in} is obtained in appendix B for the inner orbit; see Fig. 8b.

Note that both *C*_{out} and *C*_{in} vanish when the two basic frequencies are equal (i.e., ΔΩ = 0). In this case both *M*_{out}(*T*_{0}) and *M*_{in}(*T*_{0}) are equal to *M*_{1} (for definition of *M*_{1,2,3} see appendix B), which does not depend on *T*_{0}, and thus, in this case there can be no (transversal) intersections of the stable and unstable manifolds for the Poincaré map 𝒫^{T0}. As a consequence there can also be no chaotic solutions. This has been observed earlier in the monofrequency case to which it reduces. The same is true for tidal frequencies that are not both close to the same resonance or, in other words, when |ΔΩ| is too large. This follows from (34) and its inner equivalent since the denominator sinh*πk* = sinh[*π*ΔΩ/(2*ν*)] dominates the numerator when |ΔΩ| ≫ 1 (the decay rate decreases to zero when *ψ* ↑ *π* (*F* → −8/3 and ΔΩ > 0; see Fig. 8a and below). Here an upward (downward) directed arrow implies the limit is approached from below (above). For ΔΩ < 0 the critical damping curve of the outer orbit is very small, *O*(0.01), and hence not visibly different from zero on the left-hand side of Fig. 8a.

There are two other physically relevant limits: *F* ↑ 0 (i.e., *ψ* ↓ *π*/2 and *ϕ* ↑ *π*/2) and *F* ↓ −8/3 (*ψ* ↑ *π* and *ϕ* ↓ 0). In the first case (”no forcing”) one expects once again the disappearance of (possibly) chaotic solutions. Indeed we find that

when *F* ↑ 0 [by (B9) and (B11): *M*_{2} and *M*_{3} vanish exponentially fast]. In a sense *M*_{out}(*T*_{0}) and *M*_{in}(*T*_{0}) measure the same splitting distance in this case since the outer and inner homoclinic orbits merge in this limit. The difference in sign is explained by the fact that the stable manifold of the inner orbit merges with the unstable manifold of the outer orbit (and vice versa). A similar geometrical point of view gives insight in the limit *F* ↓ −8/3. In that case the inner homoclinic orbit shrinks to a point, while the outer orbit reaches a well-defined limit. Thus, nothing “singular” will happen to the outer orbit in this limit (see Fig. 8a). However, for the inner orbit, all components *M*_{1,2,3} will approach 0 as *ϕ* ↓ 0. The behavior of *M*_{1} is quite degenerate: *M*_{1} = *O*(*ϕ*^{5}) for small *ϕ,* but both *M*_{2} and *M*_{3} will approach 0 with an exponential decay rate, for a *fixed* value of ΔΩ. Thus, we conclude that *C*_{in} ↓ 0 for *F* ↓ −8/3. However, the situation changes if we take the limit *ϕ* ↓ 0 differently, namely along any straight oblique line *through* the “origin” (ΔΩ, *ϕ*) = (0, 0): ΔΩ = *qϕ,* (*q* ≠ 0, ∞), as the exponentially fast decay is then absent and the behavior of *C*_{in} is dominated by the singularity in it due to the *ϕ*-dependent factor in *M*_{1}. This implies that the critical damping magnitude may increase indefinitely in this limit. This singularity also dominates the magnitude of *C*_{in} farther away from the origin, in the interior of the ΔΩ–*ϕ* plane, and explains the two lobes visible in Fig. 8b on either side of the line ΔΩ = 0. It may be the reason why chaos can become important, even in quite strongly damped, realistic circumstances.

### b. Numerical solutions of modulation and exact equations

Proving that the modulation equations possess a chaotic invariant set in the case when the Melnikov function vanishes is as far as we go analytically. As for the (closely related) Duffing equation, more work is needed to show that this set becomes attractive in some parameter range (Wiggins 1990, 612–613). In the following we will only give some numerical support that the (periodically forced) modulation equations possess such a chaotic (strange) attractor. This, in turn, suggests that the original, quasiperiodically forced Helmholtz oscillator should likewise have a strange attractor in the corresponding parameter regime; a result that we confirm by numerical integration below.

Numerical integration of the modulation equations, (11)–(12), with a fourth-order Runge–Kutta scheme, employing a double-frequency forcing (23), shows that for two nearby frequencies within the frequency range over which multiple equilibria exist, the solution curves appear to be chaotic (Fig. 9). This is suggested both in the (slow) time domain *T,* by the strait's current speed *Y* (Fig. 9a), as well as in the phase space of current velocity *Y* versus excess volume *X* (Fig. 9b). (Recall that the original **x** and averaged quantities **X** are related by **x** ∝ ɛ**X**, while the fast *t* and slow *T* timescales are related by *t* = ɛ^{−2}*T.*) Here we take forcing amplitudes of the two external tidal components to be equal. One of these components is supposed to be at the (linear) resonance frequency, 1 (the Helmholtz frequency), while the second is at a slightly higher frequency (1.01). The latter frequency is setting the small parameter ɛ = 0.1. Within the resonance band the theory requires the forcing amplitudes to be small, *O*(ɛ^{3}), so that we take them *f*_{1} = *f*_{2} = 10^{−3}. This leads, from (27), to *F*(*t*) = *f*/ɛ^{3} = 2 cos[(ΔΩ *T* + Δ*θ*)/2] in the modulation equations (12). In striking contrast to the monofrequency case, where the Hamiltonian chaos disappears with the introduction of even the smallest amount of damping (see M), the present calculations show that chaos survives the addition of damping. Commensurate with the theoretical requirements we need to take *c* ≤ *O*(ɛ^{2}) and, in fact, take *c* = ɛ^{3} here. This amounts to *δC* = 0.1. The erratic solution curves suggest that the modulation equations may possess a strange attractor, which is indeed revealed by making a Poincaré plot, by sampling at the modulation period (Fig. 10). Its shape is very similar to the “Japanese attractor” (see Ruelle 1980), which has been encountered for the forced and damped Duffing equation. Notice, in particular, the familiar Cantor-like division of the attractor's branches.

Chaos in the modulation equations suggests the occurrence of chaos in the original equations (5) under double-frequency forcing (23). In the following numerical experiments, the same values for forcing amplitudes and frequencies, as well as damping, were taken as were used in the previous numerical integration of the modulation equations. In order to compare the slow modulation of the numerically computed orbits directly to those obtained from the modulation equations, we sample the former at the tidal period. So, in Fig. 11 Poincaré plots of damped (*c* = 0.001) flow rate *y* versus time *t* (panel a) and versus excess volume *x* (panel b) are shown. Sampling is done once every average period 2*π*/*ω*, where *ω* = (*ω*_{1} + *ω*_{2})/2. Both figures suggest that chaos is also present in the exact equations. Indeed, this is again made manifest by making a Poincaré plot, but now by sampling at the long-period timescale, 2*π*/Δ*ω. *Figure 12 shows that it reveals the presence of a strange attractor, not unlike that encountered for the modulation equations; Fig. 10.

## 5. Forcing at two widely separate frequencies

Estimating the natural frequency of a number of basins occurring in nature leads to a wide range of periods varying from a few minutes to several hours, up to periods of the tides, with a predominancy in the range of several tens of minutes (Honda et al. 1908). For this reason it is useful to consider here the response of the nonlinear Helmholtz resonator to a forcing of the form

with 0 < *b* ≪ *a*: the first term represents the dominant, low-frequency tidal component and the second term a relatively weak perturbation that is almost in resonance with the Helmholtz frequency of the basin (*ω* ≈ 1). Thus, in this setting, ɛ^{2} ≪ 1 measures the ratio of tidal to Helmholtz frequency. Note that we have for simplicity neglected all “intermediate” Fourier components of the external tide: the main purpose of this section is to show that chaos can occur in basins with a high Helmholtz frequency. The mechanism is quite similar to the one described and studied in the previous sections. We will first discuss the effect of such a forcing in qualitative terms and continue this with a more detailed treatment.

Consider first the case in which the perturbation is absent altogether, *b* = 0. Then, when damping is relatively weak because the frequency of the tide ɛ^{2} ≪ 1, the elevation within the basin will be able to follow the external tide: *ζ* ≈ *Z.* Therefore, in the presence of a high-frequency perturbation (0 ≠ *b* ≪ *a*), on the fast timescale the dominant component of the tide, *Z,* acts as a slow, quasi-adiabatic change of the mean depth of the basin. However, the basin depth (among other parameters) determines the Helmholtz frequency; see (1). Therefore, the frequency of the perturbation, assumed to be fixed in absolute measure, is slowly varying when scaled with this Helmholtz frequency. This will provoke the quasi-stationary response to drift along the response curve (see Fig. 5). When this frequency drift, induced by the apparent modulation of water depth, is large enough—covering the frequency range over which multiple equilibria exist, this will lead to a sequence of “catastrophes”: rapid changes in the range of the (near) resonant, high-frequency response, which is subject to periodical collapse and expansion. There will thus be a hysteretic change in the amplitude of the basin elevation. This change may, however, be chaotic because the moments in time at which the amplitude jumps can become randomized. This is perhaps so because the high-frequency response, if not rapidly constrained to a particular phase by strong damping, will still be in a transient adapting phase, recovering from the previous catastrophe, when it undergoes the next one.

It is ironic that this explanation, based on the possibility of undergoing a hysteretic curve requires fairly large, but not too large, damping magnitudes: if damping gets too small, the frequency range over which multiple equilibria exist becomes broader than the range of the modulating, detuned frequency, and the response is “stuck” to a particular branch; if, on the other hand, damping is too strong, the response curve has no multiple equilibria at all, again lacking the ability to shift branches.

In practice, the conditions under which chaos appears may become weaker because not only will equilibria shift location under slow variation of the mean depth, but so, and probably more importantly, will the corresponding domains of attraction. This is perhaps the reason for the apparently chaotic behavior found in Fig. 13 under a mild modulation of the forcing frequency. This figure shows the result of a numerical integration of (2) in phase space with a forcing of type (35). Both results of direct integration (Fig. 13a), as well as a sampling of these on the long, tidal timescale (Fig. 13b) are shown. The direct integration shows that the dominant resonant response on the short Helmholtz timescale vaccilates on the long timescale, between a small- and large-scale response. These small- and large-scale responses correspond to the amplified and choked regimes in the amplitude response diagram of Fig. 5 to which these solutions would tend in the absence of the tide (*a* = 0). Under modulation with a tide of amplitude *a* = 3ɛ^{2} this provokes a modulation in apparent frequency *σ,* of approximately 1.5 dimensionless units, which is not big enough to span the entire frequency range over which multiple equilibria exist. Both this figure, as well as the tidally subsampled version in Fig. 13b (with period 2*π*/ɛ^{2}), however, testify about the irregular nature of the response, the latter figure in particular bearing some resemblance to the earlier Poincaré plots. In order to verify these qualitative ideas one may proceed with a more quantitative analysis. For the sake of brievity, this will however just be sketched here.

Again consider a basin in which the area increases linearly with the depth so that, with (4) and (35), the evolution of the excess volume *υ* is described by

[see (2)]. The homogeneous problem (*b* = *c* = 0) has a “quasi-stationary” equilibium at *υ*(ɛ^{2}*t*) = *Z*(ɛ^{2}*t*) + ½*Z*(ɛ^{2}*t*)^{2}, so it is natural to make the following shift in *υ*:

which yields

where *dZ*/*dt* = ɛ^{2}*dZ*/*dT* = ɛ^{2}*Z*′. As in section 3 we rescale *υ̃* = ɛ*V, **c* = ɛ^{2}*C,* and *b* = ɛ^{3}*B* so that we can start the derivation of a modulation equation for *V*:

[see (7)]. Thus, as predicted by the above qualitative arguments, the leading-order Helmholtz frequency 1/1 + *Z*(ɛ^{2}*t*) varies slowly and periodically in time ɛ^{2}*t* = *T.* Since we intend to proceed along the lines of the preceeding sections, we need to assume that the perturbation *b* cos(*ωt* + *θ*) is almost in resonance with the Helmholtz frequency of the basin. Hence we write (again) *ω* = 1 + ɛ^{2}*σ* and assume that *Z* = *a* cosɛ^{2}*t* = ɛ^{2}*A* cosɛ^{2}*t,* so that 1/1 + *Z* = 1 − ½ɛ^{2}*A* cosɛ^{2}*t* + *O*(ɛ^{4}). This last assumption is especially motivated by the mathematical approach, but may be quite realistic in the case of the tides in fjords too. Note that the scalings of *a* = ɛ^{2}*A* and *b* = ɛ^{3}*B* are consistent with the assumption *b* ≪ *a.*

As motivated by the above arguments we now introduce a new timescale *t̃* = *t*/1 + *Z* = *t*(1 − ½ɛ^{2}*A* cosɛ^{2}*t* + *O*(ɛ^{4})), so that (39) indeed reduces to (7) with *t* replaced by *t̃* and

at leading order on a timescale of *O*(1/ɛ^{2}). The derivation of the modulation equations is now completely similar to that of (9)–(10) in section 3 (as in section 4) and results in (9)–(10) with *σ* replaced by the periodically varying function *σα*(*T*). As in section 3 we set *θ* = 0 and scale *σ* to 1 so that the modulation equations are now given by

The presence of the periodically oscillating term can indeed be interpreted as a drift along the response curve of Fig. 5. This system can again be analyzed by the Melnikov approach of section 4a if we assume that *C, **α* = *O*(*δ*) with 0 < ɛ ≪ *δ* ≪ 1. Note that in this case the drift along the response curve is also only of order *δ.* The earlier described emergence of chaos, qualitatively attributed to hysteretic changes, which can only exist for *δ* = *O*(1), is apparently overrestrictive, as a simple “jitter” of the saddle point seems to be sufficient for the emergence of crossing of stable and unstable homoclinic orbits, and thus for the presence of a chaotic invariant set. For small *δ* the Melnikov analysis proceeds along the lines of section 4a. The only difference is that the perturbation term **h**(*X, **T*) in (31) has changed. It is a straightforward procedure to write down the equivalents of (B1) and redo the calculations of section 4b for the “new” inner and outer Melnikov functions, by which equivalents of (33) can be obtained. Note that (by construction) the structures of *M*_{out,in}(*T*_{0}) will be the same as in section 4a: *C* × (a constant term *M*_{1}) + (a *T*_{0} periodic term); *M*_{1} will even be exactly as in (33) and its inner equivalent. We do not present the details of the computations here. Our main observation is that one can thus again show the existence of a Smale horseshoe: it is the same mechanism as that of the quasiperiodically forced case of sections 3 and 4, which creates chaotic solutions to systems with two widely separated frequencies.

Several aspects of the numerical integration in Fig. 13 are still unrealistic. First, the response is dominated by the local resonance and thus the elevation shows a nearly periodic change on the *fast* timescale, instead of on the (long) tidal timescale. Second, the response is much too big, covering (almost) the entire basin depth. More work is therefore needed to bring the response closer to the ordering observed in nature, as for example, in Moldefjord (see Fig. 1), where Helmholtz response, *O*(0.1 m) ≪ tidal range, *O*(1 m) ≪ total depth, *O*(80 m), while currents on the fast Helmholtz timescale are comparable to tidal currents. It should be noticed that the correct ordering can partly be obtained quite simply by choosing model parameters such that all fixed points of the modulation equations are close to the origin. With this choice the homoclinic orbits, and therefore the actual state of the oscillator, are close to the origin of the phase space. Hence, corresponding elevations will be small compared to water depth. Notwithstanding the remaining issues that need clarification, it is believed that the mechanism described in this section may act as a “building block” for the explanation of the appearance of chaotic tides in basins with high Helmholtz frequency, particularly when the tide is extended with additional tidal components and harmonics.

From a practical point of view, any observed persistence of observed secondary oscillations might point at a non-meteorological origin. When their (average) period compares to estimates of the Helmholtz frequency, the model advanced here may be appropriate. Additional support may be obtained by plotting observed elevation (volume) against a strait current. Particularly, a stroboscopic sampling thereof (at the Helmholtz and tidal periods) may elucidate the structure of the underlying attractor. By comparison to predictions from the present model, an estimate of the most elusive factor, the damping, may perhaps be inferred.

## 6. Summary

Bays and estuaries may co-oscillate with the tides in the adjacent seas (Defant 1961). The extent to which they do generally depends on the geometry of the basin. Two competing mechanisms exist that mainly determine the final state of the coastal tides: the proximity of one of the more prolific tidal frequencies to basin resonances and the amount of friction. Depending on the balance between forcing and damping the coastal tide may be amplified or choked. This twofold nature of the response may actually occur for one and the same set of geometric and frictional parameters, as was shown in the case of the Helmholtz resonator for an almost-enclosed, relatively deep tidal basin with sloping bottom (see M). The slope in the bottom provokes a nonlinear restoring term (Green 1992) that, in turn, is responsible for the occurrence of multiple (stable) equilibria near the (linear) resonance frequency. When forced at the entrance by a singlefrequency tide (located within the resonance band), depending on the initial state, the tide within the basin may either have a small or large range that, after the decay of transients, is stationary.

Here we have shown that the presence of a second tidal frequency component, also within (or near to) the resonance band, may, in contrast, lead to the occurrence of a strange (chaotic) attractor, even with relatively strong damping. The basin tide may thus exhibit a perpetual change, manifested by unstable estimates of the tidal “constants” (amplitude and phase), when based on the analysis of a finite length time series, which may be of relevance to the observation by Gutiérrez et al. (1981) that some tidal components appear “unresolvable.” It may also be relevant to chaos in consecutive tidal maxima, observed in Venice Lagoon (Vittori 1992) or in coastal tides (Frison et al. 1999).

Surprisingly, there have been many observations of “irregular” secondary oscillations, accompanying the tide, which date back long before the notion of chaos appeared in the literature. Already in 1908, Honda et al. devoted an extensive study to ascertain the presence of both regular and irregular secondary oscillations in vertical elevation records in some 50, mostly Japanese, basins, usually based on observations covering just a few days. The nature of the basin response is typically, however, a superposition of the tide and a high-frequency oscillation, rendering the previous model, which predicts secondary oscillations on the tidal timescale, less applicable. Periods of these oscillations are observed to be in the range from several minutes up to several hours, depending on the size of the basin. Depending on the basin shape, this high-frequency oscillation can either be identified as the Helmholtz or quarter-wave resonance. The amplitude of the secondary oscillations in vertical elevation is, as the name suggests, observed to be an order of magnitude less than that of the (primary) tide. This may probably explain why, until fairly recently, the interest in secondary oscillations seems to have waned, and why the topic is held as a curiosity. A recent observation by Golmen et al. (1994), in a Norwegian fjord connected to the sea by a narrow strait, has, however, put the possible relevance of secondary oscillations into a new perspective by noticing that the associated secondary, irregular currents through the connecting strait reach a magnitude *comparable* to that of the tidal current. The reason for this is that, whereas the tide is present both outside and (with some delay) inside of the strait, the secondary elevation is present basically only within the fjord. The elevation *differences,* responsible for the associated currents, may therefore, at any time be of similar magnitude when comparing tidal and secondary oscillations. Hence, also the corresponding currents may be expected to be equally important. The same process may perhaps explain similar highly irregular current observations through lagoon entrances by Kjerfve and Knoppers (1991) and by Smith (1994), and the observed preferential enhancement of current harmonics (relative to that in the elevation field) by Seim and Sneed (1988). The relevance of this observation is that secondary oscillations may greatly alter the corresponding currents and therefore seriously modify the transfer of matter and dissolved substances (the “flushing” of the bay or fjord).

As remarked, these observations require a different kind of model in which tidal and resonance frequency—in our present model the Helmholtz frequency—are widely disparate. By adding a perturbative, but resonant, second forcing term to the primary tidal forcing it was shown that modulation equations are obtained that, to some extent, mimic those derived for the previous case with two tidal frequencies. In both cases, the modulation equations are driven at the long timescale. In the case of two nearby and resonant tidal frequencies this forcing stems from the beat (difference) frequency. In the case when the tidal frequency is much less than that of the resonant perturbative forcing term the tide itself is directly providing a modulation on the long timescale. The presence of a chaotic invariant set that could be ascertained analytically in these case, unfortunately does not guarantee that this set is also attractive. Therefore further support for the presence of a strange attractor was offered by (stroboscopic) Poincaré plots stemming from the numerical integration of the nonlinear Helmholtz resonator, when appropriately forced. The presence of a chaotic response may show up in a numerical simulation (in a basin with sloping shoreline) in much the same way as it reveals itself in nature: by irregularity in the elevation and, particularly, current fields (either within one tidal period, or in between subsequent tidal periods), by broadening of spectral peaks (with poorly resolved phases), and by showing a sensitive dependence on the initial tidal phase (neighboring orbits evolving to different states). The implication is that irregularity in the current and elevation fields is predictable as a phenomenon, but not its exact development. It requires further effort to quantify the resulting enhanced exchange between sea and basin.

Thus, despite the fact that, with the introduction of satellite-altimetry-derived tidal observations, tidal predictions have reached new unprecedented accuracies, it now seems that, at the same time, the tides may to some extent and in certain locations be intrinsically unpredictable, an unpredictability that reflects the nonlinear nature of the response of such locations.

## Acknowledgments

Discussions with Marko Wilpshaar, and comments and suggestions by Jef Zimmerman, Huib de Swart, and Ferdinand Verhulst are greatly appreciated.

## REFERENCES

**,**

**,**

**,**

_{4}tide.

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### Integrating Eq. (21)

In order for the right-hand of (21) to be real, the argument of the square root, the quartic

needs to be positive somewhere, which sets a bound on *K,* whose particular value is of no relevance here. Then, there can either be one or two intervals over which the quartic is positive, corresponding to the cases that the quartic has two or four real zeros, respectively; see the upper and lower curves in Fig. 7a. For fixed forcing *F,* these curves have differing values of *K.* In the former case this value is relatively low; hence the corresponding region in the *X*–*Y* parameter plane, in between *γ*_{in} and *γ*_{out}, is labeled with a minus sign; see Fig. 7b. In the latter case, it has a relatively high value of *K,* and there are two corresponding curves in the *X*–*Y* plane, one inside *γ*_{in} and one outside *γ*_{out}, which regions are therefore labeled with a plus sign. These regions are separated from each other by separatrices *γ*_{in,out} that emanate from the third fixed point: the saddle. These separatrices, or homoclinic orbits, are characterized by the fact that the two central zeros of the quartic coalesce, in which case (21) can be solved algebraically. Equating (A1) to the quartic −(*S* − *S*_{1})^{2}(*S* − *S*_{2})(*S* − *S*_{3})/4, and requiring equivalence of each of the coefficients of the fourth-degree polynomial, leads to four equations from which *S*_{1,2,3} and *K* can be determined. Because in the *X*–*Y* plane the central, double zero of the quartic corresponds with the saddle point (*X*_{1}; 0), whose position is already determined [see (16)–(17)], we find that one of these equations is redundant. The remaining equations yield the two other zeros,

and the energy level of the separatrix

where

(*Y*_{1} = 0) and *X*_{1} = 4 sin*ν*_{1}, where *ν*_{1} = (*α* + 2*π*)/3 and *α* = −sin^{−1}(3*F*/8). Over the range of *F* values considered, −8/3 < *F* < 0, one finds −2/3 < *S*_{1} < 0.

Replacing the quartic in (21) with its expression in terms of zeros requires the solution of

where we have assumed *S*_{3} < *S* < *S*_{2}. Now *S* can still be on either side of *S*_{1}, which itself lies in between the lower and upper limits *S*_{3} and *S*_{2}, respectively. By taking *S* either smaller or larger than *S*_{1} we obtain the inner or outer separatrix, respectively.

For definiteness assume *S*_{1} ≤ *S* ≤ *S*_{2} so that we are constructing the outer separatrix. Introducing

and similarly *S*^{′}_{2} = *S*_{2} − *S*_{1} > 0 and *S*^{′}_{3} = *S*_{3} − *S*_{1} < 0, then

Introducing *σ* = 1/*S*′ > 0 yields

Multiplying by *S*^{′}_{2} and defining *σ*′ = *S*^{′}_{2}*σ* > 1 and *η* = −*S*^{′}_{2}/*S*^{′}_{3} > 0 gives

Let *σ*′ = 1 + *P*^{2} and *η* = *Q*^{2} − 1 (with *Q*^{2} > 1), then

### APPENDIX B

#### Computing the Melnikov Function

With **g**, **h** from (30) and (31), the Melnikov function (32) is, for the outer homoclinic orbit, defined as

with a suffix 0 referring to spatial coordinates along an unperturbed homoclinic orbit. Since (18)–(20) apply in general, *S*_{0} = *R*^{2}_{0}/12 − 1, *FX*_{0}/6 = *S*^{2}_{0} − *K,* and *FY*_{0}/12 = *dS*_{0}/*dT.* Thus, along these orbits, the Melnikov function simplifies to

We will denote the three integrals in (B2) as *M*_{1}(*T*_{0}), *M*_{2}(*T*_{0}), and *M*_{3}(*T*_{0}), respectively. Note that *S*_{0} approaches the saddle *S*_{1} for *T* → ±∞ along the separatrices, so that the integrands vanish in these limits (recall that *K* = −3*S*^{2}_{1} − 4*S*_{1}). We now define

so that we can rewrite the *M*_{i}(*T*_{0}):

The expressions for *M*_{2,3}(*T*_{0}) were obtained by expanding cosΔΩ(*T* + *T*_{0}) and sinΔΩ(*T* + *T*_{0}) and by partial integration; all integrands involving sin(ΔΩ*T*_{0}) “average out” due to the odd/even symmetries of *S*_{0}(*T*) and its derivative. The *J*_{n} integrals can be computed explicitly by expressing them in an integral *I*(*a, **k*) that can be evaluated by a contour integral in the complex plane (see appendix C):

where *a* = cos*ϕ.* For *J*_{1}(ΔΩ) this is a rather straightforward procedure, since, by section 3c

where *A, **a,* and *ν* can be expressed, in terms of *S*^{′}_{2} and *S*^{′}_{3}, as

It follows that −1 < *a* < 0, and therefore we introduce *ψ* by

and hence, from (43),

so that

where

Note that the limit ΔΩ → 0, or *k* → 0, is well defined. Integral *J*_{2}(ΔΩ) can be obtained by taking the *a* derivative of *I*(*a, **k*):

Once again, the limit *k* → 0 is well defined. Compiling expressions, the Melnikov function for the outer homoclinic orbit, *M*_{out}(*T*_{0}) = *M*_{1} + *M*_{2} + *M*_{3}, is given by (33). The derivation of the equivalent of *M*_{out}(*T*_{0}) for the *inner* homoclinic orbit, *M*_{in}(*T*_{0}), is similar to the above analysis. The main difference is that the roles of indices 2 and 3 in the description of the homoclinic orbit [(B6) and section 3c] are exchanged and thus that *S*_{0} < *S*_{1}. This amounts to a replacement in all formulas of *ψ* → −*ϕ,* where *ϕ* ∈ (0, *π*/2), which is significant, because it now opens the possibility of a singularity in the surface of critical damping; Fig. 8b.

### APPENDIX C

#### Evaluation of an Integral

In order to evaluate Melnikov's distance function, integral *I*(*a, **k*), Eq. (B5), needs to be obtained for 0 < |*a*| < 1. Although this integral can be found in the literature (Prudnikov et al. 1986, p. 470) we here give a short sketch of its derivation by writing it as

for *p* = exp(*x*). Hence, by separating the denominator,

where, with *a* ≡ cos*ϕ, **p*_{±} = exp(±*iϕ*). Recasting the integration in terms of *x,* the integral reads

For definiteness assume *k* > 0. Then, by extending the integral over the real axis to a complex integration in the upper half-plane over a region enclosed also by a semicircle of infinite radius, this integral can be evaluated by means of Cauchy's theorem. It determines it as 2*πi* times the sum of the residues at the poles *x*_{±} = ±*iϕ* + (2*n* + 1)*πi,* for *n* = 0, 1, 2, · · · . Here, minus and plus signs refer to the poles appearing in the first and second integral, respectively. This yields

where *S*_{N} = Σ^{N}_{n=0 }*q*^{2n+1} with *q* ≡ exp(−*πk*), whence *S*_{∞} = *q*/(1 − *q*^{2}) = 1/(2 sinh*πk*), as 0 < *q* < 1. Therefore the expression in (B5) follows. Note that this expression holds also when *k* = 0.

## Footnotes

*Corresponding author address:* Dr. Leo Maas, Netherlands Institute for Sea Research, P.O. Box 59, 1790 AB Texel, Netherlands. Email: maas@nioz.nl

^{*}

Netherland Institute for Sea Research Publication Number 3288.