Near-resonance tidal evolution of the Earth-Moon system influenced by orbital-scale climate change

TL;DR

This study models how Earth's climate variations, especially obliquity changes, influence the tidal evolution of the Earth-Moon system, revealing that climate can significantly slow down tidal evolution and cause phase shifts in orbital oscillations.

Contribution

It introduces a coupled climate-tidal model demonstrating the impact of climate-induced ocean changes on lunar orbital evolution near resonance.

Findings

Climate variations slow tidal evolution compared to constant obliquity.

Earth-Moon distance oscillates in phase with obliquity before resonance, out of phase after.

Displacement of oscillations correlates with oceanic and tidal frequency differences.

Abstract

We build a conceptual coupled model of the climate and tidal evolution of the Earth-Moon system to find the influence of the former on the latter. An energy balance model is applied to calculate steady-state temperature field from the mean annual insolation as a function of varying astronomical parameters. A harmonic oscillator model is applied to integrate the lunar orbit and Earth's rotation with the tidal torque dependent on the dominant natural frequency of ocean. An ocean geometry acts as a bridge between temperature and oceanic frequency. On assumptions of a fixed hemispherical continent and an equatorial circular lunar orbit, considering only the 41 kyr periodicity of Earth's obliquity $\varepsilon$ and the $M_2$ tide, simulations are performed near tidal resonance for $10^6$ yr. It is verified that the climate can influence the tidal evolution via ocean. Compared with the tidal evolution with constant $\varepsilon$, that with varying $\varepsilon$ is slowed down; the Earth-Moon distance oscillates in phase with $\varepsilon$ before the resonance maximum but exactly out of phase after that; the displacement of the oscillation is in positive correlation with the difference between oceanic frequency and tidal frequency.

Tables

Table 1. Table 1 : Values of Model Constants.

Constant	Value	Reference
$A$	218 W m-2	North & Kim (2017)
$B$	1.90 W m-2 K-1	North & Kim (2017)
$D$	0.67 W m-2 K-1	North & Kim (2017)
$C_l$	0.08 yr$\cdot B$	Lin & North (1990)
$C_w$	4.80 yr$\cdot B$	Lin & North (1990)
${\tilde{\alpha }}_0$	0.68	North & Kim (2017)
${\tilde{\alpha }}_2$	$-0.20$	North & Kim (2017)
$\overline{e}$	0.028707	Berger (1978)
$\overline{\epsilon }$	23.333410\dg	Berger & Loutre (1991)
$\Delta \epsilon$	$-1969.00"$	Berger & Loutre (1991)
$\gamma$	31.54068” yr-1	Berger & Loutre (1991)
${\phi }_l$	0\dg	present work
$T_f$	$-10\text{celsius}$	North & Jr. (1979)
${\rho }_{\mathrm{is}}$	0.917 g cm-3	Haynes (2017)
${\rho }_{\mathrm{sw}}$	1.037 g cm-3	Hansen (1982)
$h_{\mathrm{is}}$	2 km
${\overline{h}}_{\mathrm{sw}}$	4 km
$l$	4000 km
$\delta$	0.092	present work

Table 2. Table 2 : Initialization and Results of Numerical simulations.

Case	$\psi$	$r_i$	$r_f$	$2\pi /{\Omega }_i$	$2\pi /{\Omega }_f$	$2\pi /n_i$	$2\pi /n_f$	${\omega }_i$	${\omega }_f$
	($\text{dg}$)	($R_{\oplus }$)	($R_{\oplus }$)	(h)	(h)	(d)	(d)	(${10}^{-4}$ rad s-1)	(${10}^{-4}$ rad s-1)
B0		57.43	57.61	21.44	21.59	25.39	25.51	1.571	1.560
B$+$	$+90$	…	$+4\times {10}^{-6}$	…	$+3\times {10}^{-6}$	…	$+3\times {10}^{-6}$	…	$-2\times {10}^{-7}$
B$-$	$-90$	…	$-6\times {10}^{-6}$	…	$-5\times {10}^{-6}$	…	$-4\times {10}^{-6}$	…	$+4\times {10}^{-7}$
A0		57.80	57.98	21.74	21.89	25.63	25.75	1.549	1.538
A$+$	$+90$	…	$-12\times {10}^{-6}$	…	$-9\times {10}^{-6}$	…	$-8\times {10}^{-6}$	…	$+7\times {10}^{-7}$
A$-$	$-90$	…	$-5\times {10}^{-6}$	…	$-4\times {10}^{-6}$	…	$-3\times {10}^{-6}$	…	$+2\times {10}^{-7}$

Equations

$$C(\phi )\frac{\partial T(\phi ,t)}{\partial t}-\frac{1}{\cos \phi }\frac{\partial }{\partial \phi }[\frac{D}{\cos \phi }\frac{\partial T(\phi ,t)}{\partial \phi }]-[A+BT(\phi ,t)]=W(\phi )\tilde{\alpha }(\phi ).$$

$$C(\phi )=\{\begin{array}{ll}{C}_{\mathrm{l}},& (\phi >{\phi }_{\mathrm{l}})\\ C_{\mathrm{w}}.& (\phi <{\phi }_{\mathrm{l}})\end{array}$$

$$\tilde{\alpha }(\phi )={\tilde{\alpha }}_0+{\tilde{\alpha }}_2{P}_2(\sin \phi ),$$

$$\frac{\partial T}{\partial \phi }=0,(\phi =\pm 90\text{dg})$$

$$T^{\mathrm{s}}({\phi }_{\mathrm{f}})=T_{\mathrm{f}},$$

$$W(\phi ,a,e,\epsilon )=\{\begin{array}{ll}\frac{L_{\odot }\cos \phi }{2{\pi }^3{a}^2\sqrt{1-e^2}}[E(\frac{\sin \epsilon }{\cos \phi })+{\tan }^2\phi K(\frac{\sin \epsilon }{\cos \phi })-{\tan }^2\phi {\cos }^2\epsilon \Pi ({\sin }^2\epsilon ,\frac{\sin \epsilon }{\cos \phi })],& (|\phi |\in [0\text{dg},90\text{dg}-\epsilon ))\\ \frac{L_{\odot }\sin \epsilon }{2{\pi }^3{a}^2\sqrt{1-e^2}}[E(\frac{\cos \phi }{\sin \epsilon })+{\cot }^2\epsilon K(\frac{\cos \phi }{\sin \epsilon })-{\sin }^2\phi {\cot }^2\epsilon \Pi ({\cos }^2\phi ,\frac{\cos \phi }{\sin \epsilon })],& (|\phi |\in (90\text{dg}-\epsilon ,90\text{dg}))\\ \frac{L_{\odot }}{2{\pi }^3{a}^2\sqrt{1-e^2}}[\sin \epsilon +\frac{{\cos }^2\epsilon }{2}\ln (\frac{1+\sin \epsilon }{1-\sin \epsilon })],& (|\phi |=90\text{dg}-\epsilon )\\ \frac{L_{\odot }}{4{\pi }^2{a}^2\sqrt{1-e^2}}\sin \epsilon .& (|\phi |=90\text{dg})\end{array}$$

$$\epsilon (t)=\bar{\epsilon }+\Delta \epsilon \cos (\gamma t+\psi ),$$

$$h(\phi )=\{\begin{array}{ll}1,& (\phi \le {\phi }_{\mathrm{l}})\\ 0.& (\phi >{\phi }_{\mathrm{l}})\end{array}$$

$$h_{\mathrm{sw}}=h_{\mathrm{b}}-\frac{{\rho }_{\mathrm{is}}{h}_{\mathrm{is}}(1-\sin {\phi }_{\mathrm{f}})}{{\rho }_{\mathrm{sw}}(1+\sin {\phi }_{\mathrm{l}})},$$

$$\sigma =\frac{\pi \sqrt{g{h}_{\mathrm{sw}}}}{l},$$

$$I\Omega +M_{\mathrm{r}}n{r}^2=H.$$

$$n^2{r}^3=G(M_{\oplus }+M_{\mathrm{M}}),$$

$$I\frac{\mathrm{d}\Omega }{\mathrm{d}t}=-L.$$

$$L(\sigma ,\omega ,r)=-L_{\ast }\cdot (r_{\mathrm{p}}/r{)}^6\cdot \mathrm{Im}(Z(\sigma ,\omega )),$$

$$L_{\ast }=\frac{6}{5}\pi {\rho }_{\mathrm{sw}}{R}_{\oplus }^2\cdot G{M}_{\oplus }\cdot (M_{\mathrm{M}}/M_{\oplus }{)}^2\cdot (R_{\oplus }/r_{\mathrm{p}}{)}^6,$$

$$Z^{\mathrm{static}}=[⟨|Y{|}^{2}h⟩-|⟨Yh⟩{|}^2/⟨h⟩]/⟨|Y{|}^2⟩,$$

$$\frac{{\mathrm{d}}^2\zeta }{\mathrm{d}{t}^2}+\delta \sigma \frac{\mathrm{d}\zeta }{\mathrm{d}t}+{\sigma }^2\zeta ={\sigma }^2{\zeta }_{\mathrm{m}}^{\ast }{\mathrm{e}}^{\mathrm{i}\omega t},$$

$${\zeta }_{\mathrm{m}}=\frac{{\zeta }_{\mathrm{m}}^{\ast }}{1-(\omega /\sigma {)}^2+\mathrm{i}\delta (\omega /\sigma )}.$$

$$Z(\sigma ,\omega )=\frac{Z^{\mathrm{static}}}{1-(\omega /\sigma {)}^2+\mathrm{i}\delta (\omega /\sigma )}.$$

$$\omega (\Omega ,n)=2(\Omega -n).$$

Full text

1 Introduction

The tidal evolution of the Earth-Moon system is a classic problem, but has not yet been fully solved. The general trend of the tidal evolution has long been confirmed: since the Moon was born 4.5 Gyr ago (Halliday 2008), it has been receding from the Earth with tidal energy dissipated as heat, and Earth’s rotation has been slowing down with angular momentum transferred to the lunar orbit (Murray & Dermott 1999). When one natural frequency (also called a normal mode) of the ocean and one tidal forcing frequency, both varying over geologic time, come close to each other, the ocean tide gets excited and the dissipation is largely enhanced. This phenomenon of ”tidal resonance” speeds up tidal evolution. Currently, the $M_{2}$ (semidiurnal tide) resonance in the ocean contributes a dissipation of 2.4 TW (Munk 1997) to the total of 3.7 TW (Munk & Wunsch 1998), which is so abnormally high that extrapolating into the past unrealistically puts the Moon where it was born as recently as $\sim 2$ Gyr ago (e.g, Touma & Wisdom 1994; Bills & Ray 1999). Because of the uncertainty that tidal resonance brings, to quantitatively reconstruct the lunar orbit, the history of oceanic natural frequencies has to be acquired.

The natural frequencies of the ocean are determined by its geometry (position, shape and depth), which is associated with climate change and continental drift. In order to simulate tidal evolution driven by ocean tide, Hansen (1982) used Laplace’s tidal equations to determine the oceanic tidal torque for two ocean/land geographies, but not only the geography but also the ocean depth remained unchanged in simulations. Webb (1980, 1982a, 1982b) developed a model of average ocean, a statistical average over many hemispherical oceans centered at various positions relative to Earth’s axis, to take the change in ocean geometry due to continental drift into account, whereas the ocean shape and depth were constant. Kagan & Maslova (1994) built a stochastic model, which considered the effect of continental drift as explicit fluctuations in oceanic natural frequencies. By solving the timescale problem (e.g., Goldreich 1966), those ocean modelers demonstrated the importance of the ocean to tidal evolution of the Earth-Moon system, but their results are still qualitative without further improvement for decades, and none has involved climate. The timescale of continental drift is $10^{8}$ yr (Murray & Dermott 1999), and that of orbital-scale climate change is $\sim 10^{5}$ yr (Berger 2012). Considering only continental drift means the secular effect of climate influence can be neglected, which is based on insufficient evidence. Hence, we simulate the tidal evolution for $10^{6}$ yr, so that the influence of climate can be investigated and continental drift can be reasonably neglected.

Natural quasi-periodicities of the climate over timescales in a vast range have been discovered through geological records. On the orbital-scale of $\sim 10^{5}$ yr, the glacial-interglacial cycle dominates (Berger 2012). During a glacial, the ice sheets extend towards the equator and the sea level drops; whereas during an interglacial, the ice sheets shrink towards the poles and the sea level rises. According to the Milankovitch theory, that results from the secularly varying orbit and rotation of the Earth perturbed by the Sun and other planets, and the subsequent variation of insolation distribution on the terrestrial surface (e.g., Berger 1988). This effect on climate of the astronomical parameters (Earth’s eccentricity, obliquity and climatic precession) is called astronomical forcing. The change in sea level, which the glacial-interglacial cycle is accompanied by, alters the natural frequencies of the ocean and then the state of tidal resonance. For instance, in the last glacial maximum $\sim$19–26 kyr ago, the sea level drop of about $130$ m (Yokoyama et al. 2000; Clark et al. 2009; Lambeck et al. 2014) leads to considerably higher dissipation than at present (Thomas & Sündermann 1999; Egbert et al. 2004; Griffiths & Peltier 2009; Green et al. 2009). Therefore, it is the ocean that acts as the bridge between climate and tidal evolution.

In this work, a coupled model of climate and tidal evolution is proposed. An energy balance model (Sellers 1969; Budyko 1969) is applied to simulate the climate in response to astronomical forcing. It is a kind of conceptual model focusing on major climate components and interactions. Though very simplified, it is capable of reproducing the glacial variability with ice sheet involved (Huybers & Tziperman 2008; McGehee & Lehman 2012) and frequently used to study the climate stability (see review in North (1984) for early studies; Lin & North 1990; Wagner & Eisenman 2015). A harmonic oscillator model (Munk 1968; Murray & Dermott 1999) is applied to integrate the lunar orbit and Earth’s rotation with the oceanic natural frequency given. It is also a conceptual model, where the response of the ocean to the tidal forcing is compared to that of a harmonic oscillator, and is capable of providing a realistic timescale of tidal evolution (Hansen 1982; Kagan & Maslova 1994). It suits the case when there is one dominant oceanic natural frequency and one dominant tidal forcing frequency. In addition, a simplified ocean geometry is assumed to obtain the natural frequency from temperature field.

As a preliminary effort to study the influence of climate on tidal evolution, we aim at verifying the existence of the influence and qualitatively observing its nature and mechanism. Therefore, our idealized model and simulation time ($10^{6}$ yr) are appropriate. The period of interest is when the ocean and tidal forcing are near resonance (but not at the resonance maximum), so that the influence can be amplified. Section 2 describes the coupled model and the numerical method, and Section 3 exhibits the results of two sets of simulations in pre- and post-resonance times. A discussion about simulation results and potential improvements is in Section 4.

2 Model and Method

2.1 Climate Model

2.1.1 Steady-state temperature field

The climate can be altered by ocean/land geography. To study climate change resulting from astronomical forcing, a simple geography is used and taken to be invariant. It is assumed that a single spherical-cap continent is centered at the North Pole, extending to latitude $\varphi_{\rm{l}}$, and the rest of Earth’s surface is covered by ocean. Such a geography is similar to what Mengel et al. (1988) used. Neglecting the vertical structure of the atmosphere for this zonally symmetric planet, a one-dimensional climate model can be applied.

Considering horizontal thermal diffusion, outgoing infrared radiation as well as the solar heating being the only external forcing, an energy balance model leads to this governing equation (North & Kim 2017)

[TABLE]

The climate at time $t$ is just characterized by the temperature field $T(\varphi)$ on the surface. In the first term on the left side, $C$ is the effective heat capacity and controls the climate response to perturbations (relaxation time $\tau=C/B$). The capacity over the ocean $C_{\rm{w}}$ is larger than that over the land $C_{\rm{l}}$, and so is the relaxation time for ocean $\tau_{\rm{w}}$ (a few years) than that for land $\tau_{\rm{l}}$ (a month). According to the geography assumed,

[TABLE]

The second term involving the thermal diffusion coefficient $D$ allows the heat transport from warm areas to cool. The third term allows the infrared radiation to space and $A$ and $B$ are empirical coefficients from satellite observations. The term on the right side determines the solar radiation absorbed. The insolation function $W(\varphi)$ gives the latitudinal distribution of the solar radiation flux delivered to the surface. It is dependent on Earth’s orbital status (Sect. 2.1.2). The coalbedo $\tilde{\alpha}(\varphi)$ gives the fraction of radiation absorbed by the surface. Its mean annual form is well represented by

[TABLE]

where constants $\tilde{\alpha}_{0}$ and $\tilde{\alpha}_{2}$ are from satellite observations, and the second-order Legendre polynomial $P_{2}(\sin\varphi)=(3\sin^{2}\varphi-1)/2$. The values of the constants mentioned above and the references that provide them are listed in Table 1. Solving Equation 1 also needs the boundary condition

[TABLE]

implying that there is no net heat flux into the poles.

On the assumption of energy balance (the energy absorbed is equal to the energy lost), for every given insolation function $W(\varphi)$, there exists a corresponding steady-state solution $T^{\rm{s}}(\varphi)$, which any temperature field $T(\varphi)$ of a different profile relaxes to after a time comparable to $\tau$. The relaxation time $\tau$ of the climate system is much smaller than the astronomically driven period of mean annual $W(\varphi)$. Therefore, $W(\varphi)$ can be taken as invariant while $T(\varphi)$ is evolving towards the steady state $T^{\rm{s}}(\varphi)$.

For a steady-state temperature field $T^{\rm{s}}(\varphi)$, the iceline $\varphi_{\rm{f}}$, the edge of the permanent ice cap, is determined by

[TABLE]

where the mean annual isotherm $T_{\rm{f}}=-10\celsius$ (North & Jr. 1979). We note that although the iceline is defined, expressions of the capacity $C(\varphi)$ and coalbedo $\bar{\alpha}(\varphi)$ are not influenced by it, that is, the ice-albedo feedback is not included in the present work. The iceline position only influences the sea level (Sect. 2.2).

2.1.2 Insolation distribution

Given the timescale of interest, the mean annual version of the insolation function $W(\varphi)$ is used, with seasonal variation averaged. It is dependent on Earth’s semimajor axis $a$, eccentricity $e$ and obliquity $\varepsilon$ (Loutre et al. 2004). Assuming the Sun is a point source, Berger et al. (2010) deduced with elliptical integrals the total energy available during any time interval of one year on a given unit surface. Based on their result, we deduce the expression of the mean annual insolation function as

[TABLE]

The solar luminosity $L_{\odot}$ is a constant, and $K$, $E$ and $\Pi$ are the first, second and third complete elliptical integrals, respectively. This expression is valid for $0\dg<\varepsilon<90\dg$ and $0<e<1$.

Variation of $W(\varphi)$ results from those of the astronomical parameters $e$ and $\varepsilon$ (mainly over a timescale of $\sim 10^{5}$ yr). Based on the numerical solution for Earth’s orbit (e.g, Laskar 1988), $e$ and $\varepsilon$ can be expressed in trigonometric form as quasi-periodic functions of $t$: ${\textrm{approximation}}+\sum{\{{(\textrm{amplitude}})_{i}\cos{[{(\textrm{frequency}})_{i}t+{(\textrm{phase})}_{i}]}\}}$ (Berger 2012). In this work, we only consider the most important term of $\varepsilon$, whose period is 41 kyr, in order to avoid the compound influence of its multiple terms and the complicated effect of $e$ and $\varepsilon$ simultaneously varying. Thus, $e=\bar{e}$ and

[TABLE]

where values of the approximation $\bar{\varepsilon}$, the amplitude $\Delta\varepsilon$ and the frequency $\gamma$ given by Berger & Loutre (1991) are used, which are valid over 1–3 Myr (Berger & Loutre 1992). However, the phase $\psi$ acts as a controllable parameter in our simulations to manifest the influence of itself.

2.2 Ocean Geometry Model

An ocean geometry model is needed to connect climate and tidal evolution. An ocean function with a value of 1 over ocean and 0 over land can be defined to characterize the geography modeled in Section 2.1

[TABLE]

The iceline $\varphi_{\rm{f}}$ divides the water on Earth’s surface into two reservoirs, i.e., sea water in the ocean basin and ice sheet on the part of continent north of the iceline. It is further assumed that the depth of the ocean basin $h_{\rm{b}}$, the depth of the sea water $h_{\rm{sw}}$ and the thickness of the ice sheet $h_{\rm{is}}$ are all uniform. Sea ice is not considered. By simply taking the volume of sea water as its depth times the area of the ocean basin and the volume of the ice sheet as its thickness times the area of ice cover, conservation of mass gives rise to

[TABLE]

where $\rho_{\rm{is}}$ and $\rho_{\rm{sw}}$ are densities of ice sheet and sea water, respectively, and $h_{\rm{b}}$ is equal to the maximum of $h_{\rm{sw}}$ which corresponds to the condition $\varphi_{\rm{f}}=90\dg$.

The oceanic natural frequency $\sigma$ is determined by the geometry of the ocean. On the assumption of half-wavelength resonance, the ocean basin is simplified as a closed square and $\sigma$ can be estimated to be

[TABLE]

where $l$ is the ocean width and $g$ is the gravitational acceleration. Typical values of $h_{\rm{is}}$, $\bar{h}_{\rm{sw}}$ (mean of $h_{\rm{sw}}$) and $l$ are set (Table 1). The advantage in setting $l$ to a typical value instead of deriving it from the modeled geography is that a realistic $\sigma$ value can be obtained. Because the geography is invariant, the oceanic frequency is only dependent on climate, i.e., $\sigma$ only varies with $h_{\rm{sw}}$. On the other hand, $\sigma$ is yet needed in the calculation of tidal torque coefficient $Z$ (Sect. 2.3.2). The climate and tidal evolution are thus connected.

2.3 Tidal Evolution Model

2.3.1 Orbital parameters

A two-body system consisting of the Earth and Moon with a circular orbit in the Earth’s equatorial plane is considered. With the Moon’s rotational angular momentum neglected, conservation of angular momentum leads to

[TABLE]

The first term on the left side is Earth’s rotational angular momentum, where $I$ is Earth’s rotational inertia and $\Omega$ is Earth’s rotational speed. The second term is the lunar orbital angular momentum, where $n$ is the lunar angular orbital speed, $r$ is the Earth-Moon distance, and the reduced mass $M_{\rm{r}}=M_{\oplus}M_{\rm{M}}/(M_{\oplus}+M_{\rm{M}})$ ($M_{\oplus}$ and $M_{\rm{M}}$ are masses of Earth and Moon, respectively). The total angular momentum $H$ is constant. There are three orbital parameters, $r$, $n$ and $\Omega$, characterizing the state of tidal evolution. Besides Equation 11, $n$ and $r$ are also linked by Kepler’s third law,

[TABLE]

where $G$ is the gravitational constant. Therefore, knowledge of one of $r$, $n$ and $\Omega$ is equivalent to that of them all.

The evolution of $\Omega$ is determined by tidal torque $L$ (Sect. 2.3.2), which arises because the Earth, carrying its tidal bulge, rotates faster than the Moon orbiting it ($\Omega>n$)

[TABLE]

The tidal torque acts to decrease $\Omega$, resulting in a transfer of angular momentum from Earth’s rotation to lunar orbital motion and a dissipation of energy in Earth.

2.3.2 Tidal torque

A tide is raised on an elastic body when this body is distorted in the gravity field of another. In our model, the solid Earth and Moon are taken as rigid spheres, while the ocean is a thin deformable layer partially covering the solid Earth. Therefore, only the tidal torque exerted on the Moon by the distorted ocean is present. Furthermore, only the semidiurnal tide $M_{2}$, the dominant tidal constituent at present, is considered in this work for simplification.

Hansen (1982) derived an approximate expression of the secular variation of the tidal constituent torque after averaging on a short timescale (monthly and yearly) and neglecting terms higher than $(R_{\oplus}/r)^{8}$. For the assumed equatorial lunar orbit, that expression becomes

[TABLE]

where the constant $L_{*}$, which carries the dimensionality of the torque, is

[TABLE]

$r_{\rm{p}}$ is the present Earth-Moon distance and $R_{\oplus}$ is the radius of Earth.

Additionally, the torque coefficient $Z$ in Equation 14 is a second-degree spherical harmonic expansion coefficient of the complex tidal elevation response function. For a static ocean tide whose shape is the same as that of the tide potential, $Z$ degenerates into a static torque coefficient

[TABLE]

where the complex spherical harmonic $Y=P_{2}^{2}(\cos\theta){\rm{e}}^{{\rm{i}}2\lambda}$ for $M_{2}$ tide (the unnormalized associated Legendre function $P_{2}^{2}(\cos\theta)=3\sin^{2}\theta$, $\theta$ is colatitude and $\lambda$ is longitude), the ocean function $h$ is defined in Equation 8, and the angled brackets imply an areal integration over the global surface.

For a dynamic ocean tide, following Hansen (1982), the harmonic oscillator model is adopted to derive $Z$. A driven harmonic oscillator can be described as

[TABLE]

where $\zeta$ is the displacement from equilibrium, $\delta$ is the frictional resistance coefficient, $\sigma$ is the natural frequency of the oscillator, $\omega$ is the frequency of the external force and $\zeta_{\rm{m}}^{*}$ is the limit of maximal displacement as $\omega$ approaches 0. Its steady-state solution is $\zeta=\zeta_{\rm{m}}{\rm{e}}^{{\rm{i}}\omega t}$, where the displacement amplitude is

[TABLE]

In an analogy with ocean tides, $\zeta$ is taken as the tidal elevation in the ocean, $\sigma$ is the oceanic natural frequency, $\omega$ is the tidal forcing frequency, $\zeta_{\rm{m}}^{*}=Z^{\rm{static}}$ and $\zeta_{\rm{m}}=Z$. The tidal torque coefficient is thus derived

[TABLE]

This tidal torque coefficient $Z$ varies with $\sigma$ and $\omega$, for $Z^{\rm{static}}$ absolutely determined by geography is constant in our model. The oceanic natural frequency $\sigma$ is dependent on the state of ocean (Sect. 2.2), while the tidal forcing frequency for $M_{2}$ tide is

[TABLE]

If $\sigma$ and $\omega$ are close enough (not strictly equal for a nonzero $\delta$), ${\rm{Im}}(Z)$ will be largely enhanced and so will $L$. Tidal evolution then speeds up, and that is when a tidal resonance is considered to occur.

2.4 Method of Solution

The values of constants involved in our model are listed in Table 1. Because of the simplification of our model, whether their values are precise or not does not affect our qualitative conclusions. Numerical simulations are executed for different initial Earth-Moon distances and phases of Earth’s obliquity. The initialization of simulations will be presented in Section 3. The following is the procedure carried out for any instant after the initial time as shown in Figure 1.

In every calculation loop characterized by $t$, obliquity $\varepsilon$ is the first to be obtained. The term with the biggest amplitude 0.547\dg(Berger & Loutre 1991) in its trigonometric expansion is used to get $\varepsilon$ (Eq. 7). The corresponding period is 41090 yr. Then, the insolation function $W(\varphi)$ is derived from $\varepsilon$ using Equation 6, with $a$ fixed at its present value and $e$ fixed at its approximation over the last few million years $\bar{e}$.

An numerical method is applied to derive the steady-state temperature field $T^{\rm{s}}(\varphi)$ from $W(\varphi)$. Specifically, the differential equation of $T$ (Eq. 1) is discretized by centered finite difference method. The latitude is discretized in intervals of $\Delta\varphi=\pi/180$, while the time step is given by $\Delta t/2$, where the propagation time $\Delta t=(\Delta\varphi^{2}/2)(C_{\rm{l}}/D)$. Stepping forward in time from the initial condition $T(\varphi)=10\celsius$, the iteration does not cease until the relative error of temperature is less than $10^{-6}$, which happens after no more than $3.5\tau_{\rm{w}}=17$ yr in our simulations. Thus, a steady-state solution is considered to be reached. Because it is instantly reached compared to the time step in integration for tidal evolution ($\sim 10^{3}$ yr), the steady-state solution $T^{\rm{s}}(\varphi)$ solved with $W(\varphi)$ given at time $t$ is just taken as the temperature field at that instant. We note that the topography adopted is a hemispherical continent ($\varphi_{\rm{l}}=0\dg$), and the nonlinearity in $C(\varphi)$ does not introduce any convergence problems in our procedure.

Defining $T_{\rm{f}}=-10\celsius$, the iceline $\varphi_{\rm{f}}$ (Eq. 5) is then quickly located by linearly interpolating the adjoined latitudes where temperatures are found to just enclose $T_{\rm{f}}$. The sea water depth $h_{\rm{sw}}$ is calculated from $\varphi_{\rm{f}}$ using Equation 9, where the maximal depth $h_{\rm{b}}$ is set (in the first loop) to what ensures that the mean depth $\bar{h}_{\rm{sw}}=4$ km (equal to $h_{\rm{sw}}$ in the first loop, given $\psi=\pm 90\dg$) during simulations. The oceanic natural frequency $\sigma$ is calculated from $h_{\rm{sw}}$ using Equation 10, where the ocean width is set to $l=4000$ km so that the mean natural frequency $\bar{\sigma}=1.555\times 10^{-4}$ rad s*-1* is near the present $M_{2}$ resonance.

After $\sigma$ is obtained, a Runge-Kutta-Fehlberg method is applied to integrate the differential equation of Earth’s rotational speed $\Omega$ (Eq. 13). As pointed out in Section 2.3.1, the lunar orbital parameters $\Omega$, $n$ and $r$ are related by Equation 11 and 12, so $n$ and $r$ can be known with $\Omega$ given. In Equation 11, the total angular momentum $H=3.442\times 10^{34}$ kg m2 s*-1* which is determined by substituting $\Omega$, $n$ and $r$ with present values. The tidal frequency $\omega$ is then also known because of its dependence on $\Omega$ and $n$ (Eq. 20). Therefore, within each time step of the Runge-Kutta-Fehlberg integrator, while $\sigma$ is held constant, $L$ can be obtained with $r$ and $\omega$ by using Equation 16, 19 and 14 in turn, following the calculation of $r$ and $\omega$ from the coinstantaneous $\Omega$. In these equations, the static torque coefficient $Z^{\rm{static}}=0.5$ for a hemispherical continent, the frictional coefficient $\delta$ is set to 0.092 in order to ensure a realistic timescale of tidal evolution based on some test simulations and the dimensionality is calculated to be $L_{*}=1.998\times 10^{17}$ N m.

Given $\Omega$ at $t$, what the Runge-Kutta-Fehlberg integrator finds is $\Omega$ at the next instant $t^{\prime}$. Updating $n$, $r$ and $\omega$ using $\Omega$ at $t^{\prime}$ immediately follows. Updating the climate and ocean status at $t^{\prime}$ starts in the next loop. The above procedure is repeated for every instant until the final time.

3 Results

3.1 Pre-analysis of tidal evolution

We first present a pre-analysis of the general trend of tidal evolution based on our model, for that helps in explaining the initialization of the numerical simulations. According to our tidal evolution model, Earth’s rotational speed $\Omega$ and the lunar orbital speed $n$ can be expressed as functions of the Earth-Moon distance $r$, and so can the tidal frequency $\omega$ (Sect. 2.3). If the oceanic natural frequency $\sigma$ is constant, the resonance distance $r_{\rm{res}}$, where the tidal resonance occurs as $\omega\approx\sigma$ (for a nonzero dissipation), can be then predicted without simulations.

Figure 2 shows that as $r$ increases from $10\,R_{\oplus}$, both $\Omega$ and $n$ decrease, which means both one day $2\pi/\Omega$ and one month $2\pi/n$ lengthens. Their difference diminishes until $\Omega=n$, where Earth begins to synchronously rotate (one day is as long as one month) just like the Moon has been doing in reality (though Moon’s rotation is neglected in our model) and the tidal evolution ends. According to the constants we set, that happens at $r_{\rm{syn}}=86.9\,R_{\oplus}$. In addition, as $r$ increases, $\omega$ also keeps decreasing until $\omega=0$, when the resulting tidal torque and dissipation become zero. If the oceanic frequency is fixed at $\bar{\sigma}=1.555\times 10^{-4}$ rad s*-1*, the tidal resonance is found to occur at $r_{\rm{res}}\approx 57.7\,R_{\oplus}$, slightly smaller than the present Earth-Moon distance $r_{\rm{p}}=60.3\,R_{\oplus}$. This prediction is in accordance with the reality that the oceanic response has been near $M_{2}$ resonance currently.

Near $r_{\rm{res}}$, the tidal evolution should greatly speed up with dissipation of the total energy enhanced largely. The rapid decrease in $\omega$ at that time gives rise to the rapid pass through resonance. However, if $\sigma$ is varying, conditions become complicated. Before $r_{\rm{res}}$, a decreasing $\sigma$ delays the resonance while an increasing $\sigma$ hastens it; after $r_{\rm{res}}$, the former extends it while the latter curtails it. Even multiple passes can arise. Those conditions are beyond the scope of the present work. Therefore, although we focus on the near resonance condition, we will simulate periods a while before and after $r_{\rm{res}}$, instead of simulating the maximum period.

3.2 Numerical simulations

The initialization for our numerical simulations is shown in Table 2. Two sets of simulations, Cases B and Cases A (”before” and ”after” resonance maximum respectively), are performed with initial Earth-Moon distance $r_{\rm{i}}=57.43$ and 57.80 $R_{\oplus}$, respectively, which are slightly smaller and larger than the resonance distance $r_{\rm{res}}\approx 57.7\,R_{\oplus}$ (Sect. 3.1). The initial values of the other two orbital parameters $\Omega_{\rm{i}}$ and $n_{\rm{i}}$ are derived from $r_{\rm{i}}$ (Sect. 2.3.1). The corresponding initial tidal frequency $\omega_{\rm{i}}=1.571\times 10^{-4}$ and $1.549\times 10^{-4}$ rad s*-1* for Cases B and A, respectively, just enclosing $\bar{\sigma}=1.555\times 10^{-4}$ rad s*-1*, i.e., the oceanic frequency in the case of invariant $\varepsilon$. Each set includes three cases: Case 0 acts as a control group with constant obliquity $\bar{\varepsilon}=23.33\dg$; for Cases $+$ and $-$, $\varepsilon$ periodically varies with $\psi=+90\dg$ and $-90\dg$, respectively. The influence of the existence of climate change can thus be examined. All the other parameters are set as described in Section 2.4.

We note that our aim in this work is to study the mechanism of climate influence but not the realistic history of lunar orbit. Therefore, although not comparable to the timescale of tidal evolution $10^{9}$ yr, the simulation time $10^{6}$ yr fits this aim, and it is not a problem that the cases with varying $\varepsilon$ are initialized at the same $r_{\rm{i}}$ as the case with constant $\varepsilon$. Furthermore, it is therefore reasonable to randomly set the phase of obliquity $\psi$ at any $r$. We choose $+90$ and $-90\dg$ in order to maximize the difference between Cases $+$ and $-$.

Cases B start and end before the resonance maximum (Table 2). As shown in Figure 3a, $\varepsilon$ increases and decreases from $\bar{\varepsilon}$ at the beginning for Case B$+$ and Case B$-$, respectively, and periodically varies till the end. According to Loutre et al. (2004), the mean annual insolation $W$, being symmetrical between northern and southern hemispheres, varies in phase with $\varepsilon$ in the high latitudes but exactly out of phase in the low latitudes. These properties are completely exhibited in Figure 3c and 3d. The steady-state temperature $T^{\rm{s}}$, as a response to insolation, has the same properties as $W$. As shown in Figure 3e and 3f, located in the high latitudes, the iceline $\varphi_{\rm{f}}$ varies in phase with $\varepsilon$ for both Cases B$+$ and B$-$. Its mean $\bar{\varphi}_{\rm{f}}=65.44\dg$ and the amplitude defined as the maximal deviation from the mean is $0.23\dg$. Because the sea water depth $h_{\rm{sw}}\propto\sin\varphi_{\rm{f}}$ (Eq. 9) and the oceanic frequency $\sigma\propto\sqrt{h_{\rm{sw}}}$ (Eq. 10), they are in phase with $\varepsilon$ as well. For either Case B$+$ or B$-$, $h_{\rm{sw}}$ oscillates about $\bar{h}_{\rm{sw}}$ with an amplitude of 2.98 m (Fig. 3b), and $\sigma$ oscillates about $\bar{\sigma}$ with an amplitude of $0.001\times 10^{-4}$ rad s*-1* (Fig. 4a).

The general trend of tidal evolution with constant $\sigma$ has been interpreted in Section 3.1. We now list the initial and final values of the lunar orbital parameters in Table 2 and illustrate only the divergence of Earth-Moon distance $\Delta r$ for Cases B$+$ and B$-$ from Case B0 in Figure 4b in order to focus on the climate influence. Three features are commonly observed for Cases B$+$ and B$-$. First, $\Delta r$ varies in phase with $\sigma$ and thus in phase with $\varepsilon$. The reason is that during the pre-resonance time when $\omega>\sigma$, greater $\sigma$ means greater tidal dissipation, resulting in a leading evolution characterized by a larger $r$. Second, the equilibrium point of the $\Delta r$ oscillation is not constant but seems to decrease at least on the near side of resonance maximum. Third, the displacement of $\Delta r$ from the gradually decreasing equilibrium point expresses a positive correlation with the difference between $\omega$ and $\sigma$. One distinction between Cases B$+$ and B$-$ is that the mean of $\Delta r$ for the former is larger than the latter. We attribute this distinction to $\psi$ and the beginning behavior of $\Delta r$ it determines: the beginning increase/decrease in $\Delta r$ for Case B$+$/B$-$ contributes to a lead/drop lasting for the whole simulation time.

Cases A start later than the resonance maximum (Table 2). The evolutions of climate and ocean state ($\varepsilon$, $W$, $T^{\rm{s}}$, $\varphi_{\rm{f}}$, $h_{\rm{sw}}$ and $\sigma$) for Cases A are totally the same as those for Cases B. However, because $\omega<\sigma$ during the post-resonance time (Fig. 5a), contrary to Cases B, greater $\sigma$ means smaller dissipation and thus a smaller $r$. As shown in Figure 5b, $\Delta r$ for either Case A$+$ or A$-$ is therefore exactly out of phase with $\sigma$. The second and third features for Cases B, i.e., the general trend of decreasing and positive correlation between $\Delta r$ displacement and the difference between $\omega$ and $\sigma$, also match Cases A. In addition, $\psi$ again results in a larger mean of $\Delta r$ for the case where $\Delta r$ increases at the beginning (A$-$) than where it decreases (A$+$). It is worth mentioning that both as the case with the larger mean $\Delta r$, B$+$ holds a positive $\Delta r$ for the whole time, whereas A$-$ holds for only about $10^{5}$ yr. Considering the general decreasing trend, the fact that B$+$ holds a lead over B0 is probably a temporary effect.

In summary, the features of our climate-influenced tidal evolution (characterized by $\Delta r$ whose positive value means a lead over the evolution with climate unchanged and negative value means a lag) are

Given that iceline $\varphi_{\rm{f}}$ is in high latitudes (so that $\sigma$ is in phase with $\varepsilon$), $\Delta r$ varies in phase with $\varepsilon$ during the pre-resonance time but exactly out of phase during the post-resonance time. 2. 2.

Despite oscillation, the general trend of $\Delta r$ is decreasing. 3. 3.

The displacement of $\Delta r$ oscillation is in positive correlation with the difference between $\omega$ and $\sigma$. 4. 4.

As a whole, $\Delta r$ oscillation is shifted upwards or downwards as $\Delta r$ increases or decreases at the beginning.

4 Discussion

Based on our conceptual coupled model of climate and tidal evolution (Sect. 2), we carried out numerical simulations of the near-resonance tidal evolution for an equatorial circular lunar orbit with Earth’s obliquity $\varepsilon$ periodically varying (Sect. 3.2). Thus, the climate influence on the tidal evolution via ocean is verified. Our conclusions in terms of the influence mechanism are qualitative. The main conclusion is that compared to the case that the climate is invariant, varying climate slows down the evolution accompanied by oscillations. Furthermore, the oscillation is in phase and exactly out of phase with $\varepsilon$ before and after the resonance maximum, respectively; and can be enlarged as the difference between the oceanic frequency $\sigma$ and the tidal frequency $\omega$ increases.

The above conclusions should be applied with caution. This is not only because of the idealization and the existence of multiple parameters of the model, but also because the simulations are only done for a short instant near the resonance maximum of the whole lunar tidal evolution.

Though we focus on the mechanism of climate influence in this work, it should be pointed out that the absolute differences in final orbital parameters found between the cases with varying and invariant climates (Table 2) are insignificant indeed. However, it is still hasty to conclude that the influence of climate change can be neglected. On one hand, the timescale studied here, $10^{6}$ yr, is very short compared with the timescale of tidal evolution, $10^{9}$ yr. If the evolution keeps slowing down with climate varying, the secular accumulation may make a difference. On the other hand, the variations of the climate and ocean state produced here are not as big as in reality. The maximal drop of the sea water depth in simulations is 6 m, whereas the sea level drop during the last glacial maximum relative to the present is about 130 m (Clark et al. 2009). If a more realistic model is used, the influence will be enhanced.

One important effect that can enhance the climate influence is the ”ice-albedo feedback.” In the current model, though the ice sheet on the continent is considered, the coalbedo $\tilde{\alpha}$ has no dependence on iceline $\varphi_{\rm{f}}$. A more realistic way, for example, is to multiply $\tilde{\alpha}$ by $1/2$ in latitudes higher than $\varphi_{\rm{f}}$ (Mengel et al. 1988). In this case, as the ice cover spreads, the planetary coalbedo and thus the absorbed solar radiation diminishes, leading to a further drop in temperature accompanied by the spread of ice cover. In other words, a slight change in solar radiation can cause an abrupt climate transition (North 1984). This nonlinear feedback mechanism will be introduced in our future model.

Potential subjects of future works include improving the model, determining the quantitative correlation between climate variation and the rate of tidal evolution, and generalizing the model to other planet-satellite systems. In addition, considering that it is the normal modes of the liquid part of the Earth that can largely be excited when tidal resonance occurs, the tidal evolutions of terrestrial planets perturbed by companions in exosolar systems (e.g., Dong & Ji 2013) may also need further investigations when oceans or liquid cores (Liu & Li 2018) are present.

Acknowledgements.

We thank the referee for constructive comments. We are grateful that Professor Kwang-Yul Kim helped in climate simulations. The code of tidal evolution was developed in Nanjing University when the author Nan Wang was under the supervision of Professor Ji-Lin Zhou, and further improved in Zhejiang University. This work was funded by the Natural Science Foundation of Zhejiang Province (LR16E090001), the National Key Research and Development Program of China (Grant No. 2017YFC0305905), and NSFC-Zhejiang Joint Fund for the Integration of Industrialization and Informatization (Grant No. U1709204).