# Cloud albedo comparison project (George Datseris, Joaquin Blanco, Aiden Jönsson)

# Project goals / preliminary abstract
- $f$ is used by many people as a proxy for cloudiness. This is bad, because $f$ is arbitrary and without physical units. Furthermore, clouds are "densities", and besides spatial span they have "thickness" which $f$ completely misses.
- What are the different ways to quantify cloud abledo? Can we collect them and compare them?
- Discuss the difficulties in defining, and retrieving, $f$. This arises from the fact that the demarcation between cloud and non-cloud is nearly purely subjective. Since $f$ is necessarily parameterized in model physics, taking into account things such as convective organization, it is also necessary to provide observations of $f$ in order to calibrate our parameterizations. However, observing $f$ is nearly arbitrary since it is up to the scientist to decide where a cloud does and does not exist, and implicate that in an algorithm for producing such data. This puts us in a circular journey of defining and utilizing cloud fraction as a measure despite its difficulties.
- Come up with a **consistent** way to define and normalize cloud albedo so that it is energetically meaningful and would allow comparisons between different datasets and/or models.
- Each formulation of cloud albedo has been developed for a different specific purpose — it is fine to apply it to a purpose that it was meant for even if it has disadvantages, but explicitly stating why you apply that formulation should be the standard so as to avoid adding ambiguity and inconsistency across the literature to a term with specific physical meaning.
- Models, IPCC, Goverments, whatever, want accurate bounds for cloud radiative forcings. We can compare the uncertainty yielded in the shortwave part of the energy balance by using different datasets and different definitions of cloud albedo. How much W/m^2 difference does this give? If it is large, then model precision doesn't matter if our very own measurement accuracy is crap.
- Apply the "comparison pipeline" to a specific example, e.g. midlatitudes clouds, so that you have a presentable example.
- It's kind of novel to compare cloud albedo at large scales, no-one seems to wonder "why do you use this cloudiness definition".

# Tasks
## George
* Start a latex template
* Download/process ceres data into a monthly 1x1 degree format, September 2007-August 2016
* Create the code that compares _everything_.

## REST
* Give other datasets into monthly 1x1 degree format to George
* Order of importance:
* Cloud optical depth and fraction
* TOA radiation in and out
* Surface radiation fluxes
- ~~Provide summaries explanations for formulas [George, partly]~~
- Collect related literature and provide a short paragraph summary [all? whenever we find relevant literature]
- Collect datasets that provide relevant fields: Radiations (both TOA and SFC), cloud fraction, cloud optical depth [Aiden]
- Characterize/summarize key differences in the satellite data products, i.e. where they intersect, where instruments are shared etc. [Aiden]

# Terminology
Let's please follow the following terminology:

| Symbol | Description | Type |
|---------------|------------------------------------------------------|-------------------|
| $I$ | insolation | W/m² |
| $R$ | all-sky refl. insolation at TOA | W/m² |
| $K$ | clear-sky refl. insolation at TOA | W/m² |
| $Y$ | Seasonal component of $R$ | W/m² |
| $\alpha$ | all-sky albedo at TOA | fraction |
| $\alpha_K$ | clear-sky albedo at TOA | fraction |
| $C$ | Cloud albedo | fraction |
| $O$ | ice-free ocean area fraction | fraction |
| $E$ | ice / snow coverage | fraction
| $L$ | ice-free land fraction | fraction |
| $\mathcal{T}$ | time average (proper) | operation |
| $\mathcal{N}$ | northern hemisphere average | operation |
| $\mathcal{S}$ | southern hemisphere average | operation |
| $\mathcal{G}$ | global average, $\equiv (\mathcal{S}+\mathcal{N})/2$ | operation |
| $\mathcal{Z}$ | zonal average | operation |
| NH, SH | northern, southern hemisphere | abbrv. |

# 1. Cloud albedo definitions

## 1.1 Cloud Radiative effect
The simplest definition one can get is from [Ramanathan et al., 1989](https://science.sciencemag.org/content/243/4887/57): $$C_{CRE} = (R-K)/I$$

This is an effective albedo, because it can only be defined for I ≠ 0.

This is only a first-order estimate since clear-sky conditions also include the reflectivity of the atmosphere and the effects of aerosols, and what counts as clear-sky can be subjective (for instance, the hygroscopicity discussion in aerosol science). The uncertainty of radiance in clear-sky observations is a fundamental root of uncertainty in cloud albedo estimates, especially dangerous for this formulation.

## 1.2 Cloud contribution to planetary albedo
This formalism comes from the work of Donohoe2011. They provide a framework to separate planetary albedo into **additive contributions** from surface and atmosphere: $a=a^{ATM}+a^{SFC}$, based on approximating the entire atmosphere as a single layer with uniform reflectance, transmittance, and absorption. The Donohoe & Battisti decomposition can be done for both the all-sky and the clear-sky radiation fields. The precise decomposition formulas are:
$$
\alpha^\mathrm{ATM} = \frac{(\alpha - a_s t^2)}{1 - (a_s t)^2} \\ \alpha^\text{SFC} = a_s t^2(1-\alpha a_s)
$$
with $\alpha = F^\text{TOA}_\uparrow / F^\text{TOA}_\downarrow \equiv R/I$ the planetary albedo, $t = F^\text{SFC}_\downarrow / F^\text{TOA}_\downarrow$ the planetary transmittance and $a_s = F^\text{SFC}_\uparrow / F^\text{SFC}_\downarrow$ the surface albedo which in general is different than the surface *contribution* to the planetary albedo due to further reflections between surface and atmosphere. $F$ simply stands for shortwave radiation which needs to be available both at TOA and at SFC.

Once the decompositions are done, one can define the cloud contribution to planetary albedo as
$$
C_DB = a^{ATM}_R - a^{ATM}_K
$$
This is an effective albedo, because it can only be defined for $I\ne0$.

## 1.3 Cloud albedo parameterization due to Lacis-Hansen
From Lacis1974 we obtain a formula that can parameterize cloud albedo
$$
C_{LH} = f \frac{\sqrt{3}(1-g)\tau}{2+\sqrt{3}(1-g)\tau}
$$
In the paper, eq.(19) is the above formula *without* multiplication with $f$. Their eq.(19) is already a simplification of a more general formula, eq.(31). This simplification is used as the reflection of the cloud layer to calculate ozone heating. To derive this expression, the following two approximations are done: (1) the lower atmosphere is primarily a scattering region with negligible absorption (2) Rayleigh scattering is neglected. This is well justified since the size of cloud droplets causes mostly Mie scattering in the solar radiation frequency range. Furthermore, since we do not care about cloud absorption, it is probably safe to use the same formula if we only want to discuss cloud albedo, and not absorption by clouds.

It seems furthermore that their calculations of cloud albedo, eq.(31), are themselves based on formulas from Sagan and Pollack, 1967 which are a result of two-stream approximation.

## 1.4 Pincus optical thickness based approximation
Pincus2012 provides another formula that looks like this
$$
C_P = \frac{\tau^x}{c + \tau^x}
$$
which is an analytic approximation to look up tables. $x,c$ are fitted constants, and the authors find $x=0.895$ and $c=6.82$. @joaquin could you please add more info here?

## 1.5 Cloud fraction multiplier
Several papers approach the cloud albedo problem via the cloud fraction, which is freaking nonsense but anyways. [Cess (1976)](https://journals.ametsoc.org/view/journals/atsc/33/10/1520-0469_1976_033_1831_ccaaoa_2_0_co_2.xml) represents the change in planetary albedo with cloud fraction as
$$
\alpha = f a^{cloud} + (1-f)a^{clear}
$$
Then one can re-write this to derive
$$
C_f = a^{cloud} = \frac{1}{f}\left(\alpha - (1-f)a^{clear}\right)
$$

This assumes that $\alpha_{cloud}$ is inhomogeneous and can vary globally and with time, and is used as a tool to characterize the effects of clouds on TOA radiation by clouds by e.g. typology or region. Since it is only dependent on a subjective measure (fraction), it misses valuable information about the determinants of the effective albedo of clouds.

This is best used as a diagnostic: when you compare the dependence of the albedo impact of a cloud on planetary albedo with its cloud fraction, it should tell you how much its microphysical properties dominate the albedo impact. This can be used to compare clouds in different areas or across data sets to see how the "bulk" properties and their impacts compare. However, this is dependent on the *fraction being correct;* the ambiguity of $f$ thus makes it difficult for this diagnostic to be used. This illustrates the difficulties encountered when using $f$ to give physically meaningful information.

## 1.6 Energetic normalization
The parameterization based albedos $C_P$ and $C_{LH}$ have the problem that they are not energetically consistent, exactly because they do not use radiation to estimate their values. Datseris2021 established a process that tunes a "free" parameter of the parameterizations based on the energetically consistend $C_{DB}$. For $C_P$ this parameter is arguably $c$, for $C_{LH}$ it is $g$. The process requires $c$ or $g$ to be distributed spatially so that the temporally averaged maps of $C_{LH}$ (or $C_P$) agree with $C_{DB}$.

# Satellite data products
- Q: should we focus on monthly, 1° gridded data?
- Q: Overlap period: should we separate into sections where different CDRs overlap, or specifically focus on periods of overlap between as many CDRs as possible? (e.g. 2000-2009 should provide overlap between all of the above data sets)
- Q: How much consideration for satellite overlap should we include? I.e. would we like to diversify to minimize satellite measurement repeats among data sets? It could be worth it to disregard, since the focus here is on physical representations of cloud properties among data sets, not irradiance measurements
- To add: quick runthrough on how radiative heating fluxes are obtained, because this changes how we think/talk about clouds. This to highlight the difference between modeled, idealized SW heating fluxes and observed fluxes

After looking through the data sets, what I believe would suit us are:
- CERES SYN1deg
- MODIS Terra
- MODIS Aqua
- ISCCP-FH (MPF for monthly, FH-TOA and FH-SRF for daily)
- ~~CLARA A2 NOAA~~
- ~~CLARA A2 METOP~~

(CLARA-A2 does not have all of the necessary fields, sadly. CLARA-A3 will; but the release has been pushed back to 2022.)

All for the period September 2007 - August 2016, giving 10 years of data and plenty of diversity among satellites. (New suggested dates: 2002-2017?)

## ISCCP
Coverage: July 1983 - June 2009 (2017) for the D/B (H) series. Global coverage provided by 5 geostationary satellites (GOES-East, GOES-West, GMS, INSAT, METEOSAT) complemented by at least one polar orbiting NOAA satellite (AVHRR Global Area Coverage). The data is primarily from two standard visible (0.6 μm) and IR (11 μm) channels common to (roughly) all the satellites.

ISCCP-FH is the ISCCP flux data set containing retrieved surface and TOA fluxes, which is provided in 3-hourly time resolution and in monthly means on a 1° grid. This would be a limiting factor time-wise, providing an end date for the comparison (June 2017).

### Disentangling ISCCP-FH
- There is a "short" data type that NASA uses for some reason, which uses a conversion factor and a formula to obtain the "real" (actual) value:
`txu4cl:type = "short: convert to Real*4 = Short_value / scale_factor"`
- `REAL*4` is the Fortran way of denoting a 4-byte (32-bit) float; therefore it is just showing that the float value = short value/scale factor.
- After loading it with `xarray`, I found this note: `Note: Physical value = float(short value)/scale_factor, excluding missing_value`
- All conversion factors are found in this file: https://isccp.giss.nasa.gov/pub/flux-fh/docs/List_ISCCP-FH_product.txt
- An important note: solar incoming radiation in ISCCP is 1367 Wm^-2, while in CERES it is 1361.

To be done before this data is useable:
- Export only the variables we need from FH-MPF, which are:
- All-sky surface fluxes (up: `sxu1fl`, down: `sxd1fl`) (conversion factor: 10)
- Clear-sky surface fluxes (up: `sru1cr`, down: `srd1cr`) (conversion factor: 10)
- TOA reflected radiation (all-sky: `sxu5fl`, clear-sky: `sru5cr`) (conversion factor: 10)
- Downwelling shortwave radiation at TOA (`sxd5fl`) (conversion factor: 10)

The above step was all done with:
`cdo select,name=sxu1fl,sxd1fl,sru1cr,srd1cr,sxu5fl,sru5cr,sxd5fl ISCCP-FH.MPF.v.0.0.GLOBAL.200701-201706.nc out1.nc`

- Export only the variables we need from the monthly means of FHD:
- Cloud fraction (`cf_m__`) (conversion factor: 1000)
- Cloud optical thickness (`tau_m_`) (conversion factor: 10)

The above step was done with:
`cdo select,name=cf_m__,tau_m_ FH/ISCCP-FH.TOA.v.0.0.GLOBAL.200701-201706.nc out2.nc`

- Combine them into one data file with `cdo merge out1.nc out2.nc out.nc`
- Ensure a time dimension is included with a proper time vector; the dimension is there after using `ncecat -u time`, but needs to be filled. This is now done with: `cdo -settaxis,2007-01-15,00:00,1month out.nc out.nc`; this time axis may need to be adjusted, though.
- Convert everything to the "real" data values, so it is ready to be used elsewhere when loaded

## CERES
Coverage: March 2000 - present
Instrument: CERES instrument aboard Terra (December 1999) and Aqua (May 2002) measures in SW (0.3-5 μm) and total (0.3-200 μm) bands, as well as a window band (8.11-11.8 μm) where water vapor does not influence OLR much, used to investigate the role of water vapor in climate ([Smith et al. 2002](https://www.sciencedirect.com/science/article/pii/S0273117702802800)).

For hourly and daily data, there is only CERES SYN; for our purposes, CERES SYN1deg would probably suit us best (hourly, 3-hourly and daily temporal resolution, 1° grid boxes). There is no balanced (EBAF/Energy Balanced And Filled) CERES data at resolutions higher than monthly.

Cloud properties are based on MODIS, VIIRS and GEO. Since there is data shared between MODIS and CERES, it would be interesting/helpful to note how cloud albedo may differ in these two data sets, since the cloud albedo has a global impact on radiation balance and therefore could lead to large discrepancies in analyses of cloud radiative effects.

## MODIS
The time coverage is from December 1999 to present for Terra, and from May 2002 to present for Aqua. Instrument: 36 bands spanning 0.4-14.385 μm. Data is 5-minutely and at 1 & 5 km resolution. Cloud properties are only available for Terra and Aqua individually, and not for combined L3 products.

Downloads here: https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/science-domain/cloud

## CLARA
CLARA A2 ([v2.1](https://wui.cmsaf.eu/safira/action/viewProduktSearch?menuName=PRODUKT_SUCHE)) is provided by CM SAF at EUMETSAT, and is based on AVHRR sensors aboard NOAA and METOP (operated by EUMETSAT) satellites. The temporal resolution is monthly and daily, and global coverage is given on 0.25° grid cells. However, METOP satellites were only operational from September 2007 on, and so between 1982-2007 the CDR is based only on NOAA polar-orbiting satellites.

# Summaries of relevant papers

## Lacis174
Lacis & Hansen devised parameterizations of several absorption and reflection processes of solar radiation for use in general circulation models. To derive the parametizations they did detailed multiple scattering calculations using analytic formulas integrated on a computer. They fitted parameterizations on the multiple scattering numeric results. The variables of the parameterization are water vapor, amount of clouds, zenith angle, surface albedo and ozone distributions.

The multiple scattering method is based on a plane-parallel aproximation of atmosphere (i.e. only vertical component exists). The atmosphere is split up into a sufficient number of layers (~50), each one assummed homogenous. The monochromatic scattering properties of each layer are determined completely by its optical depth $\tau$, the single scattering albedo $\tilde{\omega}_0$ and the phase function $p(a)$ with $a$ the scattering angle. The surface is considered as just one more layer with 0 transmission.

# References

Cess, R. D. (1976). Climate Change: An Appraisal of Atmospheric Feedback Mechanisms Employing Zonal Climatology, Journal of Atmospheric Sciences, 33(10), 1831-1843. https://journals.ametsoc.org/view/journals/atsc/33/10/1520-0469_1976_033_1831_ccaaoa_2_0_co_2.xml

Ramanathan, V. L. R. D., Cess, R. D., Harrison, E. F., Minnis, P., Barkstrom, B. R., Ahmad, E., & Hartmann, D. (1989). Cloud-radiative forcing and climate: Results from the Earth Radiation Budget Experiment. Science, 243(4887), 57-63. https://science.sciencemag.org/content/243/4887/57

Donohoe2011: https://doi.org/10.1175/2011JCLI3946.1