R: Calculate Earth’s Obliquity and Precession in the Past

snvec {snvecR}

R Documentation

Calculate Earth’s Obliquity and Precession in the Past

Description

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.

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).

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().

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



# default call
snvec(tend = -1e3, ed = 1, td = 0, tres = -0.4)

[Package snvecR version 3.9.4 Index]