Is There a Tectonically Driven Supertidal Cycle?

Earth is 180 Myr into the current supercontinent cycle, and the next supercontinent is predicted to form in 250 Myr. The continuous changes in continental configuration can move the ocean between resonant states, and the semidiurnal tides are currently large compared to the past 252 Myr due to tidal resonance in the Atlantic. This leads to the hypothesis that there is a “supertidal” cycle linked to the supercontinent cycle. Here this is tested using new tectonic predictions for the next 250 Myr as bathymetry in a numerical tidal model. The simulations support the following hypothesis: a new tidal resonance will appear 150 Myr from now, followed by a decreasing tide as the supercontinent forms 100 Myr later. This affects the dissipation of tidal energy in the oceans, with consequences for the evolution of the Earth‐Moon system, ocean circulation and climate, and implications for the ocean's capacity of hosting and evolving life.


Introduction
The Earth moves through a cyclic dispersion and aggregation of supercontinents over a period of 400-500 Myr in what is known as the supercontinent cycle (Matthews et al., 2016;Nance et al., 1988;Rogers & Santosh, 2003). Pangea, the latest supercontinent, broke up around 180 Ma (Golonka, 2007(Golonka, , 1991, and it is predicted that a new supercontinent will form over the next 200-250 Myr (e.g., Duarte et al., 2018;Yoshida & santosh, 2011). The breakup of a supercontinent may lead to the formation of several internal oceans that will grow and eventually close. The lifecycle of each of these oceans is known as the Wilson cycle (Burke & Dewey, 1974;Wilson, 1966). Consequently, the completion of a supercontinent cycle through the formation of a supercontinent is generally preceded by the termination of several Wilson cycles (e.g., Burke, 2011). There is strong evidence that the tides are currently unusually large and that for most of the current supercontinent cycle, they have been less energetic than at present (Green et al., 2017;Kagan & Sundermann, 1996). The exception is the past 2 Myr, during which the continental configuration has led to a tidal resonance in the Atlantic (e.g., Green, 2010;Platzman, 1975). This (near-)resonant state has led to increased global tidal dissipation rates, which were further enhanced during glacial low stands in sea level (Arbic & Garrett, 2010;Egbert et al., 2004;Green, 2010;Green et al., 2017;Griffiths & Peltier, 2008;Wilmes & Green, 2014). An ocean basin can house resonant tides when the width of the basin, L is equal to a multiple of half wavelengths, = √ gHT (T is the tidal period, g is gravity, and H is water depth) of the tidal wave. Because today's resonant basin, the Atlantic, is currently opening, we would expect it to move away from the resonant state as it continues to widen, if we assume that the water depth and tidal period remain constant. Later on during the supercontinent cycle, we would expect either another basin to become resonant, that is, the Pacific if the Atlantic continues to open, or the Atlantic to become resonant again if it starts to close. This leads us to ask two questions: (i) when and where does this second resonance occur, if at all? and (ii) is there a supertidal cycle, that is, a cycle in the tidal amplitudes, associated with the supercontinent cycle? Here the questions will be answered by new simulations of the possible evolution of the tide over the next 250 Myr. This is done by implementing the tectonic scenario in Duarte et al. (2018) as the bathymetric boundary condition in the tidal model described by Green et al. (2017). Duarte et al. (2018) describe one possible scenario for the formation of the next supercontinent, Aurica, 250 Myr into the future. Aurica is predicted to be fairly circular and located in the present-day (PD) equatorial Pacific Ocean (see Duarte et al., 2018, and our Figure 1). It is formed by the closure of both the PD Atlantic and Pacific Oceans, which can only happen if a new ocean opens up. In the scenario, a bisection of Eurasia leads to the formation of a new ocean basin via intracontinental rifting. The motivation for the double-basin closure is that both the Atlantic and the Pacific oceanic lithospheres are already, in some regions, 180 Ma old (although the Pacific oceanic basin is much older; Boschman & van Hinsbergen, 2016;Golonka, 1991;Müller et al., 2008), and oceanic plates older than 200 Ma are rare in the geological record (Bradley, 2011). Consequently, it can be argued that both the Pacific and Atlantic must close to form the new supercontinent.
Changes in the tides, and the associated tidal dissipation rates, on geological timescales have had profound implications for the Earth system. Herold et al. (2012) and Green and Huber (2013) show that the changed location of abyssal tidal dissipation during the Eocene (55 Ma) can explain the reduced meridional temperature gradients seen in the proxy record for sea surface temperature that coupled climate models have struggled to reproduce (see Herold et al., 2012, for a summary). Furthermore, the reduced tidal dissipation during the Mesozoic and Cenozoic eras reported by Green et al. (2017) had implications for lunar recession rates and hence for interpreting cyclostratigraphy and long-term climate cycles (Waltham, 2015), for the evolution of the Earth-Moon system (Green et al., 2017), and for the evolution of life (Balbus, 2014).
The disposition of the continents on Earth over geological timescales consequently has a direct and major impact in the evolution of the Earth-Moon system, and tidal dissipation should be included in global oceanand climate models, especially over long timescales (Green & Huber, 2013). The overall aim of this paper is to evaluate if there is a supertidal cycle linked to the supercontinent cycle. We do this by expanding the work of Green et al. (2017) by adding tidal simulations 250 Myr into the future using the tectonic predictions in Duarte et al. (2018) as bathymetric boundary conditions. Consequently, this paper will increase our fundamental understanding of the Earth system, and it will, if the hypothesis is correct, lead to a first-order predictability of when large supertides may occur in Earth's history. To obtain this knowledge, we want to cover a full supercontinent cycle to see if there is a supertidal cycle. The logical thing to do is to expand the supercontinent cycle we are currently in into the future, because the first part of it has already been covered and shown to be tidally less energetic than PD (Green et al., 2017). In the next section we describe the tidal model and the bathymetric time slices used to obtain the results in section 3. Section 4 closes the paper with a discussion, conclusions, and an outlook into further work.

Tides
We use OTIS-the Oregon State University Tidal Inversion Software-to simulate the evolution of the future tides. OTIS is a portable, dedicated, numerical shallow-water tidal model, which has been used extensively for both global and regional modeling of past, present, and future ocean tides (e.g., Egbert et al., 2004;Green, 2010;Green & Huber, 2013;Green et al., 2017;Pelling & Green, 2013;Wilmes & Green, 2014). It is highly accurate both in the open ocean and in coastal regions (Stammer et al., 2014), and it is computationally efficient. The model solves the linearized shallow-water equations (e.g., Hendershott, 1977): Here U is the depth-integrated volume transport (i.e., tidal current velocity u times water depth H), f is the Coriolis vector, g denotes the gravitational constant, is the tidal elevation, SAL denotes the tidal elevation due to self-attraction and loading (SAL), and EQ is the equilibrium tidal elevation. For simplicity we used a constant SAL correction with = 0.1 (Egbert et al., 2004). F represents energy losses due to bed friction and tidal conversion. The former is represented by the standard quadratic law: where C d = 0.003 is a drag coefficient, and u is the total velocity vector for all the tidal constituents. The conversion, F w = C|U|, includes a conversion coefficient C, which is here defined as (Green & Huber, 2013;Zaron & Egbert, 2006) C(x, y) = (∇H) 2 N bN 8 2 (4) Here = 50 is a scaling factor, N b is the buoyancy frequency at the seabed,N is the vertical average of the buoyancy frequency, and is the frequency of the tidal constituent under evaluation. The buoyancy frequency is given by N = N 0 exp(−z∕1, 300), where N 0 = 5.24 × 10 −3 s −1 and based on a least squares fit to PD climatology values (Zaron & Egbert, 2006). The future stratification is obviously unknown, and to estimate potential effects of altered stratificaiton, we did a set of sensitivity simulations in which C was doubled or halved. As in other tidal simulations this had a relatively minor effect on the global tides, and we will not discuss these results further (see, e.g., Egbert et al., 2004;Green & Huber, 2013).

10.1002/2017GL076695
The model solves equations 1 and (2) using the astronomic tide-generating force as the only forcing (represented by EQ in equation (1). An initial spin-up from rest of over 7 days is followed by a further 5 days of simulation time, on which harmonic analysis is performed to obtain the tidal elevations and transports. Here we focus on the M 2 and K 1 constituents only.

Bathymetry Data 2.2.1. PD Bathymetries
The PD bathymetry is the same as in Green et al. (2017): see our Figure 1, top left panel. To avoid open boundaries, the equilibrium tide was used as forcing at 88 ∘ when appropriate. Tests with a vertical wall at the poles (not shown) did not change the results. All simulations were done with 1/4 ∘ horizontal resolution. The PD control simulation was compared to the elevations in the TPXO8 database, and the root-mean-square errors (RMSEs) were computed from the difference between modeled and observed elevations. TPXO8 is an inverse tidal solution for both elevation and velocity based on satellite altimetry and the shallow-water equations and is commonly taken as the truth for tidal elevations (see Egbert & Erofeeva, 2002, and http://volkov.oce.orst.edu/tides/tpxo8_atlas.html for details).
To evaluate the sensitivity of our solutions to the lack of detail in the future bathymetries, we constructed a simplified PD bathymetry having the same (lack of ) detail as the future bathymetries (see Figure 1, top right panel). This case is denoted PD reduced in the following, and it is the simulation we use as a benchmark for the evolution of the tide. In PD reduced any water currently shallower than 200 m was set to 200 m. PD oceanic ridges were smoothed out and set to have a peak depth of 2,500 m and a total width of 5 ∘ over which the ridge approaches the depth of the deep ocean linearly, whereas subduction zones were set to be 1 ∘ wide and 6,000 m deep with a triangular cross section. The remaining ocean was set to a depth computed to conserve the ocean's PD total volume. The same values were used in the construction of the future bathymetries shown in Figure 1.

Future Bathymetries
We used GPlates for the kinematic tectonic modeling of the future scenario (see Duarte et al., 2018;Qin et al., 2012, and https://www.gplates.org/ for a description). The continental polygons provided in the GPlates data repository were used as the starting point for the PD ocean, as in Matthews et al. (2016). The drift paths of the continental plates were constrained for the first 25 Myr by the drift velocities in Schellart et al. (2007). For the remaining 225 Myr, we used the PD globally averaged plate velocity of 5.6 cm year −1 (see Duarte et al., 2018, for a summary) but applied deviations from the average based on the observations in Zahirovic et al. (2016). The plate and land boundaries from the model were output as digital gray scale images, which were used to build the model bathymetries based on the details given for the PD sensitivity bathymetry.

Dissipation Computations
The computation of tidal dissipation rates, D, was done following Egbert and Ray (2001) and thus given by Here W is the work done by the tide-generating force and P is the energy flux given by where the angular brackets mark time averages over a tidal period.

PD Sensitivity
The PD control simulation (Figure 2, top left) has an RMSE of 11 cm when compared to TPXO8; the same computation for the reduced M2 tide gives 23 cm. The K1 RMSE are 2 cm for PD and 10 cm for PD reduced, respectively. As discussed above, we did a series of sensitivity simulations for both bathymetries in which the tidal conversion coefficient was changed within a factor of 2, and the RMSE and dissipation rates did not change significantly (not shown). Green and Huber (2013) and Green et al. (2017) did an extensive series of sensitivity simulations and came to the same conclusion. Consequently, we have confidence in the robustness of our results, and we have a well-constrained error bound on the simulations. The PD sensitivity simulation reveals a less energetic global tide (Figure 2, top right), with reduced M2 tidal amplitudes in the Atlantic and the emergence of fairly large M2 tides along the Siberian shelf and around Antarctica. The new tides along the northern coast of Eurasia are due to the sub-Arctic seas being deeper, allowing the tide to propagate into the Arctic Basin. The large PD Atlantic tides are reduced because of the water-world-like ocean and reduced shelf sea area, leading to a more equilibrium-like tide (see Egbert et al., 2004, for a discussion). The weaker M2 tide in the PD reduced scenario means that we are potentially underestimating the M2 amplitudes for all future scenarios. For K1 we see a different pattern in Figure 3 (top row): our synthetic bathymetry appears to produce a larger K1 amplitude than the PD bathymetry. This is likely because the changed water depth allows the K1 tide to be nearer resonance in some areas, such as around Greenland and Indonesia. It is thus possible that we are overestimating K1 in the future scenarios. Note that this is a reversed response to that in Green et al. (2017), where their simplified bathymetry gives an enhanced M2 tide. The bathymetries in Green et al. (2017), however, have more topographic detail, especially in shallow water, than the ones used here. For clarity, we will describe our globally averaged or integrated metrics in relative terms by normalizing by the respective values from the simplified PD bathymetry.

Tidal Amplitudes
The global M2 tidal amplitude increases slightly over the next 50 Myr (refer to Figures 2 and 4a for the following discussion) due to an enhanced tide in the North Atlantic and Pacific at 25 Myr, followed by a very large Pacific tide at 50 Myr. This is because the equatorial Pacific becomes half-wavelength resonant at these ages. At 100 Myr, there are large tides in the newly formed Pan-Asian Ocean (in the Asian rift) and in the Indian Ocean. This signal persists to 150 Myr when the Atlantic comes back into resonance to form the next tidal maximum. After 150 Myr, there is a decline of the global amplitudes as the new supercontinent starts to come together and the resonant properties of the basins are lost. When Aurica has formed fully at 250 Myr, we only see large tides locally, in embayments with a geometry allowing for local resonances. K1 follows a different pattern to M2, with a global tidal maximum when M2 hits a minimum at 100 Myr (the maximum K1 amplitude is then about 5 m). The average K1 amplitude remains relatively constant between 150 and 200 Myr, before a very sharp decline as the next continent forms (Figures 3 and 4a). It appears that K1 does not have two resonances in this tectonic scenario, whereas M2 does, because it becomes resonant again when the Atlantic closes, as well as in what is left of the Pacific at 150 Myr ( Figure 2). This makes sense from a basin size perspective: because the Atlantic continues to open for a while before closing again, K1 will never have an opportunity to become resonant in the Atlantic, whereas M2 will be. Because of the changing size of the Pacific, it will be resonant for the K1 tide at 100 Myr only.

Dissipation and Earth-Moon Evolution
Overall, the global M2 dissipation rates for the remainder of the supercontinent cycle is 84% of the present values, or 2.2 TW (Figure 4b). This expands the results in Green et al. (2017) 250 Myr into the future and strongly suggests that Earth is presently in an M2 tidal maximum. It also suggests that the maximum has a width of 50 Myr or less and that there will be another M2 maximum occurring during the cycle around 150 Myr from now, i.e., 100 Myr before the formation of the next supercontinent. K1, in contrast, will be resonant only once in the current cycle, at 100 Myr. This is in agreement with results for the late Silurian (430 Ma), which showed more energetic tides than during the Early Devonian (400 Ma; H. Byrne, personal communication, 2017;Balbus, 2014). Pangea, the previous supercontinent, formed around 330 Myr ago and started breaking up some 180 Myr ago. It thus seems plausible that Earth's oceans go through tidal maxima some 150-200 Myr after supercontinental break up (i.e., at present) and around 100 Myr before a supercontinent forms (i.e., during the Silurian, before Pangea, and 150 Myr into the future for Aurica).
Following the lunar recession model in Waltham (2015) and summarized in Green et al. (2017), the recession rate, a∕ t, can be written as a t = fa −5.5 (8)  (8) and (9). The distance is normalized by the PD distance, a 0 .
where f is the tidal factor given by Here m ′ = mM∕(m + M) is the reduced mass of the Moon (M = 5.972 × 10 24 kg and m = 7.348 × 10 22 kg are the masses of the Earth and the Moon, respectively), and Ω = 7.2923 −5 s −1 ( = 2.6616 × 10 −6 ) s −1 is the rotation rate of the Earth (Moon). Using the dissipation rates in Figure 4b, interpolated to every 1 Myr using linear interpolation to produce a smoother curve, we obtain the result in Figure 4c. We have also, for comparison, computed the lunar distance assuming a continuous PD dissipation rate (dashed). These results further highlight the conclusions in Green et al. (2017): appropriate tidal dissipation rates should be used in investigations involving lunar recession rates or distances, especially over long periods of time. Consequently, the PD recession rate is anomalously high because of the current tidal resonance in the Atlantic and PD tides are a poor proxy for past or future tides over large parts of the supercontinent cycle.

Discussion
Our results support previous ideas that the tides are at their lowest when the Earth is in the supercontinent configuration. The dissipation is then less than 40% of the PD value in our simulations. The tenure of a supercontinent varies, but both Pangaea and Rodinia, the two most recent supercontinents, maintained their formation for over 100 Myr (Rogers & Santosh, 2003). This means that dissipation rates could remain at this very low level for long periods of time-much longer than the timescale of its resonant peaks, which here are less than 50 Myr (see below for a tighter constraint).
This project aimed to evaluate if there is a supertidal cycle. The results strongly suggest that the answer is yes: there is a repeated gradual change between states of high and low tidal dissipation levels over the period of Aurica forming. However, there is more than one supertidal cycle within the supercontinent cycle. Combined with the results in Green et al. (2017), who goes back to Pangea 252 Myr ago, we suggest that the oceans will go through two M2 supertidal cycles and at least one K1 cycle during the current supercontinent cycle. Consequently, the global tides are weak for long periods of time and then pass through several quite narrow (on geological timescales) resonances. This is because there are several Wilson cycles involved in one supercontinent cycle, and as the basins open and close there can be several supertidal cycles associated with the Wilson cycles. This also means that the supertidal cycle is not necessarily in phase with the supercontinent cycle. The mechanism behind the supertidal cycle is tidal resonance, which is set up by the continental configurations: peak resonance occurs when the continental configuration results in an ocean basin of a length that is an exact multiple of half wavelengths of the M2 tidal wave. Theoretically, one would therefore be able to predict when each basin may be resonant, without being able to provide any details of the actual magnitude.
To lowest order, one can assume that the tide will be large when the natural frequency of a basin is within, say, 20% of the tidal period (see Figure 11 in Egbert et al., 2004, for a theoretical estimate). For the present, this would give a period window of about 3 hr in which the basin is close enough to resonance to support a large tide. If the ocean is 4,000 m deep and we are looking at a half-wavelength resonance, we get a range of the width of the basin in which it is resonant of about 1,100 km. With a continental drift rate of 6 cm year −1 , the width of the resonant peak would then be some 18 Myr, implying that Earth is currently in the beginning of the tidal maximum. There is further support for this in our results here and in Green et al. (2017). They show that the tides were weak 2 Myr before present, and the 25 Myr time slice in the present paper still shows a rather large tide. This is an interesting idea worth pursuing in a future paper, which would look into the time span of the resonances in more detail by simulating more time slices between now and 25 Myr.
Ocean basin closure is a result of consumption of oceanic plate at subduction zones within the basin and sea floor spreading in a neighboring basin. This means that there are long periods where the direction of closure is mostly fixed, with two continental plates being pulled or pushed together. The observation that there are multiple peaks in tidal dissipation makes sense in this context. There will be multiple modes of resonance for each ocean basin as it reaches the dimensions that are resonant for smaller or larger multiples of the tidal wavelength. The implication of these observations is that the length of a supertidal cycle is directly related to the length of the supercontinent cycle. Consequently, the period of the supertidal cycle is set by how quickly the continental configurations move from one resonant mode to another. There are of course other factors that contribute to the total tidal dissipation, such as sea level changes and variations in the extent of continental shelves, but through this study we have a clear indication that changing the position of the continents alone is enough to elicit significant changes in the energy of the tidal system. However, if an ocean basin is close to resonance it is much more sensitive to relative sea level changes and/or continental shelf configurations than when its not in a near-resonant state. This was the case for the Last Glacial Maximum (21-18 kyr) where Green et al. (2017) find the largest M2 amplitude in their 252 Myr time series. This exceptionally large tide is explained by a low stand in sea level, exposing the dissipative shelf seas (Egbert et al., 2004;Wilmes & Green, 2014).