On the Tidal History and Future of the Earth–Moon Orbital System

Earth’s rotation rate and the evolution of the Earth–Moon system have been controlled by tidal dissipation in Earth’s ocean. Attempts to model the tidal history have shown incomplete compatibility with observations and unclear isolation of the most important controlling factors. Here it is shown that a relatively simple model with no explicit description of the continents (their effects are instead parameterized) can accurately reproduce the available observations (lunar distance and month, Earth’s rotation and deceleration rates) describing the tidal evolution over the past 2.5 billion yr and also evade the paradox of an early Moon falling within the Roche limit. Notably, the model reproduces an observed dissipation peak 400 million yr ago. The model involves fitting two input parameters (the effective ocean depth h and the nondimensional dissipation timescale T˜ , the latter of which can be related to the more typically used quality factor Q). The best-fit values (h = 2.3 km and T˜=40 , corresponding to Q ≈ 20–28) derived empirically from the model correspond well with the plausibly expected values (h ≈ 2.5 km and Q ≈ 20) derived using independent reasoning and observations. The model shows very clearly that the tidal evolution has been primarily controlled by the rate of resonance in the ocean response, not the evolving amplitude of the tidal forces or forcing frequency. More closely, it also shows that the effect of Earth’s rotation on the Lamb parameter (or nondimensional wave speed) has driven the evolution, and that changes in frequency or other parameters have been much less important.


Introduction
Earth's rotation rate and the orbital distance of the Moon evolve due to the frictional drag of ocean tides. As the Earth rotates with respect to the Moon, an idealized ocean response would have a tidal bulge propagating so as to stay directly below the Moon. Because the "tidal" part of the gravitational force is a spatially differential component, there would also be an antipodal second bulge deforming the global ocean to look something like an American football with its major axis aligned with the Earth-Moon axis. The effect of frictional drag on this propagating bulge would be to give the bulge a time lag. The lagged bulge is no longer aligned with the Earth-Moon axis, and the result is that the Moon is pulled forward in its orbit. A similar idealized tidal response is expected due to the Earth's rotation with respect to the tidal pull of the Sun.
This idealized tidal response is simplistic, however, for several reasons. First, it should be expected that there are additional lags in the response if the fluid wave speeds are not fast enough to keep up with the propagation of the tidal forces. Second, the Earth's rotation creates a tendency to scatter the response toward smaller spatial scales. Third, the existence of continents creates a spatially inhomogeneous response and complicates the propagation pathways. Fourth, resonant amplification of the tidal response can occur, and these resonant modes depend sensitively on the ocean and orbital parameters.
Given the expected complexity of the tidal response and its dependence on ocean parameters that are uncertain in the past, it could seem doubtful that there can be much accuracy in modeling the tidal history of the Earth-Moon system. The purpose of this paper is to show that a relatively simple model of the Earth's tidal response is in fact capable of explaining the collection of important observations regarding the Earth-Moon history.
There is a long history of studies aimed at understanding the evolution of the Earth-Moon system and the role of ocean tides (MacDonald 1964;Gerstenkorn 1967aGerstenkorn , 1967bMunk 1968;Jeffreys 1973;Lambeck 1975Lambeck , 1977Darwin 1881;Hansen 1982;Webb 1982;Williams 1993Williams , 2000Kagan & Maslova 1994;Kagan 1997;Bills & Ray 1999;Egbert et al. 2004;Gao & Xiao 2008;Coughenour et al. 2009Coughenour et al. , 2013Na & Lee 2014;Williams et al. 2014;Green et al. 2017;Van Putten 2017;Meyers & Malinverno 2018;Blackledge et al. 2020;Motoyama et al. 2020). As a starting point on this important topic, the reader is referred to MacDonald (1964) and Munk (1968) for historical background, the concise review in Bills & Ray (1999), and the more recent papers cited above for updates. Notably, Motoyama et al. (2020) described a methodological approach similar to the one used here. A recent review of the observational data is described in Coughenour et al. (2009).
The data that will be used to evaluate the model performance in this study are identical to the data used in Bills & Ray (1999) but extended to include the more recent results. The observed data (shown in Figure 1) comprise estimates of the present and past values for the Earth's length of day (LOD), the length of the month (LOM), the orbital separation of the Moon (d 0 ), and the tidal dissipation rate. (Note that LOM and d 0 are related by Kepler's third law.) In some cases, the time derivatives of these data are calculated and available for model evaluation as well. corresponding to Q ≈ 20-28). The model history curves reproduce well the observed data points (see text for description) for lunar separation d 0 (panel (B)) and length of month (LOM; panel (F)). The resonance peak in tidal dissipation (indicated by the spike in ∂ t LOD about −400 Ma shown in panel (D)) also agrees. The model is in reasonable agreement with the data points for length of day (LOD; panel (A)) and the current dissipation rate shown by the diamonds (panels (C) and (E)).
In Figure 1, the red circles correspond to data in Denis et al. (2011), the black stars to data in Meyers & Malinverno (2018), the magenta squares to data in Bills & Ray (1999;after Sonett et al. 1996), the green asterisks to Coughenour et al. (2009), and the black square to data in Walker & Zahnle (1986); the blue diamond represents the current dissipation rate taken as the sum of the dissipations due to the M 2 and S 2 tidal components as listed in Table 1 of Egbert et al. (2004).
The current rate at which d 0 is increasing (about 3.8 cm yr −1 ) is obtained from several precise methods, one involving measurement using laser ranging off reflectors left on the Moon (Williams et al. 2014;Williams & Boggs 2016). This rate must, however, be high in comparison with the average over Earth's history because a simple backward extrapolation of this rate to early times leads to several strong inconsistencies, not least of which is a catastrophically close lunar orbit (or Moon much younger than expected). In addition to fitting the observational data points shown in Figure 1, a constraint for the model is that the orbital history should not show a d 0 less than the Roche limit (taken here to be 3 Earth radii).
There is also evidence from a variety of marine sedimentary and fossil assemblages that while the average rate of increase of day length over the last 640 million yr appears to be roughly similar to the present 2 ms century −1 (5.6 hr Gyr −1 ), the average rate over the 2 billion yr preceding this was about five times less (e.g., Denis et al. 2002). This then implies a "paradox" in that the tidal dissipation rate was lower during a time when the Moon was presumably closer and the associated tide-raising forces stronger. The resolution to this paradox is that the dissipation of rotational energy is presently large because the ocean's response to the tidal forces has entered a more resonant state and therefore draws more energy from the Earth's rotation than has been typical over at least the past 2 billion yr. Looked at more closely (e.g., Figure 1, panel (d)), time differencing of the sparse observations suggests that the current dissipation rate in fact passed through a peak about 400 million yr ago and has been declining since. A model of the tidal history must then show dissipation rates that are higher today than over most of the past and yet lower than during a peak in dissipation some 400 million yr ago. This abrupt and nonmonotonic behavior in the history of tidal dissipation appears to implicate changes in the level of resonance in the tidal response.
The reason for this enhanced resonance, however, is not clear. One early explanation proposed that this resonance is due to continental breakup and the spectral shift in normal modes associated with the smaller basin geometries (e.g., Webb 1982;Kagan & Maslova 1994), and this role of "continentality" in controlling the resonance has been extended to exoplanet ocean tides (Blackledge et al. 2020). The effects of the continents are surely important because of the barriers on propagation. The continental shelves are also important because increased flow speeds in the shallow water leads to higher dissipation rates and because of the potential back effects of the shelf tides on tides in the deep ocean (Arbic et al. 2009). The variation in this continentality effect on the tidal response is difficult to judge for the distant past, and perhaps most importantly, the methods previously used have provided only a very limited sampling of the full set of plausible ocean tidal solutions. The ocean tidal solutions that have been studied have been provided by an analytical approach allowing a variation in the placement of a hemispherical continent (Hansen 1982;Webb 1982), numerical integrations with continental configurations selected from the recent geological past (Brosche & Sunderman 1978), stochastic studies under idealized conditions (Kagan & Maslova 1994), or a spherical-harmonic solution routine as in Motoyama et al. (2020). These studies are informative but undersample the set of possible scenarios and are insufficient for providing or characterizing a description of the long history of the ocean's tidal response. The primary support for this hypothesis that continentality is responsible for pushing the ocean into a resonant state is that normal modes calculated for the realistic ocean indeed show a frequency overlap with that of the semidiurnal forcing (Platzman 1983;Sanchez 2008). However, these basin modes do not have the degree-two form required to couple with the tidal forcing; therefore, no separate role of continentality in the level of resonant forcing is established, even for the present ocean. When considering tides over geological history, the present resonant state is explained only as a relatively brief accident of continental drift, rather than a robustly predictable result of feedback (such as in the next hypothesis). Finally, it is important to consider that the additional basin modes provided by continentality are not directly forced by the astronomical tidal forces, so it is not immediately obvious how their existence extracts more work from the tidal forces. That is, resonant forcing requires not only a match in temporal frequency between force and normal mode but also a match in the spatial form. Laplace developed the spherical-harmonic expansion for compactly representing the projection of a gravitational potential field over a sphere, so it is not surprising that only the lowest, degree-two components need to be considered in the tidal forcing. When one considers the total work performed on the ocean (equating this to total dissipation), it is clear that only the similar degree-two component of the ocean's response couples with the forcing to perform work. Hence, the astronomical forces must first excite these global-scale, degree-two tides. Once excited, these large-scale tides can scatter energy into basin modes (for example, the resonant, large-amplitude tides in the Bay of Fundy are excited not by the Moon or Sun but by the tides in the Atlantic basin). Hence, while the additional basin modes provided by continentality are expected to indirectly influence the amount of total work and dissipation, this expectation does not follow from a simple frequency match between basin normal modes and forcing. Although the Bay of Fundy has a normal-mode frequency match with the lunar semidiurnal forcing, it would not be resonantly excited and add to dissipation if there was not a significant tidal excitation in the North Atlantic.
An alternative or additional hypothesis for the recent enhanced resonant ocean tidal response is aligned with statements made in earlier work (e.g., Hansen 1982) and strongly supported by the results in this study. In this case, it is not the continents but rather the rotation rate that is most important. It is the Earth's decreasing rotation rate that has had the largest systematic effect of pushing the ocean tidal response toward enhanced resonance. The continents may still have primary control over other factors, such as the parameter Q ("quality factor"), but the systematic effects of the continents on the ocean's normal modes are either secondary or average out over geological time, as suggested in Hansen (1982). It is sometimes misconstrued that the Q factor controls the amount of dissipation. In fact, Q is simply a dimensionless timescale describing the attenuation of the tidal energy, the level of which can depend far more sensitively on shifts of parameters affecting the normal modes and thereby resonance.
The question of the relative importance of continentality versus spin/orbital parameters in systematically controlling the evolution of tidal dissipation over long timescales is important to resolve because if continentality dominates, then there may be little expectation for reliably calculating the full Earth-Moon tidal history. Similarly, if continentality dominates, then conjectures on the ocean tides of exoplanets and the orbital evolution (e.g., Blackledge et al. 2020) will be difficult to evaluate. Interestingly, Motoyama et al. (2020) showed tidal modeling for the first half of Earth's history (assuming this precedes continent formation) that may have a more reliable basis than what is possible for the second half. The removal of continents from the model in Motoyama et al. (2020) is justified by applying the model only to the early time that may predate continental formation. The model in this study will remove any explicit description of continents by attempting to parameterize their effects, with the justification that when considering the longest timescales, continentality effects might be effectively averaged and represented by generic model parameters/terms fitted to observations. Aside from treating the full history, an advantage of this approach is that rather than forward-integrating in time from uncertain initial conditions, as in Motoyama et al. (2020), backward integration from the present state is possible. In this paper, it will be shown that a relatively simple model of ocean tidal dissipation including no explicit representation of continental history is capable of explaining the observations.

Methodology and Approach
In the subsections below, the method and model assumptions for calculating the tidal response are discussed, as is the description of the tidal forcing and the method for integrating the Earth/Moon/Sun tidal system backward in time from the present. Section 2.2 describes the approach in the modeling, most notably that the goal is to provide a model with the most reduced set of fitted parameters and model process terms. Models with more degrees of freedom can fit the data better but risk losing explanatory power if overfitting has occurred.
In the descriptions below, a few standard terms from geophysical fluid dynamics are used. First, tidal flows, or tides, are the fluid response to tidal forces. The response must be carried by fluid wave modes, the two most important of which in this application are rotational-gravity and Rossby waves. The rotational-gravity waves refer to gravity waves modified by rotation (here "gravity" refers to the restoring force and not the field in undulation). For wavelengths much longer than the fluid thickness, gravity waves become nondispersive and include as familiar examples swell and tsunamis. These (long) gravity waves necessarily become modified by the Coriolis force (rotation) because the timescales considered (the tidal periods) are comparable to the rotation period. Rossby waves (also called "Rossby-Haurwitz" or "planetary waves") have as a restoring force the variation of the Coriolis parameter with latitude. Long Rossby waves, as in this application, are also nondispersive and propagate only in the retrograde direction. Familiar examples of Rossby waves can be seen in the propagating patterns of low-/high-pressure cells on a weather map. In the early literature, the rotational-gravity and Rossby waves can also be associated with the "class I" and "class II" oscillations, respectively.

Calculation of the Fluid Tidal Response
To calculate the tidal response, the Tidal Response Of Planetary Fluids (TROPF) software package is used. The TROPF package, including a comprehensive manual describing the formulation and code, is available for download (Tyler 2019a). The manual also provides a large suite of configured generic examples covering most cases involving tides in a binary system. By "generic," what is meant is that solutions are obtained for the nondimensionalized equations, and processes (e.g., dissipation) are presented as generic classes of mathematical forms with free-ranging coefficients. By considering the range of input parameters and processes, millions of tidal solutions are calculated, and solution domains showing the full range of tidal behavior consistent with the equations are presented. The tidal solutions used in this study are essentially rescaled versions of the generic solutions described in the TROPF manual (see the section Nonsynchronous Arbitrary Rotation). For this study, the tidal responses needed could then in fact be obtained from the generic solution library available in the TROPF package. Here, however, only the input configuration scripts are adopted from the TROPF package, and the tidal solutions are recalculated instantly during the temporal integration of the tidal history. This is practical because TROPF is optimized for speed (as a benchmark, the 1.1 million tidal solutions calculated to spherical-harmonic degree 500 and used to produce Figure 2 take less than 10 minutes to calculate on a 2.4 GHz, 8-core MacBook Pro).
TROPF includes several solution methods. Only the finitevolume method allows for the explicit representation of continents. This study concerns the calculation of tides over Earth's long history where prescription of the continents is not reliable. Because of this, the effects of the continents will be represented parametrically, in which case, the fast sphericalharmonic methods in TROPF can be used. Here we use the is assumed. The overlaid trajectory of the modeled tidal evolution (see Figure 1) shows that the ocean crosses through resonant states three times: (1) very early in the history, near −4.2 Gyr; (2) starting about −500 Myr to peak at −325 Myr and decline at present; and (3) a final resonance that will be encountered about 4.1 Gyr in the future. This figure shows that the amount of resonance in the tidal evolution has been controlled not by the evolution of w , but rather c e 2 (or the Lamb parameter).
primary method (described in the MATLAB® file tropf.m) discussed in the manual. This faster and more versatile method uses a higher-order formulation than one previously used (tropf_sas.m) formulated through the classical Laplace tidal equations with added dissipation and forcing. The previous method was developed by the author as an extension of the approach in Longuet-Higgins (1968) and has been used in a range of applications (e.g., Tyler 2014). The method in tropf_sas.m is also identical to that used in Motoyama et al. (2020). The added versatility of the latest method (tropf.m) is, however, not required in this application concerning a thin, uniform-density ocean. Both methods are based in spherical harmonics, and, as shown in the validation sections of the manual, the two methods provide solutions agreeing to within machine precision. In the case of this study, it is expected that the method used here can be regarded as simply a faster version of the method previously used by this author, as well as in Motoyama et al. (2020). The newer method in tropf.m has also been used in Tyler (2019bTyler ( , 2020, as well as the many configured examples in the manual. The TROPF tidal solution codes are formulated in terms of fully nondimensionalized parameters and variables. This has the efficiency of reducing the inputs needed to the most reduced set. It also provides insight by separating scaling factors from the fundamental nondimensional numbers controlling the behavior of the solution. This separation is especially useful in this study involving parameters that evolve with time because it is important to understand what part of the evolving tidal response is due simply to changes in the scaling factors (e.g., changes in the force amplitude) and what part is due to changes in the behavior of the tidal response (e.g., the level of resonance). In this section, the "fluid tidal response" is regarded as the nondimensional form, and the scaling terms are described in the next section.
The code (tropf.m) providing the nondimensional tidal solution requires the prescription of a number of input parameters falling into three different categories. The description and choices for this application are as follows.
1. Methodological parameters: the truncation degree (N = 500) of the spherical-harmonic expansion; the choice of the nondimensionalized rotation rate (W = 1 ). 2. Tidal force parameters: frequency w (set by the evolving spin/orbit parameters; see Section 2.3); spherical-harmonic order (s = 2); vector G n s describing coefficients for each spherical harmonic (all zero except degree two with prescribed value i/6 such that the projection on the sphere has a maximum amplitude of unity); vectors K n s , d n s , e n s describing other types of (nontidal) forcing set to zero. 3. Fluid response parameters: dissipation constants a d and a r associated with the divergent and rotational components of the flow (more generally, these constants can be entered as vectors describing different values for each degree and therefore harmonic and biharmonic dissipation forms) set to be a a w , where T is the dissipation timescale, treated in this study as a free parameter (i.e., one that is uncertain and must be considered over a range); the potentially complex squared slowness vector taken to have equivalent entries of the scalar parameter n 2 , taken here to be n = c 1 e 2 2˜, where c ẽ is the (real) wave speed treated here as a free parameter.
The justification for these parameter choices is as follows. In the methodological parameters, the truncation degree is much higher than required for this application but consistent with other examples in the manual and also easily feasible given TROPF's speed. The rotation rate is chosen to be unity, as in the examples in the manual (extremely slow rotation cases might use a different choice for numerical stability). The tidal force parameters are set by the description of the force (see Section 2.3), with the only choice left being to choose a unit maximum amplitude for the nondimensionalized tidal potential.
The appropriate choices for the fluid response parameters are less clear. Setting the two dissipation coefficient vectors to be simply equal scalars (i.e., a a = d r˜) is equivalent to specifying a linear drag form. As described in the manual, turbulent viscosity and other dissipation forms can also be as easily prescribed using coefficient vectors rather than the scalars here, but the simple linear drag form seems to be a suitable initial assumption following the approach set out in Section 2.2. Because the coefficients are prescribed to be equal, they can be simply represented with one dissipation timescale T . Similarly, the squared slowness n 2 can be prescribed as a vector (to represent dispersion) and can be complex (the imaginary component representing dissipation by parameterized vertical propagation of energy), but only the simpler real, scalar form seems warranted for this initial study. In this case, n = c describes the ratio of the fluid wave speed c e to the equatorial rotation speed at radius r (Ω is the Earth's rotation rate). In this application involving the thin-shell barotropic tidal response, r can be taken to be a constant representing the average Earth radius, the ocean density can be taken as uniform, and is the shallow-water wave speed (g is gravity and h is the thickness of the uniform ocean shell). One can regard c ẽ as a fundamental parameter expressing the ratio of the speed c e at which the fluid can adjust to disturbances to the equatorial rotation speed 2Ωr (the factor 2 is due to considering the velocity difference across an equatorial diameter). While TROPF is formulated in terms of the parameter c ẽ such that more general cases can be considered (e.g., stratified oceans), for the application treated in this study,c e 2 can be equated with the so-called Lamb parameter (e.g., Longuet-Higgins 1968), which is (in this case) also associated with the eigenvalues of the Laplace tidal equations describing the fluid response.
In this case of forced rather than free oscillations, the frequency of the response w is regarded as being set by the forcing frequency. The response parameter choices then reduce to a prescription of T and c ẽ (or h, using Equation (2)), which will be treated as free parameters. Some further discussion is included below to help in understanding the expected ranges for these parameters.
The dissipation process assumed is of the "kinetic energy" class described in the manual (e.g., Sections 7.1.2, 7.2.1, and 7.3.1), where the dissipation rate is proportional to the average kinetic energy density (as described above, this is consistent with assuming that dissipation follows a linear drag law). The associated nondimensional timescale T can be regarded as the number of tidal oscillations it would take to deplete the kinetic energy if the forcing were turned off. A more typical parameter describing the nondimensional dissipation timescale is the quality factor Q. The notion of Q can become confusing because different definitions have been applied and tides are not an example of a simple harmonic oscillator. In some cases, Q is referred to as a dissipation parameter but calculated from the phase lag of the total work performed (with an implicit assumption of total dissipation and work balance). In electrical engineering, Q can be defined in terms of the width of resonance peaks. Here we use the definition (e.g., as in Egbert et al. 2004 where ω is the (dimensional) forcing frequency, D is the total dissipation rate, and E kinetic + E potential refers to the total timeaveraged energy of the tidal flow. We also assume that total work W and dissipation D balance and can be referred to as the total power P, such that = = W D P. When adopting Equation (3), one sees an additional problem in using Q, at least as a model parameter characterizing dissipation processes. With Equation (3), Q is related to both kinetic and potential energy densities and therefore involves the dynamical partition of energy between gravity and rotational modes (or class I and class II oscillations). Most ocean dissipation processes prescribed in terrestrial tidal models (including this study) are, however, related to the flow and therefore kinetic energy but not the potential energy. Hence, Q is related to not just the dissipation process but also the relative importance of rotation in the dynamics. When comparing two bodies with different rotation rates, Q is therefore not a good parameter for characterizing or discussing differences in their dissipation process or efficiency. Referring to Q as the "dissipation" or "specific dissipation" parameter is therefore confusing for multiple reasons.
That said, Q is indeed useful in describing the damping of a tidal response once calculated. The objection here is only in using Q as a fundamental parameter characterizing the dissipation process. For this reason, in TROPF modeling, Q is considered to be a postprocessed diagnostic variable and not a parameter for characterizing or prescribing the dissipation processes. Within the nondimensional variables used in TROPF, we may also write (following Section 3.3 in the manual) a definition of Q in terms of the nondimensional variables is the nondimensional forcing frequency, Ω is Earth's rotation rate, and   P , , k p˜˜a re the nondimensional time/globe averaged energy and power densities.
Within TROPF, the input parameters controlling dissipation can be understood to be a collection of nondimensional dissipation timescales T j( ) associated with different dissipation processes. Each T j( ) may also be a coefficient vector if the timescale depends on the spatial scale (degree). When one prescribes an equivalent T j( ) for the dissipation of kinetic energy (both the rotational and divergent components) and potential energy, the dissipation rate is proportional to the total energy density, and one can derive the simple relationship = T Q 2 , 6 ( ) such that < T 1 describes the overdamped case, = T 1 describes the critically damped case, and > T 1 describes the underdamped case. Linking the dissipation rate to the flow potential energy can describe the damping of tides by a thick overlying ice layer (see Tyler 2020) but is probably not a good assumption for the case of Earth's ocean tidal history. Instead, dissipation here is linked only to the kinetic energy, as described above, and the effective Q is calculated using Equation (3) once the solution is obtained.
No corrections for self-attraction or loading effects are included in this initial work, but it is expected that improved accuracy can be obtained by including these effects in future modeling. The combined effects of self-attraction and loading can result in a factor of about 0.7 reduction in the effective tidal potential applied (Cartwright 1999), but this depends on the present continentality and Earth interior parameters, and this topic needs to be more carefully explored to decide how it should be parameterized for deep-time integrations such as in this study.
In summary, the calculation of the nondimensional tidal solution requires prescribing a number of input parameters, most of which are straightforward or must follow consistently with the evolving spin/orbit parameters. Once the tidal forcing and generic form for the dissipation process are chosen, the tidal response is dependent on just two parameters that are not certain: dissipation timescale T and ocean thickness h. Within the forcing and model assumptions, the tidal response and therefore the evolution of the Earth-Moon system are then completely controlled by T and h.

Modeling Approach
The proximate goal in the modeling is to reproduce the tidal history with the fewest number of free parameter assumptions, adding more of these degrees of freedom to the model only if the results do not reasonably reproduce the observations. The two parameters c ẽ and T are probably the most reduced set of tunable parameters, since it is not obvious how these should be prescribed even for the present ocean. Evidently, c ẽ and T should represent effective global values incorporating the effects of continents. Note that in Equation (1), c ẽ depends on both an effective ocean thickness h and the Earth's rotation rate Ω, which is one of the parameters that evolve. So the choice of free input parameters to determine is really h and T .
It is expected that h and T can vary over time with changing sea level and continental formation and configuration. But perhaps on the very long timescale of the Earth's history, these time variations average and allow a bulk (or slowly changing) prescription of h and T . The results to be discussed show that the observations can be reasonably matched assuming constant values for h and T , so this paper focuses on this and does not go to the second step of allowing these parameters to change with time.
The task is then to determine which combination of constant h and T is best at explaining the observed data. This will be done by calculating a large set of tidal histories associated with a range of h and T combinations. A "best-fit" example will then be presented and discussed.
In the next two sections, the formulation for the prescription of the tidal force and the method for integration backward through time are described. The parameter values used for the present-day Earth, Moon, and Sun are listed in Table 1.

Tidal Force
Consider the fluid tides excited in a body B s due to the gravity field of a body B p . In the application that will be specifically addressed, we consider the tides raised by the Moon (B p ) in the ocean of Earth (B s ) following simplifying assumptions. (We also simultaneously consider the two-body Earth-Sun system.) In the case of zero obliquity and eccentricity, the degree-two component of the gravitational tide potential given by Section 2.1.20 of the TROPF manual becomes where θ and f are the colatitude and longitude in a B s -fixed frame rotating with frequency Ω, the orbital frequency is Ω o , t is time, P 2 2 refers to the degree and order two Associated Legendre Function, and

= -
The first factor in parentheses in Equation (8) is simply the gravitational potential at the center of B s due to the mass M p of the body B p located at a distance d 0 away (G = 6.6741 × 10 11 m 3 kg −1 s −2 is the gravitational constant). The second factor involves the small ratio r 0 /d 0 (where r 0 is the reference radius of B s ). This factor appears in the first "tidal" component of a McLauren expansion of the gravitational potential. It is a tidal component of the gravitational force because it involves a differential gradient across the body. Lower-order terms in the expansion involve potential fields that are either constant or have a constant gradient (therefore, there is no differential tidal stretch across the body). The third factor can be regarded as a correction factor to account for the possibility that the layer of fluid considered is not at the nominal radius r 0 but rather the radius r. (In this study involving a thin ocean, we take r = r 0 .) The fourth factor is specific to the tidal constituent considered.
An alternative form of Equation (7) is How do we specify Φ 2 in TROPF? TROPF is formulated with a generic, nondimensional forcing potentialG that is related to the dimensional form G by = -G s 1 G G , where G s is an arbitrary scaling factor (or function independent of the horizontal coordinates). For the application here, = -F 2 G , and we select G s = − A, following Equation (8).
The associated input required to representG in TROPF is the vector of spherical-harmonic coefficients G n s , as shown in the expansions (TROPF, Sections 4.2.12 and 4.2.16). For the application here, we select = - , and all other elements in the vector are zero. Once the nondimensional solutions are obtained from TROPF, dimensionalized solutions are produced using our prescribed G s and using the scaling factors described in Table 1 of the TROPF Manual.

Method for Time Integration of the System
For N B tidally interacting bodies in circular, coplanar orbits, ) are the angular momentum and energy of the orbit of system B α , B β ; and τ ( α, β) is the torque on B α due to its tides raised by B β . The rotation rate and moment of inertia of B α are Ω (α) and a I s ( ) . The orbital frequency and separation of system B α , B β are 1.991 × 10 −7 2.6617 × 10 −6 d 0 orbital semimajor axis ( is the reduced mass, and G is the gravitational constant. The moment of inertia for Earth is calculated as described in Webb (1982).
The system in Equations (10) W ( ) and thereby the torques τ ( α, β) ultimately responsible for evolving the systems. Although this formulation allows for the collective effects of tides on the Earth, Moon, and Sun in evolving the system, for the study here, the fluid thicknesses are chosen to vanish on the Moon and Sun, so the only tidal fluid evolving the system is Earth's ocean.

Results
As described, the approach in the modeling is to find the two free parameters h and T (effective ocean depth and nondimensional timescale for dissipation) that best explain the observations of the Earth-Moon tidal history. Given that the observations for the lunar separation d 0 (see Figure 1(b)) have both error bars and the longest time span of records, the fit will seek a tidal history matching these data. The lunar month (Figure 1(f)) is related to d 0 by Kepler's third law and is therefore not an independent validation. All other observed data in Equation (1) do, however, provide some independence in validation. Notably, whereas d 0 directly depends only on the ocean tides raised by the Moon, Earth's rotation rate depends on tides raised by both the Moon and Sun. Hence, agreement with the observations in Figures 1(a) and (d) is additional validation not concomitant with a fit forced on d 0 . Agreement between the model and the current tidal dissipation rate (see Figures 1(c) and (e)) is the weakest requirement because there is little reason to expect that the current tidal dissipation rate represents a longer-term average. Indeed, the current continental configuration may be anomalously effective at blocking westward tidal propagation when compared to previous configurations or times before the formation of continents.
The results (see Figures 3 and 4) show that there is a fairly narrow range around the choices h = 2.3 km and = T 40 that satisfy the observations of d 0 remarkably well. Indeed, with these parameter choices, the model history passes within the error bars of all of the data and also evades placing the Moon inside the Roche limit. The modeled results evade postulating a late formation for the Moon and also fit a history of observations where one can infer a mean, a trend, and at least two inflections (i.e., a nonzero second derivative with respect to time). There is, therefore, little chance that this model with only two free parameters has overfitted the observed data. Even if the model had more tunable parameters, it is unlikely that the data could be spuriously fit because of the physical constraints the model imposes on the types of curves that can be produced. The model also fits the observations better and/or more completely than any previous Earth-Moon history published.
The model also provides reasonable agreement with the other observed data, even though this is not required by the fit. Notably, the LOD and its inflections agree fairly well (Figure 1(a)). The data and model agree on the timing of the peak change in the LOD (Figure 1(d)) about 400 million yr ago. As will be described below, this peak was due to the passage through the first subharmonic resonance.
As seen in Figures 1(c) and (e), the current dissipation rate 2.811 GW, taken to be the sum of the M 2 and S 2 dissipation rates given in Table 1 of Egbert & Ray (2003), is also close enough to the model curve to seem reasonable, especially if it is supposed that the current continental configuration may cause an anomalously higher dissipation rate.
Do these parameters (h = 2.3 km, = T 40 ) found empirically in this study seem reasonable? If we assume that over geological history, the total volume of the ocean has been constant at the present value of about 1.3 × 10 18 m 3 , the removal of continents on a spherical earth would lead to an average thickness h = 2.5 km. Adding continents increases the average ocean depth and therefore, with Equation (2), the wave speed. But the propagation paths should also be typically longer to reach around continents, so this has a canceling effect. How close is this canceling expected to be? Let us first consider the effects of continentality on the wave speed c e , and then we can use Equation (2) to express this in terms of an effective h. Let f c represent the fractional area of the continents on the globe and assume the ocean volume remains as above. If the continents are primarily at high latitudes, then they may offer little obstruction to the fluid waves attempting to propagate the degree-two form at low latitudes (and only the degree-two component of the response is required for calculating the total work). Let us consider the alternative, then, where the continents are all at low latitudes. Considering circular continent(s) separated so as to maximize their obstruction potential, the propagation distance is increased by a factor π/2 over the continent-free case. Combining these two effects, where the first square bracket describes the factor due to increased propagation distance (where the two sets of parentheses respectively describe components that are unaffected and affected by the continental obstruction), and the second square bracket describes the factor due to increased water depth. If we choose f c = 0.3, then c e is increased by a factor of 1.07, and the effective h is a factor of 1.13 over its continent-free value. If continents form over Earth's history, then f c was smaller in the past, and the effects of continentality on increasing the effective h were even smaller. Hence, it seems that the canceling of the two effects of continentality is indeed close following the assumptions here. Therefore, we take as an expected reasonable value h ≈ 2.5 km.
Similarly, T should depend on the continental arrangement, and this may depend on the distribution of basin modes and other factors that are hard to parameterize with unknown continental distributions and ocean dissipation parameters. We have, at least for the current ocean, an estimate that comes from the Q associated with the semidiurnal tides of Q ≈ 18-20 (Egbert et al. 2004). Calculation of Q, using Equation (3), from the model results of this fitted tidal scenario show that over the history, the Q associated with the lunar semidiurnal tide has varied from 24 to 29, while the Q associated with the solar semidiurnal tide has remained close to 20 (see Figure 5(d)). The fitted effective depth h = 2.3 km is remarkably close to the postulated expected value of h ≈ 2.5 km. The dissipation timescale is also reasonably close, especially if the current continental configuration is expected to lead to a Q lower than over most of the history. In short, the empirically found parameters (h = 2.3 km, = T 40 ) indeed seem reasonable. To fit the observations, the range of parameters around the fitted h = 2.3 km and = T 40 is remarkably narrow. In Figures 3  and 4, it can be seen than varying h by just ±100 m or T by ±10 leads to a significant departure in matching the observations. Interestingly, this can, however, lead to a better match for the nonfitted LOD and dissipation rates.

Discussion
To within the confidence of the fitted model, one can now provide a description of the predominant processes that have controlled the evolution of the Earth-Moon system.
The history of tidal dissipation (Figure 1(c)) shows that over most of the history, tidal dissipation has been weaker than at present, despite the decline in the amplitude of tidal forces due to the recession of the Moon (from Equation (8) and the evolution of d 0 , one can show that the tidal forces decrease by a factor of 2.5 over the history). Interestingly, the history was punctuated by a peak about 400 million yr ago when dissipation was an order of magnitude higher than during the 2 or 3 billion yr before or at present. From Figures 1(c) and (e), one can determine that the peak was associated with the Earth rotating with a period of about 21.5 hr. One can see (panel (a)) that this resonance peak is responsible for the observed rapid increase from a 19 to a 24 hr day over the last 500 million yr. A second peak is possible very early in the Earth's history. Note that because the model is integrated backward in time, the results for the recent peak and history are not contingent on accepting as valid the model tidal behavior at this very early period.
The production of the model peak and its agreement with data strongly suggest that this peak is due to resonant enhancement, but why were the conditions for resonance better at −0.4 Gyr? Examination of the model results can make this clear.
The tidal response depends on parameters that enter in two ways. In one, the parameters simply control the scaling factors. In the other, the parameters control the fundamental nondimensional numbers controlling the eigenvalues of the governing equations and therefore the nondimensionalized solutions. Because the latter controls the level of resonance, the nondimensional solutions can be very sensitive to the changes in the input parameters. In calculating the nondimensional tidal response, the only input parameters that evolve with time are the nondimensional forcing frequency w and wave speed c ẽ . By Equation (7), the forcing frequency for the lunar tides is where Ω is Earth's rotation rate and Ω o is the lunar orbital frequency. By Figure 5(a), one sees that the dimensional forcing frequency ω (associated with either the lunar or solar forcing) has evolved significantly, showing a factor of 2 change over the history. However, the fundamental frequency affecting the solution is the nondimensional forcing frequency w , which has remained fairly constant (panel (b)). Hence, the peak in dissipation at −0.4 Gyr is not due to a change in the fundamental forcing frequency.
The ocean fluid responds to the tidal forces with a dimensional wave speed c e , which, in Equation (2), has remained constant over time because both g and h are held constant. However, as seen in Figure 5(c), the more fundamental c ẽ has increased substantially (≈70%) as the Earth's rotation rate Ω has declined.
These results indicate that the levels of resonance and dissipation have evolved over Earth's history primarily due to the effect of rotation on the nondimensional wave speed c ẽ . As described, in this application,c e 2 is equivalent to the Lamb parameter. When the ocean's c e 2 matches one of the eigenvalues of the governing equations, resonant enhancement is achieved.
To show more clearly that it is the evolution of c ẽ and not the changes in forcing frequency that have controlled the tidal response, consider Figure 2. This figure (similar to that presented and discussed in the manual, Chapter 3) summarizes the results from 1.1 million tidal solutions. It presents the nondimensional tidal power as a function of w and c e 2 . (Here = T 40 is prescribed to be that of the fitted value; the effect of decreasing T in this figure is to thicken and reduce the peaks.) The resonance peaks are clearly evident, and most can be described as due to resonantly excited class I rotational-gravity waves. The peak extending to high c e 2 can be described as excitation of a class II Rossby wave. The tidal history has crossed two class I resonance peaks, one (the second class I subharmonic) early on and a second (the first class I subharmonic) more recently. Integrating the model forward in time 4.5 Gyr (see Figure 6), the model expects that the ocean tidal dissipation will continue to decrease for most of the planet's remaining future, with the possibility that the primary class I resonance peak, as well as dissipation stronger than at any previous time, will be encountered in about 4 billion years.
In future work, a larger scope can be considered, where one attempts to use further observational and theoretical constraints. Notably, while within the model framework used here, the tidal response of the Moon need not be considered, the integration method is (partially) set up to allow this. This study has assumed the lunar fluid to have zero thickness and therefore no response, but this can be immediately generalized once a model for the lunar fluid is chosen. The integration method shown does not, however, integrate the changing obliquity and eccentricity, and these are, however, needed to calculate the synchronous tidal forces in the lunar fluid. While adding the lunar fluid tidal response is not expected to significantly alter the model agreement with the observational parameters discussed in this paper, there is an anomalous eccentricity change rate (Williams & Boggs 2016) inferred from observations that could be included in a more complete study.

Conclusions
Earth's rotation rate and the evolution of the Earth-Moon system have been controlled by tidal dissipation in Earth's ocean. Given the importance of understanding this history, there have been many attempts to model this evolution, with results showing incomplete compatibility with observations and unclear isolation of the most important controlling factors. Notably, there has been debate over the relative importance of continentality and the evolving spin/orbital parameters in controlling the tidal dissipation rate. If an accurate prescription of the evolving continents and ocean basins is required, then it might be difficult to form reliable models of the tidal history. This paper has addressed the question of whether it is possible to reliably model the tidal evolution of the Earth-Moon system, given the limited knowledge of past ocean parameters and continental configurations. It is found that a relatively simple model with only two free parameters, h and T , can be fitted to the observations remarkably well. The best-fit values (h = 2.3 km and = T 40 , corresponding to Q = 20-28) are also consistent with the values (h ≈ 2.5 km and Q ≈ 20) expected as reasonable using arguments independent of the model.
Through examination of the model solution, the factors most controlling the tidal evolution are determined. The results show very clearly that dissipation rates have been controlled primarily by the level of resonance of the ocean's response and not the evolution of the amplitude of the tidal forces or other factors. This was shown by demonstrating that the tidal power history follows that of the nondimensional solutions and not that of the scaling factors. The nondimensional solutions evolve with wave speed = W c c t e ẽ˜( ( )) and forcing frequency w w = W W t t , õ˜( ( ) ( )). The wave speed c ẽ increases as Earth's rotation Ω decreases. The forcing frequency varies with both Earth's rotation rate and the orbital frequency, but, as seen clearly in Figure 2, the evolution path has retained w very close to −1. The evolution is almost entirely due to the increase of c ẽ with decreasing rotation rate. Because, in the specific tidal formulation assumed here, c ẽ becomes equivalent to the root inverse of the Lamb parameter, one can alternatively state that the Earth-Moon history has been primarily controlled by the ocean's Lamb parameter, itself evolving primarily with the despinning of the Earth.
In summary, the modeled tidal evolution of the Earth-Moon system has been controlled by the effect of Earth's despinning (due to both the Moon and Sun) on the ocean's nondimensional wave speed c ẽ . The dimensional wave speed c e has remained the same, but the nondimensional form c ẽ is what affects the eigenvalues of the governing equation, and this has increased with the decline of Earth's rotation rate Ω. The results also show that the tidal history has depended much less on changes in forcing frequency w or tidal force amplitudes.
The model presented includes the semidiurnal tidal forces on the ocean from both the Moon and the Sun. The effect of the tides raised by the Sun has primarily been to add a smaller, additional torque despinning the Earth. Interestingly, the Sun still has an effect on the Earth-Moon history because the additional despinning of the Earth has lead to a faster migration toward higher c ẽ than would be obtained by the lunar tides alone.
Given that two free model parameters are fit to the data showing a mean, a trend, and at least two inflections, it is not expected that the agreement of the model with the observations Figure 6. Tidal future of the Earth-Moon system for the fitted input parameters (assuming h = 2.3 km and = T 40 ). The rate of tidal dissipation (and the rate of Earth's despinning) will continue to decline for the next couple of billion years. After 4 billion years, a dissipation peak (associated with an LOD of about 42 hr) will be entered that is larger than any previous one encountered.