The modelled climatic response to the 18.6-year lunar nodal cycle and its role in decadal temperature trends

The 18.6-year lunar nodal cycle arises from variations in the angle of the Moon's orbital plane. Previous work has linked the nodal cycle to climate but has been limited by either the length of observations analysed or geographical regions considered in model simulations of the pre-industrial period. Here we examine the global effect of the lunar nodal cycle in multi-centennial climate model simulations of the pre-industrial period. We find cyclic signals in global and regional surface air temperature (with amplitudes of around 0.1 K) and in ocean heat uptake and ocean heat content. The timing of anomalies of global surface air temperature and heat uptake is consistent with the so-called slowdown in global warming in the first decade of the 21st century. The lunar nodal cycle causes variations in mean sea level pressure exceeding 0.5 hPa in the Nordic Seas region, thus affecting the North Atlantic Oscillation during boreal winter. Our results suggest that the contribution of the lunar nodal cycle to global temperature should be negative in the mid-2020s before becoming positive again in the early 2030s, reducing the uncertainty in time at which projected global temperature reaches 1.5 ^∘C above pre-industrial levels.

1 Introduction

The lunar nodal cycle arises from variations in the angle of the Moon's orbital plane relative to plane of the Earth's Equator (lunar declination), between 18.3 and 28.6^∘, over a period of 18.6 years (Pugh, 1987). A potential connection of this cycle to climate is through the modulation of ocean tides (Loder and Garrett, 1978), the dissipation of which is a major driver of vertical diffusion in the world's oceans (Pease et al., 1995; de Lavergne et al., 2020). The change in lunar declination results in an 18.6-year-period modulation of all lunar and lunisolar tidal constituents and potentially the resulting tidally driven diffusion. The amplitude of the modulation varies depending on tidal constituent, but for the dominant semidiurnal and diurnal constituents (M₂ and K₁) the modulation is small (3.7 % and 11.5 % respectively; Table 1). Previous research has attempted to identify the effect of this signal in climate observations, but since the total modulation of the tide is small, demonstrating a significant effect on global temperature is extremely hard (Ray, 2007). Regional climatic records, such as sea level in regions with large tides, or multi-century proxies, have been shown to exhibit an 18.6-year cycle (Currie, 1984; Yndestad et al., 2006; Yasuda et al., 2006; Gratiot et al., 2008; Agosta et al., 2013; Hamamoto and Yasuda, 2021).

Modelling studies of this phenomenon are relatively rare: simple studies considering modulated stratification in the ocean have suggested an effect on global temperature (Loder and Garrett, 1978). More complex studies involving ocean circulation models (Osafune and Yasuda, 2013) and most recently coupled ocean–atmosphere models have demonstrated an effect of the nodal cycle on the circulation of the Pacific Ocean (Osafune et al., 2020) and suggested a link to variability in the Pacific basin, particularly the Pacific Decadal Oscillation (PDO) (Tanaka et al., 2012; Osafune et al., 2014).

Here we perform millennial length runs of a coupled atmosphere–ocean global circulation model (AOGCM) to quantify the effect of a parameterisation of the lunar nodal cycle on climate. The flexibility of the parameterisation allows for sensitivity tests to be conducted. We investigate the effect of the lunar nodal cycle on long-term trends, with a particular view to understanding its role on the so-called slowdown on global warming in the early part of the 21st century.

2 Method

Our research employs the FORTE2 climate model (Blaker et al., 2021), which uses the primitive equations of meteorology and oceanography on a sphere. The atmospheric component is the IGCM4 and the ocean component is MOMA (Joshi et al., 2015; Webb, 1996). The IGCM4 is run in its full 35-layer stratosphere-resolving configuration, with a horizontal resolution approximating to 2.8^∘. MOMA is run with a 2^∘ horizontal resolution and 15 vertical levels. The background vertical diffusion in the ocean component of the model, a large part of which is accounted for by tidal dissipation, is then modulated using a simple parameterisation that assumes all tidal energy is dissipated locally at the 2^∘ grid scale.

The nodal cycle parameterisation is constructed using the geographical distribution of rms (root mean square) current velocity magnitude for the eight largest tidal constituents (M₂, S₂, N₂, K₂, K₁, O₁, P₁ and M_f), calculated from the TPXO7.2 inverse model (Egbert and Erofeeva, 2002), multiplied by their nodal amplitudes as defined by Pugh (1987; see Table 1). The constituent sum of modulated rms velocity magnitude is divided by the constituent sum of unmodulated rms velocity magnitude, giving the relative modulation of tidal currents at each ocean grid point; such a normalisation is necessary because the parameterised tide in FORTE2, as in most AOGCMs, has constant amplitude in space. Note that modulations of M₂ and N₂ are 180^∘ out of phase with the other tidal constituents, so in regions with a strongly semidiurnal tidal regime (e.g. around New Zealand) the amplitude of the nodal cycle parameterisation may be negative. S₂ and P₁ are pure solar tides so are not directly modulated by the nodal cycle; however, they do contribute to total unmodulated rms velocity magnitude and so affect the relative modulation of tidal currents.

Table 1 Characteristics of the eight tidal constituents using in the parameterisation (Pond and Pickard, 1983; Pugh, 1987). Subscript 2 denotes semidiurnal tides, subscript 1 denotes diurnal tides, and subscript f denotes fortnightly tides. Negative nodal amplitude indicates that the modulation of the constituent is 180^∘ out of phase. Italics indicate a pure solar tide that is not directly modulated by the lunar nodal cycle.

The geographical shape of the function, shown in Fig. 1, is determined by the relative strength of each tidal constituent at a given location and the constituent modulation amplitude. This is multiplied by a normalised 18.6-year sinusoidal cycle to yield a spatially and temporally varying modulation function. The phase of the modulation is such that, at most grid points, tidal currents are at a maximum at 4.75 years into the cycle (e.g. June 2006). The Pacific and Arctic oceans feature modulations of approximately 5 % in amplitude, while the Atlantic Ocean has comparatively little modulation of tidally driven diffusion (typically <2 %). Although the tidal modulation is largest (exceeding 10 %) in the Arctic and Southern oceans, high-latitude water columns are typically only weakly, or negatively, temperature stratified (i.e. the near-surface vertical gradient of temperature is either small or negative). So, counterintuitively, the effect of tidal modulation on climate in these regions might actually be small.

Tides are known to play a controlling role in the energetics of the global ocean, dissipating well over half of the kinetic energy in the oceans, with the greatest dissipation occurring near the ocean floor (Munk and Wunsch, 1998; Egbert and Ray, 2000; St. Laurent et al., 2002). Given the uncertainties in the vertical contribution of tides to the background diffusion, two idealised perturbation runs have been performed – one in which the nodal cycle parameterisation is applied uniformly with depth to the vertical diffusion (“constant”) and one in which it is applied such that its amplitude linearly decreases from 1 at a depth of 5000 m to 0 at the ocean surface (“scaled”) – to mirror the effect of tidal dissipation. The scaled run should be seen as an underestimate of the near-surface effects of the lunar nodal cycle, with the constant run being an overestimate.

Figure 1 (a) Variation in time of the modulation (with reference year for illustrative purposes on the top axis). The tidal modulation in the model is (a) T(t) multiplied by (b) M(x,y), multiplied by 1.0 for the constant run, or a scaled function in the scaled run – see Eq. (1). (b) Geographical distribution of the modulation of tidally driven diffusion by the 18.6-year lunar nodal cycle.

The nodal cycle modulation is applied to the vertical diffusion with a period of 19 FORTE2 years, such that the total diffusion has the form

$$\begin{array}{}\text{(1)}& K^{\prime }=K\cdot T\left(t\right)M(x,y)\cdot S\left(z\right),\end{array}$$

where K is the standard background diffusion in FORTE2 (Blaker et al., 2021), T(t) is the sinusoidal function of Fig. 1a, M(x,y) is the geographically varying function in Fig. 1b, and S(z) is unity for the constant run, or the scaled function described above in the scaled run. Given the length of the year in FORTE2 is 360 d, such an approximation results in a nodal cycle whose length in days is within 0.7 % of the observed cycle. FORTE2 is run for three configurations: pre-industrial control (as in Blaker et al., 2021), scaled, and constant, for 2300 years, with years 1520–2280 being analysed, i.e. 760 years or 40 full cycles.

Figure 2 (a) Globally averaged variation of ocean temperature anomalies (K) vs. ocean depth and phase (years) in the scaled run. The phase coordinate describes the temporal modulation of the vertical diffusion shown in the inset at the top (T(t) in Eq. 1) and in Fig. 1a; i.e. maximum diffusion is in year 5 and minimum diffusion is in year 14. (c) As panel (a) but for the constant run. (b) As panel (a) but for global salinity (PSU). (d) As panel (c) but for global salinity.

3 Results

The global averaged ocean temperature stratification has warm waters in the upper ocean and cooler waters at depth. As the amplitude of the tidally driven diffusion increases in the first phase of the nodal cycle, the global mean vertical temperature gradient is reduced with surface waters cooling and deeper waters warming. The surface temperature anomalies are larger than those at depth, as the vertical temperature gradient is largest in the upper ocean, through the permanent thermocline. Figure 2 shows the evolution of 19-year global ocean temperature and salinity anomalies with depth as a function of the amplitude of the lunar nodal cycle diffusion modulation. In both the scaled case (panel a) and the constant case (panel c), the top 100–150 m of ocean displays a cooling in phase with maximum vertical diffusion in years 4–6. In the absence of any feedback from the atmosphere, the global mean sea surface temperature cold anomaly might be expected to peak at the same time as the subsurface warm anomaly, which is halfway through the nodal cycle in years 9–10. This is when the tidally driven diffusion changes from its enhanced phase to a reduced phase. The response for global salinity is small. Negative salinity anomalies exist in the upper 250 m of the ocean of amplitude of approximately 0.1 PSU in years 8–10 with positive anomalies in years 0–2 (Fig. 2b and d). The salinity anomalies below 250 m are an order of magnitude smaller than those at the surface, reflecting that much smaller vertical salinity gradients are small at depth.

Figure 3 (a) Variation in time of the modulation of diffusion T(t) (see Eq. 1). (b) Globally averaged ocean heat content anomaly (10²² J) vs. tidal modulation phase. The mean ±2 standard errors are shown for the constant run in thin blue and for the scaled run in thin red; sinusoidal best fit curves of global temperature anomalies versus phase are shown for the constant run in thick blue and for the scaled run in thick red; the mean ±2 standard errors of heat content anomaly in the control integration of FORTE2 are shown for reference in black in order to demonstrate the size of the signal compared to internal variability in the control run. (c) As for (b) but for the globally averaged surface ocean heat flux anomaly (W m⁻²) vs. tidal modulation phase.

As shown in Fig. 3c, the atmosphere almost immediately responds to the anomalously cool sea surface temperatures by fluxing heat into the ocean during years 3–7, causing an increase in total ocean heat content between years 3 and 10 (Fig. 3b). The uptake of heat by the ocean results in a global ocean heat content anomaly approximately in quadrature with the surface heat flux; i.e. maximum heat content is in years 9–10 (Fig. 3b), while the maximum surface flux is in years 4–5, or approximately 4.5 years or 90^∘ out of phase with the maximum heat content. The deeper (below approximately 1000 m) temperature anomalies are largely isolated from the surface forcing and are approximately in quadrature with the nodal cycle (Fig. 2). Thus the deep ocean response to the nodal cycle actually lags the response at the surface (see later). As the tidally driven diffusion reduces in the second half of the cycle, the situation described above is reversed.

Figure 4 As Fig. 3 but for (b) globally averaged surface temperature T_surf (K) vs. tidal modulation phase. (c) As for middle panel but for surface temperature (K) in the Arctic region (70–90^∘ N).

Figure 4 shows the effect of the lunar nodal cycle on the model global mean surface temperature T_surf expressed as a function of the phase of the cycle. A best fit of a 19-year harmonic to T_surf shows the phases at which minimum and maximum cooling occur. Minimum global temperatures are reached within a year of the maximum diffusion occurring at year 4.5: such behaviour should be contrasted with the modelled response to transient solar or volcanic forcing, where a lag of approximately 2–3 years is present between maximum forcing and response (Gray et al., 2013). The amplitude of response in the scaled run is 0.03 ± 0.02 K, while the amplitude of the response in the constant run is 0.06 ± 0.02 K. The response in the scaled run is statistically significant only at times of maximum/minimum T_surf, while the larger response in the constant run is significant at more times. Figure 4c shows the response in the Arctic region. Here the pattern is very noisy, but there is some indication of a shift in phase of temperature, which is examined in more detail below.

Figure 5 (a) Geographical variation of the amplitude of the sinusoidal trigonometrical fit to surface temperature (K) whose globally averaged counterpart is shown in Fig. 4b in the scaled run; i.e. yellow colours are where the fitted curve has an amplitude of >0.5 K. Grey shaded areas show where the amplitude is less than 2 standard errors of annually averaged T_surf in the control integration and are used to denote areas where the response is likely to be noise; note nonlinear contour interval. (b) As panel (a) but for the constant run.

We now analyse the geographically varying response of FORTE2 to the lunar nodal cycle. Figure 5 shows the geographical variation of the amplitude of the response of T_surf to the nodal cycle. Figure 5 exhibits quite large responses of amplitude 0.1 K in the northwest Pacific Ocean in both the scaled run (panel a) and the constant run (panel b), consistent with previous work (Tanaka et al., 2012), and the Nordic Seas. Generally, the response in the constant run is larger and statistically significant in more areas than the scaled run, consistent with the larger surface forcing of the nodal cycle in the latter. Interestingly, in both the scaled run and the constant run, the subpolar North Atlantic Ocean displays a significant response of amplitude 0.05–0.1 K. The fitted response is inconsistent with the relatively small tidal forcing in the Nordic Seas (Fig. 1b) and suggests a significant feedback. In addition, the constant run displays a large response on the southern flank of the Southern Ocean of amplitude >0.3 K. Both regions are areas with significant sea ice present, and the robustness of the results in these areas is discussed below.

Figure 6 (a) Geographical variation of the phase of the sinusoidal trigonometrical fit to surface temperature (K) shown in Fig. 5. At each point, the colour denotes the phase in years (see Fig. 1a) where the fitted T_surf reaches its most negative value. Note the cyclic colour interval, since a phase of 19 years is equivalent to a phase of 0 years. Blue colours correspond to phases associated with a globally averaged minimum in T_surf (Fig. 4b). As with Fig. 5, grey shaded areas show where the amplitude is less than 2 standard errors of annually averaged T_surf in the control integration and are used to denote areas where the response is likely to be noise. (b) As panel (a) but for the constant run.

Figure 6 shows the phase in years at which minimum T_surf is reached, with the phase being defined as in Fig. 1a. Note that the globally averaged T_surf exhibits a minimum at years 4–6 (Fig. 4b). Most areas display a minimum T_surf at years 4–6 (blue-purple colours), consistent with Fig. 4. A notable difference is the subpolar northeast Atlantic Ocean, much of the Arctic Ocean, and parts of the Southern Ocean, where a minimum in T_surf occurs in years 14–18, completely out of phase with the global response, in both the scaled run and the constant run (red/orange colours). This polar response can be understood in terms of the local stratification: in the Nordic Seas and Southern Ocean, the ocean temperature maximum occurs at mid-depth rather than at the surface because salinity is the dominant stratifying property. As the internal ocean heat flux is upwards above the temperature maximum, increased (reduced) vertical diffusion associated with the nodal cycle leads to higher (lower) surface temperatures.

The geographically varying phases suggest a potential for geographically varying temperature and circulation responses, especially in the case of the North Atlantic Oscillation (NAO) and Southern Annular Mode (SAM), which are simulated quite well by FORTE2 (Blaker et al., 2021). Figure 7 shows the geographical variation of the amplitude of the response of November–March mean sea level pressure (MSLP) to the nodal cycle. The response is much noisier than for T_surf (Fig. 5), with most areas displaying responses that are not statistically significant (grey shaded). However, there are some regions where significant responses do occur in both the scaled run and the constant run. In particular the northeast Atlantic and Europe exhibit amplitudes exceeding 0.5 hPa, which is consistent with the large amplitude of the T_surf response shown in Fig. 5. The size of these responses in the context of variability forced by other mechanisms is discussed below. There is some indication of a wavelike response at 50^∘ S in the scaled run – however, since there is little sign of such a response in the constant run, which has a larger T_surf response in this region (Fig. 5b), it is unlikely that this response is robust.

Figure 7 As Fig. 5 but for November–March (NDJFM) mean sea level pressure anomaly (hPa).

Figure 8 shows the phase in years at which minimum MSLP is reached, with the phase being defined as in Fig. 1a. In the Atlantic/European region, a significant signal does exist in both the scaled run and the constant run, but the phasing is different. In the scaled run, minimum MSLP occurs over Scandinavia in years 4–6, whereas in the constant run, which has a larger nodal forcing, a statistically significant minimum MSLP occurs over a wider region in years 10–14, somewhat coincident with the phase at which minimum T_surf occurs over the northeast Atlantic Ocean (Fig. 6b). The response implies that in the opposite phase of the nodal cycle, i.e. in years 1–5, positive T_surf anomalies over the northeast Atlantic Ocean and maxima in November–March MSLP (i.e. a negative phase of the NAO) would occur, which is consistent with ideas that link the loss of Arctic sea ice with colder European land temperature in winter (Stroeve et al., 2012). We also note that in the constant run the phasing in years of the minimum in T_surf in northwest Europe is in years 4–6, almost in antiphase to the minimum in T_surf of the adjacent northeast Atlantic Ocean (see Fig. 6b).

Figure 8 As Fig. 6 but for November–March (NDJFM) mean sea level pressure anomaly (hPa).

There is some indication of a dipole across the Pacific Ocean, suggesting an El Niño–Southern Oscillation (ENSO) response, which has been associated with the lunar nodal cycle (Loder and Garrett, 1978; Yasuda, 2018), but the response is not statistically significant. We have analysed surface temperatures in the Nino 3.4 region but find no significant signal. A similar statistically insignificant result is found for the time variation of the Atlantic Meridional Overturning Circulation (AMOC). The size of these responses in the context of other forcings is discussed below.

The long-term effects of the lunar nodal cycle are now examined, since the effect of the ocean circulation would be expected to “redden” a 19-year periodic forcing signal into lower frequencies, measurable on longer timescales. Figure 9 shows decadal T_surf anomalies in each run. There is an increase in the standard deviation of running decadal-mean temperature in the constant run but no clear increase in the scaled run. The presence of the lunar nodal cycle will add a small positive or negative tendency to warming decadal trends in the 21st century.

Figure 9 Histogram of decadally averaged surface temperature anomalies (K) in the control (black), scaled (red), and constant (blue) runs.

Figure 10a shows the assessed projections for global temperature for the five different emissions scenarios used in the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC AR6; Lee et al., 2021). The assessed 5 %–95 % uncertainty range is indicated with shading for the two more extreme SSP (Shared Socioeconomic Pathway) scenarios (SSP5–8.5 and SSP1–1.9; the magnitude is similar for the other scenarios). To highlight the effect of the lunar nodal cycle on these assessed projections, we add a sine wave with peak amplitude of 0.04 K, with the correct lunar nodal cycle timing, to each of the curves (panel b). The lunar nodal cycle is expected to act as a slight cooling influence on the climate in the mid-2020s, delaying the arrival of the 1.5 ^∘C temperature threshold in SSPs with higher carbon emissions (shown in red), but is a warming influence in the early to mid-2030s, hastening the arrival of the 1.5 ^∘C temperature threshold in SSPs with lower carbon emissions (shown in blue). The net effect is to reduce the spread in time at which the world is projected to reach 1.5 ^∘C above pre-industrial levels from 5 to 3 years.

4 Discussion and conclusions

The timing of the lunar nodal cycle is of special interest when considering so-called hiatus and surge decades. A purported slowdown of global temperatures at the start of the 21st century has been much discussed with mechanisms such as volcanic aerosol forcing (Santer et al., 2014) and stratospheric water vapour (Solomon et al., 2010) invoked as part of the explanation, although updated observational datasets show less of a global temperature slowdown than previously identified (e.g. HadCRUT5 and others). Anomalously high heat uptake by the world's oceans (Meehl et al., 2011; Guemas et al., 2013) and circulation changes in the Pacific Ocean (Kosaka and Xie, 2013) are also suggested. Figures 2 and 3 suggest a potential role for the lunar nodal cycle in driving decadal variations in warming rates, with the scaled run implying an average flux of ∼ 0.07±0.07 W m⁻² into the world's oceans over the period 2002–2011. While the uncertainty in the value is clearly large, its magnitude suggests that it cannot be discounted as a significant driver of multidecadal variability of global temperature, given that, for example, the additional heat uptake into the oceans through the surface during hiatus-type periods is approximately 0.7 W m⁻² (Drijfhout et al., 2014).

Figure 10 (a) IPCC AR6 assessed trends in GSAT (global surface air temperature) for different Shared Socioeconomic Pathways (SSPs). (b) As panel (a) but with a lunar nodal cycle of 0.04 K amplitude, chosen to be approximately an average of the constant and scaled runs.

Figure 4 suggests that the contribution of the lunar nodal cycle should be a global cooling of 0.03–0.06 K over the period 2020–2029 and a warming of 0.03–0.06 K over the period 2030–2039. Given the magnitude of such changes, and the results shown in Fig. 10, we suggest that a parameterisation of the lunar nodal cycle should be implemented in 1D integrated assessment models (IAMs) in order to have them better represent the effect of this repeatable and predictable source of climate variability on the impacts of climate change. Although it is known that inclusion of the cycle affects projections of future regional sea level change, for example in the North Sea (Baart et al., 2012), we find that the global modulation of global sea level is 1–2 mm (not shown), because of counteracting influences of hot and cold anomalies in the ocean (Fig. 2). Such a value is far less than the currently observed global sea rise of 3 mm yr⁻¹ (Dangendorf et al., 2019), suggesting that the impact of the lunar nodal cycle on global sea level rise is small.

The geographical response of the model to the lunar nodal forcing can be better understood by putting it in context with other modes of variability. Figure 5 shows that the response of the north Atlantic Ocean has an amplitude of order 0.1 K. For context, this is about 20 %–30 % of the size of sea surface temperature (SST) anomalies associated with Atlantic Multidecadal Variability (Omrani et al., 2022). The results shown in Fig. 7 can be better understood by being put in the context of other sources of variation in the North Atlantic region. The lunar nodal response is smaller but certainly noticeable when compared to the natural variability of the NAO in FORTE2 and observations, which have peak-to-peak amplitudes of 3 and 4 hPa respectively (Blaker et al., 2021). For added context, the response of the NAO to observed Atlantic decadal SST variability is 2–3 hPa (Årthun et al., 2021), while the response to solar variability is 3–4 hPa (Gray et al., 2016), suggesting that the lunar nodal cycle has a much smaller but noticeable effect on Atlantic European winter climate and the NAO.

The Arctic response is almost in antiphase to the rest of the world, reflecting the reversed temperature gradient in the upper ocean, i.e. the lack of a permanent stratifying thermocline. Maximum Arctic temperatures are modelled as occurring in years 5–9 of the cycle, which in reality correspond to 2007–2011, consistent with enhanced Arctic warming during this time (Stroeve et al., 2012). A caveat in interpreting the above results, as well as results suggesting a large response on the northern edge of Antarctica, is that the sea ice representation of FORTE2 is simplified, consisting of one slab (Blaker et al., 2021). Future work regarding the nodal cycle in the subpolar and polar oceans should be carried out with a more realistic sea ice model, with other forcings included to assess potential nonlinear combinations of response. All other things being equal, similar warm Arctic anomalies might be expected during 2026–2030. Figure 7 also implies that the NAO is likely to be more negative than average at the same time.

A caveat in this work lies in the nature of the tracer vertical diffusion scheme which is being modulated. Here we use a simple profile that represents the sum of all diffusion processes that has been tuned to give a good representation of the global thermocline structure. State-of-the-art coupled climate models use a variety of more sophisticated vertical diffusion parameterisations in combination to represent a number of different processes, including wind mixing, tidal mixing, and internal gravity wave scattering (MacKinnon et al., 2017; de Lavergne et al., 2020). Only the tidal and internal tidally induced diffusion is enhanced by the lunar nodal cycle. Thus, if this accounts for one-half of the mixing that we apply (a conservative estimate), then we would expect the magnitude of the response to be halved. An important point to note with regard to parameterisation of the nodal cycle is that tidal forcing is not necessarily at the same time and place where tidal dissipation takes place, implying a limit to the spatial resolution that a parameterisation might employ.

We have implemented a simple, flexible parameterisation of the lunar nodal cycle into an AOGCM; examined its effects on multicentennial-length runs; and assessed its potential effects on 21st century climate. Our results lend further weight to the idea that the phenomenon should be parameterised in decadal-scale forecasts made using global circulation models (e.g. Osafune et al., 2014), as well as in integrated assessment models, given the potential effect of the lunar nodal cycle on future climatic trends.