Calculate Earth’s Obliquity and Precession in the Past

snvec() computes climatic precession and obliquity (or tilt) from an astronomical solution (AS) input and input values for dynamical ellipticity (\(E_{d}\)) and tidal dissipation (\(T_{d}\)). It solves a set of ordinary differential equations.

Usage

snvec(
  tend = -1000,
  ed = 1,
  td = 0,
  astronomical_solution = "full-ZB18a",
  os_ref_frame = "HCI",
  os_omt = NULL,
  os_inct = NULL,
  tres = -0.4,
  atol = 1e-05,
  rtol = 0,
  solver = "vode",
  quiet = FALSE,
  output = "nice"
)

Arguments

tend

Final timestep in thousands of years (kyr). Defaults to -1000 kyr.

ed

Dynamical ellipticity \(E_{d}\), normalized to modern. Defaults to 1.0.

td

Tidal dissipation \(T_{d}\), normalized to modern. Defaults to 0.0.

astronomical_solution

Character vector with the name of the desired solution. Defaults to "full-ZB18a".

os_ref_frame

Character vector with the reference frame of the astronomical solution. Either "HCI" for heliocentric inertial reference frame or "J2000" for ecliptic J2000 reference frame. Defaults to "HCI" for HNBody output.

os_omt

Longitude of ascending node of the solar equator relative to ECLIPJ2000.

os_inct

Inclination of the solar equator relative to ECLIPJ2000.

tres

Output timestep resolution in thousands of years (kyr). Defaults to -0.4. To determine the sign, think of from 0 to tend by timestep tres.

atol

Numerical absolute tolerance passed to deSolve::ode()'s atol. Defaults to 1e-5.

rtol

Numerical relative tolerance passed to deSolve::ode()'s rtol. Defaults to 0.

solver

Character vector specifying the method passed to deSolve::ode()'s method. Defaults to "vode" for stiff problems with a variable timestep.

quiet

Be quiet?

If TRUE, hide info messages.
If FALSE (the default) print info messages and timing.

output

Character vector with name of desired output. One of:

"nice" (the default) A tibble with the columns time, eei, epl, phi, cp.
"full" A tibble with all the computed and interpolated columns.
"ode" A matrix with the output of the ODE solver.

Value

snvec() returns different output depending on the outputs argument.

If output = "nice" (the default), returns a tibble with the following columns:

time Time in thousands of years (kyr).
epl Calculated Obliquity \(\epsilon\) (radians).
phi Calculated Precession \(\phi\) (radians) from ECLIPJ2000.
cp Calculated Climatic precession (-) as \(e\sin\bar{\omega}\).

where \(\bar{\omega}\) is the longitude of perihelion relative to the moving equinox.

If output = "all" (for developers), additional columns are included, typically interpolated to output timescale.

sx, sy, sz The \(x\), \(y\), and \(z\)-components of Earth's spin axis unit vector \(\vec{s}\) in the heliocentric inertial reference frame.

See the source code for descriptions of all the intermediate computational steps.

If output = "ode", it will return the raw output of the ODE solver, which is an object of class deSolve and matrix, with columns time, sx, sy, and sz. This can be useful for i.e. deSolve::diagnostics().

Details

This is a re-implementation of the C-code in the supplementary information of Zeebe & Lourens (2022). The terms are explained in detail in Zeebe (2022).

Reference Frames of Astronomical Solutions

NASA provides their asteroid and planet positions in the ecliptic J2000 reference frame, while long-term astronomical solution integrations are often performed in the heliocentric inertial reference frame (HCI) or in the inertial reference frame. This is to align the reference frame with the spin vector of the Sun, making J2 corrections intuitive to implement.

Obliquity is typically given in the ecliptic reference frame, so snvec converts all outputs to J2000 if the os_ref_frame is equal to "HCI" and does no transformations if it is already in "J2000".

For this, it uses \(\Omega_{\odot} = 75.5940\) and \(i_{\odot} = 7.155\) as in Zeebe (2017). You can overwrite these defaults with os_omt and os_inct if desired.

ODE Solver

Note that the different ODE solver algorithm we use (Soetaert et al., 2010) means that the R routine returns an evenly-spaced time grid, whereas the C-routine has a variable time-step. This means we need to explicitly set the stepsize tres.

References

Zeebe, R.E. (2017). Numerical Solutions for the Orbital Motion of the Solar System over the Past 100 Myr: Limits and New Results. The Astronomical Journal, 154(5), doi:10.3847/1538-3881/aa8cce .

Zeebe, R. E., & Lourens, L. J. (2019). Solar System chaos and the Paleocene–Eocene boundary age constrained by geology and astronomy. Science, 365(6456), 926–929. doi:10.1126/science.aax0612 .

Zeebe, R. E., & Lourens, L. J. (2022). A deep-time dating tool for paleo-applications utilizing obliquity and precession cycles: The role of dynamical ellipticity and tidal dissipation. Paleoceanography and Paleoclimatology, e2021PA004349. doi:10.1029/2021PA004349 .

Zeebe, R. E. (2022). Reduced Variations in Earth’s and Mars’ Orbital Inclination and Earth’s Obliquity from 58 to 48 Myr ago due to Solar System Chaos. The Astronomical Journal, 164(3), doi:10.3847/1538-3881/ac80f8 .

Wikipedia page on Orbital Elements: https://en.wikipedia.org/wiki/Orbital_elements

Karline Soetaert, Thomas Petzoldt, R. Woodrow Setzer (2010). Solving Differential Equations in R: Package deSolve. Journal of Statistical Software, 33(9), 1–25. doi:10.18637/jss.v033.i09 .

Examples

# \donttest{
# \dontshow{
# set the cachedir to a temporary directory
pth <- withr::local_tempdir(pattern = "snvecR")
withr::local_options(snvecR.cachedir = pth)
# }
# default call
snvec(tend = -1e3, ed = 1, td = 0, tres = -0.4)
#> This is snvecR VERSION: 3.9.3.9000 2024-04-24
#> Richard E. Zeebe
#> Ilja J. Kocken
#> 
#> Integration parameters:
#> • `tend` = -1000 kyr
#> • `ed` = 1
#> • `td` = 0
#> • `astronomical_solution` = "full-ZB18a"
#> • `os_ref_frame` = "HCI"
#> • `os_omt` = defaulting to 75.594
#> • `os_inct` = defaulting to 7.155
#> • `tres` = -0.4 kyr
#> • `atol` = 1e-05
#> • `rtol` = 0
#> • `solver` = "vode"
#> ℹ started at "2024-04-24 21:04:56.95667"
#> Final values:
#> • s[1][2][3]: 0.404184487124565, -0.0537555129057148, and 0.913036138471423
#> • s-error = |s|-1: -5.51290422495798e-05
#> Final values:
#> • obliquity: 0.413060472710089 rad
#> • precession: -0.562357122261026 rad
#> ℹ stopped at "2024-04-24 21:05:01.232802"
#> ℹ total duration: 4.28
#> # A tibble: 2,501 × 4
#>     time   epl    phi     cp
#>    <dbl> <dbl>  <dbl>  <dbl>
#>  1   0   0.409 0      0.0163
#>  2  -0.4 0.410 0.0975 0.0168
#>  3  -0.8 0.411 0.195  0.0171
#>  4  -1.2 0.412 0.292  0.0170
#>  5  -1.6 0.413 0.389  0.0168
#>  6  -2   0.414 0.486  0.0163
#>  7  -2.4 0.414 0.582  0.0156
#>  8  -2.8 0.415 0.679  0.0146
#>  9  -3.2 0.416 0.775  0.0134
#> 10  -3.6 0.417 0.871  0.0120
#> # ℹ 2,491 more rows
# }