Exponent of effective_tau_cirrus

[Wessel99]

Description

Hello pycontrail team,

Looking at the 'effective_tau_cirrus' function in radiative_forcing.py, it seems that in the code the product (tau_cirrus_eff * delta_sc_aps) is outside the exponent, while the source paper places it inside the exponent. Which one is correct?

Best regards,
Wessel

Details

• Version: 0.47.2
• Module: radiative_forcing.py
• OS: windows

Additional Notes

def effective_tau_cirrus(
tau_cirrus: npt.NDArray[np.float_],
mue: npt.NDArray[np.float_],
delta_sc: npt.NDArray[np.float_],
delta_sc_aps: npt.NDArray[np.float_],
) -> npt.NDArray[np.float_]:
r"""
Calculate the effective optical depth of natural cirrus above the contrail, e_sw.

Refer to Eq. (11) of Schumann et al. (2012).

Parameters
----------
tau_cirrus : npt.NDArray[np.float_]
    Optical depth of numerical weather prediction (NWP) cirrus above the
    contrail for each waypoint.
mue : npt.NDArray[np.float_]
    Cosine of the solar zenith angle (theta), mue = cos(theta) = sdr/sd0
delta_sc : npt.NDArray[np.float_]
    Habit-specific parameter to account for the optical depth of natural
    cirrus above the contrail
delta_sc_aps : npt.NDArray[np.float_]
    Habit-specific parameter to account for the optical depth of natural
    cirrus above the contrail

Returns
-------
npt.NDArray[np.float_]
    Effective optical depth of natural cirrus above the contrail,
    ``n_waypoints x 8`` (habit) columns.
"""
tau_cirrus_eff = tau_cirrus / (mue + 1e-6)
return np.exp(tau_cirrus * delta_sc) - (tau_cirrus_eff * delta_sc_aps)

[zebengberg]

Thanks @Wessel99 ! I agree that the code differs from Schumann's paper, and that this most likely a bug. @roger-teoh, can you confirm this?

[roger-teoh]

Hi @Wessel99. Thank you so much for checking our codes in detail. We really appreciate this as you successfully identified another bug in the codes. Upon further investigation, both CoCiP Fortran and pycontrails are affected by a bug when calculating E_SW, c.f. Eq. (11) in Schumann et al. (2012), and I have just alerted Professor Schumann about this.

@zebengberg Yes, all the terms should be inside the exponent. The following is the correct version:

return np.exp((tau_cirrus * delta_sc) - (tau_cirrus_eff * delta_sc_aps))

[Wessel99]

Glad to be of help! Thank you for clarifying which version is correct.

[roger-teoh]

Hi @Wessel99 , we just had a follow up discussion with Professor Schumann and identified another bug in the code and the paper.

Professor Schumann identified a print error in Eq. (11) of Schumann et al. (2012), where the positions of delta_sc_aps and delta_sc should be swapped. The correct equation is presented below:

exp(deltasc’* tauc – deltasc * tauc_eff)

We are working to incorporate this change into pycontrails, and only expect very minor changes in the computed contrail RF < 0.1%.

[zebengberg]

Fixed in #100