Unexpected bumps in insolation (daily_insolation function) observed at high eccentricities (e > ~0.4)

[sid-bhatnagar]

It seems that the daily_insolation function works as expected for 2 kinds of external forcing ('long_peri' and 'obliquity'), but not for the whole range of 'ecc'. Small eccentricities (up to e ~ 0.4) work as expected, but for higher eccentricities, unexpected bumps in the daily average insolation occur approximately halfway between the perihelion and aphelion (for e.g., see contour plot (output_3_0.png) in the attached zip file).
For a slightly more detailed description of the problem, please refer to the attached short version of the Jupyter notebook (zip containing .md and output files).
Insolation_to_send_to_climlab.zip

[brian-rose]

Thanks for this feedback!

It may be a while before I have a chance to look at this carefully. Do you have a sense for whether there's a numerical issue, or is it rather that the Berger expansion formulas become inaccurate for e > 0.4?

(My first guess would be the second option)

[sid-bhatnagar]

Yes, my starting guess is similar, that the Berger expansion (in the MATLAB code, as per the documentation) is inaccurate for high eccentricities. The equation which calculates the insolation here, Fsw, has a term that is not in the Berger (1978) equation (eq. 10) - essentially, the constant (86.4/rho**2) in the 1978 paper becomes the term (1+ecc.*cos(lambda-omega)).^2 / (1-ecc.^2).^2 here.

Since this is the only difference in the calculation of insolation between Berger (1978) and the code, this could be a good starting point for the investigation.

P.S: Just as an initial check, I plotted the above term for the range of eccentricities keeping lambda = omega (planet at perihelion). The corresponding plot doesn't seem to have any issues - it is a purely increasing non-linear function. Associated code:

lambdaa = 281.37 #True longitude of Earth at perihelion (measured from vernal equinox)
omega = 281.37 #LOP (measured from vernal equinox)

ecc = np.arange(0.0,0.99,0.01)

constant = [( ((1+(i*np.cos(lambdaa-omega)))**2) / ((1-(i)**2)**2) ) for i in ecc]
plt.plot(ecc,constant)

[sid-bhatnagar]

Hello! Just checking in to see if this issue has been rectified. Please let me know when you can.

Just to confirm, I am using the latest version of the package, 0.7.12.

Thanks!

[brian-rose]

Hi, no, I have not had any time to look at this. Your help is welcome on a revision!

[sid-bhatnagar]

Okay, thanks for the reply. Will attempt a revision!

[brian-rose]

Excellent, thanks!

FYI the formulas are implemented in climlab in the source file insolation.py

[sid-bhatnagar]

Thanks for the file! Will have a look.