An interactive ocean surface albedo scheme (OSAv1.0): formulation and evaluation in ARPEGE-Climat (V6.1) and LMDZ (V5A)

Ocean surface represents roughly 70 % of the Earth's surface, playing a large role in the partitioning of the energy flow within the climate system. The ocean surface albedo (OSA) is an important parameter in this partitioning because it governs the amount of energy penetrating into the ocean or reflected towards space. The old OSA schemes in the ARPEGE-Climat and LMDZ models only resolve the latitudinal dependence in an ad hoc way without an accurate representation of the solar zenith angle dependence. Here, we propose a new interactive OSA scheme suited for Earth system models, which enables coupling between Earth system model components like surface ocean waves and marine biogeochemistry. This scheme resolves spectrally the various contributions of the surface for direct and diffuse solar radiation. The implementation of this scheme in two Earth system models leads to substantial improvements in simulated OSA. At the local scale, models using the interactive OSA scheme better replicate the day-to-day distribution of OSA derived from ground-based observations in contrast to old schemes. At global scale, the improved representation of OSA for diffuse radiation reduces model biases by up to 80 % over the tropical oceans, reducing annual-mean model–data error in surface upwelling shortwave radiation by up to 7 W m⁻² over this domain. The spatial correlation coefficient between modeled and observed OSA at monthly resolution has been increased from 0.1 to 0.8. Despite its complexity, this interactive OSA scheme is computationally efficient for enabling precise OSA calculation without penalizing the elapsed model time.

1 Introduction

The surface radiation budget has long been recognized as fundamental to our understanding of the climate system (IPCC, 2001, 2007, 2013). The flow of radiative energy through the Earth system and the radiative interactions between the atmosphere and the ocean remain one of the major sources of uncertainty in climate predictions (Allen et al., 2009; Frölicher, 2016; Gillett et al., 2013; Myhre et al., 2015; Otto et al., 2013).

In the atmosphere, the spatiotemporal variations in incoming solar radiation and its atmospheric absorption drive the hydrological cycle as well as the flow of air masses. In the oceans, the fraction of solar radiation penetrating the subsurface is controlled by the ocean surface albedo (OSA). The corresponding amount of heat stored into the ocean constitutes an important term in the ocean energy surface balance and affects in turn the whole climate system. On short (daily to seasonal) timescales, solar radiation absorbed into the upper-ocean layers affects the stability of the ocean mixed layer, the sea surface temperature and may, in turn, influence the geographic structure of large-scale atmospheric convection (Gupta et al., 1999). Over longer timescales, the fraction of energy entering the ocean contributes to increase the ocean heat content, which is a key term for determining climate sensitivity from observations (Otto et al., 2013).

OSA interacts with a multitude of biophysical processes occurring in the first meters of the ocean. In particular, it governs the amount of solar radiation entering the upper-most layer of the ocean which interacts with marine biological light-sensitive pigment like chlorophyll and other materials in suspension (e.g., Morel and Antoine, 1994; Murtugudde et al., 2002). OSA also influences a number of biogeochemical processes such as photosynthesis or photolysis, which respond to incoming solar radiation within the upper-most layer of the ocean. Conversely, penetrating ultraviolet solar radiation can also have detrimental impacts on the marine biota (e.g., Li et al., 2015; Smyth, 2011). Consequently, OSA influences marine primary productivity directly, and hence also ocean ecosystems and ocean carbon uptake (e.g., Behrenfeld and Falkowski, 1997; Nelson and Smith, 1991; Siegel et al., 2002).

Despite its importance, OSA is a parameter that receives insufficient attention from both an observational and modeling point of view. Most of the available data are indirectly retrieved from satellite observations of the top-of-atmosphere radiative budget (Wielicki et al., 1996), with relatively few direct observations of surface radiative fluxes. Nonetheless, the OSA processes are relatively well understood so OSA can be parameterized at the global scale.

Both empirical and theoretical approaches indicate that the solar zenith angle (SZA) is the single most prominent driving parameter for OSA. However a wide range of other parameters such as the partitioning of incoming solar radiation between its direct and diffuse components, the sea surface state (often approximated through the surface wind), the concentration of suspended matter and plankton light-sensitive pigment in the surface ocean and the extent and physical properties of whitecaps also affect OSA. All of these contributions vary spectrally and OSA thus depends on the spectral distribution of the incoming solar radiation at the surface (Jin et al., 2002, 2004).

Over recent decades, several schemes have been proposed to model OSA (e.g., Cox and Munk, 1954; Hansen et al., 1983; Kent et al., 1996; Larsen and Barkstrom, 1977; Preisendorfer and Mobley, 1986; Ohlman et al., 2000). Some schemes depend only on the solar zenith angle while others additionally depend on quantities like wind speed or cloud optical depth, inducing substantial differences in OSA patterns and variability. Li et al. (2006) investigated the impact of various OSA schemes in the Canadian atmosphere climate model, AGCM4. The authors show that the difference in clear-sky upwelling shortwave radiation between schemes can reach 20 W m⁻² at the top of the atmosphere and more than 20 W m⁻² at the surface.

Most of the schemes assessed in Li et al. (2006) do not resolve spectral variations in OSA, thus excluding the possibility of representing subtle processes and couplings in Earth system models as suggested by complex ocean radiative transfer models (e.g., Ohlman et al., 2000; Ohlman and Siegel, 2000). Indeed, changes in whitecaps and ocean color, whether due to climate variability or climate change, can modify the OSA, with potential impacts on photochemistry in the atmosphere and biological activity in the upper-most layer of the ocean (Hense et al., 2017).

In this study, we propose a new interactive OSA scheme well adapted to the current generation of Earth system models which may benefit from and provide benefits to the coupling between Earth system model components like surface ocean waves or marine biogeochemistry. This study provides details of its implementation into two atmospheric models and discusses its performance on daily to seasonal timescales.

The outline is as follows. In Sect. 2, we introduce the formulation of the interactive OSA scheme which is derived from old schemes published in literature over recent decades. In Sect. 3, we analyze the importance of the various components of this scheme using a stand-alone version of the OSA scheme. Section 4 describes the experimental design and the two state-of-the-art atmospheric models that are used in Sects. 5, 6 and 7 to evaluate the interactive OSA scheme against available observations. Section 8 concludes the present study.

2 Interactive ocean surface albedo parameterization

Albedo is a measure of the reflectivity of a surface and is defined as the fraction of the incident solar radiation that is reflected by the surface. It depends not only on the properties of the surface but also on the properties of the solar radiation incident on that surface. Technically speaking albedo can be computed from the knowledge of the spectral and directional distribution of the incident solar radiation L (λ, θ, ϕ) and the bidirectional reflectance distribution function (BRDF), ρ (λ, θ_i, ϕ_i, θ_r, ϕ_r), which links the reflected radiation of direction (θ_r, ϕ_r) to that of incident radiation of direction (θ_i, ϕ_i). Here, λ represents the wavelength, and (θ, ϕ) the zenith and azimuthal angles.

While the atmospheric incident radiation, L (λ, θ, ϕ), can be solved using a radiative transfer model and the BRDF can be modeled from the knowledge of the surface ocean properties, the complexity and the computational cost of such models are prohibitive for climate applications. Thus, estimation of OSA in climate models has to rely on several simplifying assumptions. In particular, incident solar radiation is usually characterized by a downward direct flux (for which SZA is known) and a diffuse downward flux (for which a typical angular distribution can be assumed).

In the present work, most of the analytic formulations employed are derived from the azimuthally averaged radiative transfer equation (Chandrasekhar, 1960), enabling a straightforward estimation of the OSA for direct and diffuse radiation. This implies that zenith solar angle is the only directional parameter involved in the parametrization.

The suite of processes involved in our scheme is displayed in Fig. 1. The incident solar radiation (either direct and diffuse) is first influenced by the presence of foam (composing the whitecaps), which exhibits different reflective properties from seawater. Then, the reflective properties of the uncapped fraction of the sea surface are determined separately for direct and diffuse incident radiation. Finally, the subsurface – or ocean interior – reflectance of seawater is computed for both direct and diffuse incident radiation.

Hereafter, we describe the various components of OSA according to the nature of the incident solar radiation (direct or diffuse) and the processes involved in its reflection.

2.1 Treatment of whitecaps

The first contribution to the new interactive OSA scheme is the whitecap cover. Indeed, the whitecap albedo, α^WC(λ), could significantly increase the OSA at high wind speeds (e.g., Frouin et al., 2001, 1996; Gordon and Wang, 1994; Stramska, 2003).

Figure 1 Pathways of solar radiation over oceans as described in the new interactive scheme. Whitecaps, surface Fresnel, ocean interior influence of the reflection and the refraction of both direct and diffuse radiation.

The whitecaps originate from turbulence induced by the breaking of waves, which generates foam at the sea surface (Deane and Stokes, 2002; Melville and Matusov, 2002). In the absence of an ocean wave model such as WAM (e.g., Aouf et al., 2006; Ardhuin et al., 2010) which would provide more accurate whitecap coverage (WC) based on significant wave height (Bell et al., 2013; Woolf, 2005), we used the formulation of WC published in Salisbury et al. (2014). Their expression is based on recent spaceborne observations with a 37 GHz channel radar. It parametrizes WC as a function of the 10 m wind speed, w, in unit of m s⁻¹:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{WC}\left(w\right)=\mathrm{3.97}\times {\mathrm{10}}^{-\mathrm{2}}{w}^{\mathrm{1.59}}.\end{array}$$

As mentioned in Salisbury et al. (2014), this approximation of WC is valid for w ranging between 2 and 20 m s⁻¹, which corresponds well to the range of w values simulated by the current generation of Earth system models. The formulation employed here does not account for the temperature dependence of wave breaking in agreement with other parameterizations for WC (Stramska, 2003) because its effect is weaker than that of surface winds on the whitecap coverage.

In order to solve the spectral dependence of the whitecap albedo, we use the relationship proposed by Whitlock et al. (1982). However, we rely on previous work indicating that the whitecap albedo of ordinary foam, α^WC(λ), tends to be twice as small as that of fresh and dense foam (Koepke, 1984). We consequently apply a 1∕2 coefficient to the formulation proposed by Whitlock et al. (1982) for α_WC(λ) from 400 to 2400 nm, as follows:

$$\begin{array}{ll}{\displaystyle }{\mathit{\alpha }}^{\mathrm{WC}}\left(\mathit{\lambda }\right)=& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\times {\displaystyle \frac{\mathrm{1}}{\mathrm{100}}}(\mathrm{60.063}-\mathrm{5.127}\ln r_{\mathrm{WC}}(\mathit{\lambda })\\ {\displaystyle }& {\displaystyle }+\mathrm{2.779}{\left(\ln r_{\mathrm{WC}}\left(\mathit{\lambda }\right)\right)}^{\mathrm{2}}-\mathrm{0.713}{\left(\ln r_{\mathrm{WC}}\left(\mathit{\lambda }\right)\right)}^{\mathrm{3}}\\ \text{(2)}& {\displaystyle }& {\displaystyle }+\mathrm{0.044}{\left(\ln r_{\mathrm{WC}}\left(\mathit{\lambda }\right)\right)}^{\mathrm{4}}),\end{array}$$

where r_WC is the absorption coefficient of clear water in m⁻¹. We use r_WC(λ) as published in Whitlock et al. (1982) from 400 to 2400 nm. Outside the 400 to 2400 nm range, we chose to set r_WC(λ) to zero due to the lack of available data in the literature. Tabulated values of α^WC(λ) computed using this assumption are provided in Table S1 in the Supplement.

In the present scheme, we assume that whitecaps reflect equally direct and diffuse incoming solar radiation. We also assume that our formulation, which is based on observations, is not affected by the small contribution from the subsurface reflectance.

2.2 Treatment of the uncapped surface

2.2.1 Fresnel surface albedo for direct radiation

We describe the contribution of Fresnel reflection at the ocean surface, which is a major component of the OSA. Fresnel reflection is assumed to depend at a given wavelength (λ) on the solar zenith angle (θ), the refractive index of seawater (n) and the two-dimensional distribution of the ocean surface slopes (f). As mentioned earlier we neglect the dependence on the azimuthal angle of the incident radiation.

We follow Jin et al. (2011) and express the direct surface albedo (${\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}})$ as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}\left(\mathit{\lambda },\mathit{\theta },w\right)=r_{\mathrm{f}}\left(n\left(\mathit{\lambda }\right),\mathit{\mu }\right)-{\displaystyle \frac{r_{\mathrm{f}}\left(n\left(\mathit{\lambda }\right),\mathit{\mu }\right)}{r_{\mathrm{f}}\left(n_{\mathrm{0}},\mathit{\mu }\right)}}f\left(\mathit{\mu },\mathit{\sigma }\right),\end{array}$$

where μ=cos (θ), n is the spectral refractive index of seawater, r_f is the Fresnel reflectance for a flat surface and f is a function that accounts for the distribution of multiple reflective facets at the ocean surface. Tabulated values for n(λ) are shown in Table S1. n₀=1.34 corresponds to the refractive index of seawater averaged from 300 to 700 nm (i.e., visible spectrum) for which the f function is estimated.

The interaction of incident shortwave radiation with the multiple reflective facets of the ocean surface at various angles and directions is difficult to model. It is nonetheless possible to represent statistically the distribution of the slope of reflective ocean facets with a probabilistic function. The probabilistic function provided by Cox and Munk (1954) assumes a Gaussian distribution of mean slope facet as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }p\left(\tan \mathit{\vartheta }\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}}}\exp \left({\displaystyle \frac{-{\tan }^{\mathrm{2}}\mathit{\vartheta }}{{\mathit{\sigma }}^{\mathrm{2}}}}\right),\end{array}$$

where ϑ is the facet angle, i.e., the angle between the normal to the facet and the normal to the horizontal ocean surface, and σ is the width the distribution of the facet angle. The parameter σ, also called the surface roughness, is modulated by the influence of surface (i.e., 10 m) wind speed (w) as ${\mathit{\sigma }}^{\mathrm{2}}=\mathrm{0.003}+\mathrm{0.00512}w$. The formulation by Cox and Munk (1954) assume that (1) shading influence of ocean facets is neglected and (2) ocean surfaces never behave as a theoretical Fresnel surface (requiring σ=0). These approximations can impact ${\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}$ calculation at high SZA and/or in absence of wind. In any case, this formulation (based on wind speed only) ignores the effect of the wind direction on the wind sea and the effect of swell, which both affect the distribution of slopes. This latter set of assumptions can also be revised in the foreseeable future when climate models include an interactive ocean wave model.

In order to account for various impacts of multiple ocean surface facets on ${\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}$ including both multiple scattering (increasing surface reflection) and shading effect (reducing reflection), Jin et al. (2011) proposed expressing f as a polynomial function. This function intends to parameterize the mean contribution of multiple reflective facets at the ocean surface to ${\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}$ using only the parameters μ and σ. This polynomial function is expressed as follows:

$$\begin{array}{ll}{\displaystyle }f\left(\mathit{\mu },\mathit{\sigma }\right)=& {\displaystyle }(\mathrm{0.0152}-\mathrm{1.7873}\mathit{\mu }+\mathrm{6.8972}{\mathit{\mu }}^{\mathrm{2}}-\mathrm{8.5778}{\mathit{\mu }}^{\mathrm{3}}\\ {\displaystyle }& {\displaystyle }+\mathrm{4.071}\mathit{\sigma }-\mathrm{7.6446}\mathit{\mu }\mathit{\sigma })\times \exp (\mathrm{0.1643}\\ {\displaystyle }& {\displaystyle }-\mathrm{7.8409}\mathit{\mu }-\mathrm{3.5639}{\mathit{\mu }}^{\mathrm{2}}-\mathrm{2.3588}\mathit{\sigma }\\ \text{(5)}& {\displaystyle }& {\displaystyle }+\mathrm{10.054}\mathit{\mu }\mathit{\sigma }).\end{array}$$

Coefficients of f have been fitted using several accurate calculations of ${\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}$ using a radiative transfer model (Jin et al., 2006, 2005).

2.2.2 Fresnel surface albedo for diffuse radiation

The amount and distribution of incident diffuse radiation strongly depend on the amount and characteristics of cloud and aerosols. It is therefore difficult to derive an analytical formulation for ${\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{S}}$ from a BRDF that would be applicable to all atmospheric conditions. We therefore choose to use the simple expression for the diffuse surface albedo (${\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{S}})$ under cloudy sky proposed in Jin et al. (2011):

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{S}}\left(\mathit{\lambda },w\right)=-\mathrm{0.1479}+\mathrm{0.1502}n\left(\mathit{\lambda }\right)-\mathrm{0.0176}n\left(\mathit{\lambda }\right)\mathit{\sigma },\end{array}$$

with σ and n defined as previously.

2.3 Contribution of the ocean interior reflectance to surface albedo

In this section, we describe the contribution of the ocean interior reflectance to the ocean surface albedo, α^W. It is caused by solar radiation penetrating the ocean but eventually returning to the atmosphere after one or multiple reflections within the seawater volume. Below the ocean surface, solar radiation interacts not only with seawater but also with material in suspension in the water like the marine biological pigment or detrital organic materials (DOM). Previous studies show that DOM can influence radiative properties of the open ocean (e.g., Behrenfeld and Falkowski, 1997; Dutkiewicz et al., 2015; Kim et al., 2015). However, we chose to solely account for the influence of the marine biological pigment which is characterized by its chlorophyll content because the influence of DOM on ocean surface albedo is expected to be small compared to the surface chlorophyll. Furthermore the abundance of chlorophyll in seawater has been monitored from space for decades (e.g., Behrenfeld et al., 2001; Siegel et al., 2002; Yoder et al., 1993) by ocean color measurements. Such observations can provide a climatology to use in the climate model in absence of an ocean biogeochemical module.

Over the uncapped ocean surface, the fraction of direct radiation penetrating into the upper-most layer of the ocean $\mathrm{1}-{\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}$ interacts with the seawater, which has reflectance, R₀. Upwelling radiation can be reflected downward at the air–sea interface with reflectance, r_w. Therefore, the contribution of multiple reflections of penetrating radiation to the ocean albedo takes the form of the following Taylor series:

$$\begin{array}{ll}{\displaystyle }{\mathit{\alpha }}^{\mathrm{W}}\left(\mathit{\lambda },\mathit{\theta },\mathrm{Chl}\right)=& {\displaystyle }\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}\right)\left(\mathrm{1}-r_{\mathrm{w}}\right)R_{\mathrm{0}}(\mathrm{1}+r_{\mathrm{w}}{R}_{\mathrm{0}}\\ \text{(7)}& {\displaystyle }& {\displaystyle }+(r_{\mathrm{w}}{R}_{\mathrm{0}}{)}^{\mathrm{2}}+(r_{\mathrm{w}}{R}_{\mathrm{0}}{)}^{\mathrm{3}}+\mathrm{\cdots }),\end{array}$$

which can be arranged as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathit{\alpha }}^{\mathrm{W}}\left(\mathit{\lambda },\mathit{\theta },\mathrm{Chl}\right)={\displaystyle \frac{\left(\mathrm{1}-r_{\mathrm{w}}\right)R_{\mathrm{0}}}{\mathrm{1}-r_{\mathrm{w}}{R}_{\mathrm{0}}}}\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}\right).\end{array}$$

We employ the formulation of r_w and R₀ proposed by Morel and Gentili (1991).

These authors express r_w as a function of surface roughness σ, that is,

$$\begin{array}{}\text{(9)}& {\displaystyle }{r}_{\mathrm{w}}=\mathrm{0.4817}-\mathrm{0.0149}\mathit{\sigma }-\mathrm{0.2070}{\mathit{\sigma }}^{\mathrm{2}}\end{array}$$

R₀ represents an apparent optical property of seawater, which can be written as follows:

$$\begin{array}{}\text{(10)}& {\displaystyle }{R}_{\mathrm{0}}\left(\mathit{\lambda }\mathit{\eta }\mathit{\mu }\mathrm{Chl}\right)=\mathit{\beta }\left(\mathit{\eta },\mathit{\mu }\right){\displaystyle \frac{\mathrm{0.5}{b}_{\mathrm{w}}\left(\mathit{\lambda }\right)+b_{\mathrm{bp}}(\mathit{\lambda },\mathrm{Chl})}{a_{\mathrm{w}}\left(\mathit{\lambda }\right)+a_{\mathrm{bp}}(\mathit{\lambda },\mathrm{Chl})}},\end{array}$$

where a_w(λ) and b_w(λ) are the absorption and backscattering coefficients of seawater (in m⁻¹); a_bp(λChl) and b_bp(λChl) are absorption and backscattering coefficients of biological pigments (i.e., chlorophyll).

β(η,μ) is function of seawater and biological pigment backscattering and can be written as

$$\begin{array}{ll}{\displaystyle }\mathit{\beta }\left(\mathit{\eta },\mathit{\mu }\right)=& {\displaystyle }\mathrm{0.6270}-\mathrm{0.2227}\mathit{\eta }-\mathrm{0.0513}{\mathit{\eta }}^{\mathrm{2}}\\ \text{(11)}& {\displaystyle }& {\displaystyle }+(\mathrm{0.2465}\mathit{\eta }-\mathrm{0.3119})\mathit{\mu },\end{array}$$

where μ=cos (θ) and $\mathit{\eta }=\frac{\mathrm{0.5}{b}_{\mathrm{w}}\left(\mathit{\lambda }\right)}{\mathrm{0.5}{b}_{\mathrm{w}}\left(\mathit{\lambda }\right)+b_{\mathrm{bp}}(\mathit{\lambda },\mathrm{Chl})}$.

Backscattering of biological pigment, b_bp, is computed using the formulation proposed in Morel and Maritorena (2001) which uses chlorophyll concentration, [Chl], as a surrogate of biological pigment concentration as follows:

$$\begin{array}{ll}{\displaystyle }{b}_{\mathrm{bp}}\left(\mathit{\lambda }\right)=& {\displaystyle }\mathrm{0.416}\left[\mathrm{Chl}{]}^{\mathrm{0.766}}\right(\mathrm{0.002}+{\displaystyle \frac{\mathrm{1}}{\mathrm{100}}}\left(\mathrm{0.50}\right.\\ \text{(12)}& {\displaystyle }& {\displaystyle }\left.-\mathrm{0.25}\ln \left[\mathrm{Chl}\right]\right){\left({\displaystyle \frac{\mathit{\lambda }}{\mathrm{550}}}\right)}^{\mathrm{0.5}\left(\ln \left[\mathrm{Chl}\right]-\mathrm{0.3}\right)}),\end{array}$$

with λ expressed here in nm and [Chl] in mg m⁻³. This formulation is valid for [Chl] ranging between 0.02 and 2 mg m⁻³.

Figure 2 Comparison of chlorophyll absorbance (m⁻¹) as a function of wavelength (nanometers, in log scale) from Morel and Maritorena (2001) in blue and that of Frigaard et al. (1996) in red, which is used in the new interactive OSA scheme.

The absorption of biological pigment, a_bp(λChl), is also computed using Morel and Maritorena (2001) formalism:

$$\begin{array}{ll}{\displaystyle }{a}_{\mathrm{bp}}\left(\mathit{\lambda }\right)=& {\displaystyle }\mathrm{0.06}{a}_{\mathrm{chl}}\left(\mathit{\lambda }\right)[\mathrm{Chl}{]}^{\mathrm{0.65}}+\mathrm{0.2}(\mathrm{0.00635}\\ \text{(13)}& {\displaystyle }& {\displaystyle }+\mathrm{0.06}\left[\mathrm{Chl}{]}^{\mathrm{0.65}}\right)e^{\mathrm{0.014}\times (\mathrm{440}-\mathit{\lambda })},\end{array}$$

where a_chl(λ) is the absorption of chlorophyll in m⁻¹ and λ and [Chl] as previously defined.

Previous estimates of a_chl(λ), a_w(λ) and b_w(λ) used by Morel and Maritorena (2001) cover values for wavelengths ranging between 300 and 700 nm. Therefore, we have combined and interpolated several sets of tables of coefficients in order to solve consistently α^W, ${\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}$ and ${\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{S}}$ across the same range of wavelengths (i.e., from 200 to 4000 nm). a_w(λ) has been derived from tables provided by Smith (1982) and Irvine and Pollack (1968), which span 200 to 800 and 800 to 4000 nm, respectively. a_chl(λ) has been derived from values published in Frigaard et al. (1996) which differ from those in Morel and Maritorena (2001; Fig. 2). b_w(λ) is estimated from seawater backscattering coefficients published in Morel and Maritorena (2001) that have been interpolated from 300–700 to 200–4000 nm with polynomial splines. Tabulated values for a_chl(λ), a_w(λ) and b_w(λ) are given in Table S1.

The difference in the contribution of the ocean interior reflectance to the ocean surface albedo for direct and diffuse radiation essentially stems from the incident direction of incoming radiation. In the case of ocean interior reflectance for direct incoming radiation, ${\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{W}}$, μ = cos (θ) whereas in the case of ocean interior reflectance for diffuse, ${\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{W}}$, μ = 0.676. This value is considered an effective angle of incoming radiation of 47.47^∘ according to Morel and Gentili (1991). Hence ${\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{W}}\left(\mathit{\lambda },\mathrm{Chl}\right)={\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{W}}\left(\mathit{\lambda },\arccos \left(\mathrm{0.676}\right),\mathrm{Chl}\right)$.

2.4 Computation of OSA

With the various components of OSA being now parameterized, the OSA for direct and diffuse radiation is estimated as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{OSA}}_{\mathrm{dir}}\left(\mathit{\lambda },\mathit{\theta },w,\mathrm{Chl}\right)=\left({\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}\left(\mathit{\lambda }\mathit{\theta }w\right)+{\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{W}}\left(\mathit{\lambda }\mathit{\theta }\mathrm{Chl}\right)\right)\\ \text{(14)}& {\displaystyle }& {\displaystyle }\left(\mathrm{1}-\mathrm{WC}\left(w\right)\right)+\mathrm{WC}\left(w\right){\mathit{\alpha }}^{\mathrm{WC}}\left(\mathit{\lambda }\right),{\displaystyle }& {\displaystyle }{\mathrm{OSA}}_{\mathrm{dif}}\left(\mathit{\lambda },\mathit{\theta },w,\mathrm{Chl}\right)=\left({\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{S}}\left(\mathit{\lambda }w\right)+{\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{W}}\left(\mathit{\lambda }\mathrm{Chl}\right)\right)\\ \text{(15)}& {\displaystyle }& {\displaystyle }\left(\mathrm{1}-\mathrm{WC}\left(w\right)\right)+\mathrm{WC}\left(w\right){\mathit{\alpha }}^{\mathrm{WC}}\left(\mathit{\lambda }\right).\end{array}$$

Since detailed atmospheric radiative transfer (e.g., Clough et al., 2005; Mlawer et al., 1997) are now part of the current generation of Earth system models, most radiative codes resolve radiation from near-ultraviolet (∼ 200 nm) to near-infrared (∼ 4000 nm) wavelengths. Here, we design our scheme to compute both the spectral and broadband OSA. To this effect, the scheme computes the OSA from λ₁=200 to λ₂=4000 nm with a resolution of 10 nm. The contribution of each wavelength interval dλ to OSA is weighted by its amount of solar energy under the standard solar spectra ASTM E-490 AM0 (Shanmugam and Ahn, 2007), E(λ) assumption as follows.

$$\begin{array}{}\text{(16)}& {\displaystyle }\mathrm{OSA}\left(\mathit{\theta },w,\mathrm{Chl}\right)={\displaystyle \frac{\underset{{\mathit{\lambda }}_{\mathrm{1}}=\mathrm{200}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}=\mathrm{4000}}{\int }}E\left(\mathit{\lambda }\right)\mathrm{OSA}\left(\mathit{\lambda },\mathit{\theta },w,\mathrm{Chl}\right)\mathrm{d}\mathit{\lambda }}{\underset{{\mathit{\lambda }}_{\mathrm{1}}=\mathrm{200}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}=\mathrm{4000}}{\int }}E\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}}\end{array}$$

Tabulated values for E(λ) are given in Table S1.

Finally for the total incoming radiation, the OSA can be written as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{OSA}\left(\mathit{\theta },w,\mathrm{Chl}\right)=(F_{\mathrm{dir}}{\mathrm{OSA}}_{\mathrm{dir}}\left(\mathit{\theta },w,\mathrm{Chl}\right)\\ \text{(17)}& {\displaystyle }& {\displaystyle }+F_{\mathrm{dif}}{\mathrm{OSA}}_{\mathrm{dif}}\left(\mathit{\theta },w,\mathrm{Chl}\right))/\left(F_{\mathrm{dir}}+F_{\mathrm{dif}}\right),\end{array}$$

where F_dir and F_dif are the downward surface fluxes of direct and diffuse radiation, respectively. OSA $\left(\mathit{\theta },w,\mathrm{Chl}\right)$ is then computed for each model ocean grid cell at each model time step. it should be noted that the SZA used in LMDZ is the average of the SZA during the daytime fraction of the time step.

3 Contribution of various OSA components

In this section, we analyze the geographical structure of OSA which is decomposed as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{A}_{\mathrm{dir}}\left(\mathit{\theta },w,\mathrm{Chl}\right)=\\ \text{(18)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\underset{{\mathit{\lambda }}_{\mathrm{1}}=\mathrm{200}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}=\mathrm{4000}}{\int }}E\left(\mathit{\lambda }\right)\left({\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}(\mathit{\lambda },\mathit{\theta },w)+{\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{W}}(\mathit{\lambda },\mathit{\theta },\mathrm{Chl})\right)\mathrm{d}\mathit{\lambda }}{\underset{{\mathit{\lambda }}_{\mathrm{1}}=\mathrm{200}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}=\mathrm{4000}}{\int }}E\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}},{\displaystyle }& {\displaystyle }{A}_{\mathrm{dif}}\left(\mathit{\theta },w,\mathrm{Chl}\right)=\\ \text{(19)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\underset{{\mathit{\lambda }}_{\mathrm{1}}=\mathrm{200}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}=\mathrm{4000}}{\int }}E\left(\mathit{\lambda }\right)\left({\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{S}}(\mathit{\lambda },\mathit{\theta },w)+{\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{W}}(\mathit{\lambda },\mathrm{Chl})\right)\mathrm{d}\mathit{\lambda }}{\underset{{\mathit{\lambda }}_{\mathrm{1}}=\mathrm{200}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}=\mathrm{4000}}{\int }}E\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}},\text{(20)}& {\displaystyle }& {\displaystyle }{A}_{\mathrm{WC}}={\displaystyle \frac{\underset{{\mathit{\lambda }}_{\mathrm{1}}=\mathrm{200}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}=\mathrm{4000}}{\int }}E\left(\mathit{\lambda }\right){\mathit{\alpha }}^{\mathrm{WC}}\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}{\underset{{\mathit{\lambda }}_{\mathrm{1}}=\mathrm{200}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}=\mathrm{4000}}{\int }}E\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }}}=\mathrm{0.174},\end{array}$$

where A_dir and A_dif are the broadband ocean surface albedos for direct and diffuse radiation in the absence of whitecap albedo; A_WC is the broadband albedo of whitecaps.

A_dir, A_dif and A_WC have been estimated from offline calculations using ERA-Interim forcing fields from 2000 to 2009 at monthly frequency (Dee et al., 2011) and chlorophyll climatology from SeaWiFS (Siegel et al., 2002). Compared to A_dir and A_dif, A_WC is constant in space; therefore its geographical structure arises from whitecap coverage (WC).

Figure 3a shows that A_dir displays a strong meridional gradient with high values over high-latitude oceans and low values over the tropical oceans. It confirms that the solar zenith angle is the prominent driver of A_dir. This albedo exhibits nonetheless geographical structures over the tropical oceans which are linked to the easterly wind regimes, suggesting that surface wind variability may have a small but noticeable influence on the ocean surface albedo for direct radiation.

Compared to A_dir, A_dif does not exhibit such a large meridional gradient (Fig. 3b). A_dif shows values close to 0.06. It displays nonetheless values > 0.06 over the subtropical gyres and values < 0.06 over the North Atlantic and the Southern Ocean in response to the 10 m wind speeds. Those patterns are related to surface wind patterns but also to the geographical structure of oligotrophic gyres with low chlorophyll values which reinforce the contribution of the ocean interior reflectance to surface albedo for the diffuse incoming radiation.

Figure 3 Maps of ocean surface albedo for (a) direct and (b) diffuse radiation in the absence of whitecaps, and map of whitecaps surface coverage (c). Estimates are derived from offline calculation using ERA-Interim forcing fields (Dee et al., 2011) from years 2000 to 2012 and SeaWiFS chlorophyll climatology (Siegel et al., 2002) over years 1998–2007. Whitecap albedo is constant in space and is equal to 0.174.

Figure 3c provides further insight into the regional influence of WC which display a broadband albedo of 0.174. Offline calculation of WC shows that whitecaps influence albedo for direct and diffuse radiation where westerly winds blow regularly, for example, in the Southern Ocean, the North Atlantic and the North Pacific. A weaker but noticeable influence is also found over the tropical oceans.

While A_WC is larger than A_dir and A_dif, the convolution of broadband albedo of the whitecaps and their coverage results in the maximum contribution of 0.003 to the broadband albedos for direct and diffuse radiation. However, its strong albedo makes the whitecaps an important factor at interannual and climate timescales. Indeed, this component of OSA for direct and diffuse radiation responds to the interannual variability of 10 m wind speed and also to climate change. Indeed, the contraction of Southern Ocean westerly winds (e.g., Böning et al., 2008) might induce subtle regional fluctuations in OSA that can feedback on the climate response.

4 Materials and methods

4.1 Observations

To assess model reliability to simulate realistic OSA, we compare fields to available observations. For those observations, we estimate the time-averaged OSA from the ratio between the time-averaged upwelling and downwelling shortwave radiative fluxes provided in those data sets.

At local scale, we use the CERES Ocean Validation Experiment (COVE; Rutledge et al., 2006) ground-based measurements. This ocean platform instrument, located at Chesapeake Bay, has provided continuous measurements of several radiative fluxes since 2001. In this study, we use measurements of upwelling and downwelling global (i.e., direct and diffuse) shortwave radiation averaged across several instruments (https://cove.larc.nasa.gov/instruments.html).

At global scale, we perform model evaluation with retrievals from the CERES satellite radiation measurements (Wielicki et al., 1996). CERES data provides estimates of global shortwave radiation at top of the atmosphere and at the surface. In the present study, we focus on surface estimates since our analyses aims at assessing the representation of ocean surface albedo.

4.2 Models

4.2.1 LMDZ v5A

LMDZ is an atmospheric general circulation model developed at the Laboratoire de Météorologie Dynamique. The version of this atmosphere model, known as LMDZ v5A, is described in detail in Hourdin et al. (2013); it is part of the main IPSL climate model used for CMIP5 and described in Dufresne et al. (2013; IPSL-CM5A). The atmospheric resolution is 96 × 95 horizontally and 39 layers vertically. The old OSA scheme in this version of LMDZ is based on the formulation of Larsen and Barkstrom (1977). It is parameterized in terms of μ as follows:

$$\begin{array}{}\text{(21)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}\left(\mathit{\theta }\right)={\displaystyle \frac{\mathrm{0.058}}{\mathit{\mu }+\mathrm{0.30}}}.\end{array}$$

Consequently, OSA varies between 0.0446 for a sun at zenith and 0.193 for a sun at the horizon. Direct and diffuse radiation are not distinguished, and only a broadband albedo is used in the visible spectrum (${\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{S}}={\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}})$.

In LMDZ, the partitioning between direct and diffuse light is derived from the presence of clouds in the atmosphere model grid cell.

We assess simulated OSA using an atmosphere-only simulation with prescribed radiative forcing (greenhouse gases, aerosols, land-cover change) and fixed sea surface temperature as recommended by CMIP5 (Taylor et al., 2012). LMDZ has been integrated from 1979 up to 2012 under this protocol.

Similarly to the observations, simulated OSA at a given frequency is derived from the ratio between the time-averaged upwelling and downwelling shortwave radiative fluxes at that frequency.

4.2.2 ARPEGE-Climat v6.1

ARPEGE-Climat v6.1 derives from ARPEGE-Integrated Forecasting System (IFS), the operational numerical weather forecast models of Météo-France and the European Centre for Medium-range Weather Forecasts (ECMWF). Compared to the version used in Voldoire et al. (2013), several improvements in atmospheric physics have been implemented. They consist of a new vertical diffusion scheme which solves a prognostic turbulent kinetic energy equation following Cuxart et al. (2000), an updated prognostic microphysics representing the specific masses of cloud liquid and ice water, rain and snow, as detailed in Lopez (2002), and a new convection scheme known as the Prognostic Condensates Microphysics Transport (PCMT; Guérémy, 2011; Piriou et al., 2007). ARPEGE-Climat v6.1 is implicitly coupled to the surface model called SURFEX (Masson et al., 2013), which considers a diversity of surface formulations for the evolution of four types of surface: land, town, inland water and ocean. The old OSA formulation implemented in SURFEX follows Taylor et al. (1996). This scheme enables the computation of ${\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}$ as a function of μ:

$$\begin{array}{}\text{(22)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{dir}}^{\mathrm{S}}\left(\mathit{\theta }\right)={\displaystyle \frac{\mathrm{0.037}}{\mathrm{1.1}{\mathit{\mu }}^{\mathrm{1.4}}+\mathrm{0.15}}}.\end{array}$$

Since this schema does not enable computation of ${\mathit{\alpha }}_{\mathrm{dif}}^{\mathrm{S}}$, it is set to a constant value of 0.066.

Like LMDZ, the partitioning of direct and diffuse radiation depends on the cloud cover in the atmospheric model grid cell.

Simulations performed with ARPEGE-Climat also consists of AMIP simulations as in LMDZ, except for sea surface temperature which relies on data recommended by CMIP6 (Eyring et al., 2016a). This simulation also spans 1979 to 2012.

Figure 4 Analytical solution for (a) direct and (b) diffuse ocean surface albedo as used in Taylor et al. (1996) and Larsen and Barkstrom (1977), and computed solution for the new interactive ocean surface albedo scheme, as a function of solar zenith angle. Analytical solution for direct and diffuse ocean surface albedo as derived from Payne (1972) formulation is also represented in both panels, because this parameterization is currently used in a number of state-of-the-art atmospheric and ocean models. Hatching depicts potential variations related to changes in 10 m wind speed and surface chlorophyll.

Analyses are complemented using another simulation of ARPEGE-Climat in which the resolved dynamics is nudged towards that of ERA-Interim. Nudging consists of restoring the model wind divergence and vorticity and the surface pressure towards those from ERA-Interim. The restoring timescale is 12 h for the wind divergence and surface pressure and 6 h for the wind vorticity. This simulation is employed hereafter as a kind of reference of what could be expected by the OSA parameterization if the wind spatiotemporal properties were “realistic”. In this case, only the direct-to-diffuse incident radiation partitioning remains tied to ARPEGE-Climat. This simulation replicates the chronology of the observed day-to-day variability of 10 m wind speed and hence is expected to be closer to the ground-based observations.

For those models simulations, simulated OSA at a given frequency is determined from the ratio between the time-averaged upwelling and downwelling shortwave radiative fluxes at that frequency.

5 Comparison of analytical calculation

In order to better understand changes in simulated OSA, we compare first the analytical solution of old and new interactive OSA schemes used in the two atmospheric models for both direct and diffuse radiation (Fig. 4). We also compare old and new interactive OSA schemes to Payne (1972) OSA scheme that is currently used in a number of atmospheric and ocean models such as NEMO (Madec, 2008).

Figure 4a shows that old OSA schemes for direct radiation differ in terms of response to solar zenith angle. Indeed, for a given solar zenith angle, the scheme used in LMDZ (Larsen and Barkstrom, 1977) leads to a greater OSA than that used in ARPEGE-Climat (Taylor et al., 1996). The shape of the response to variations in solar zenith angle suggests that the scheme used in ARPEGE-Climat leads to a slightly stronger meridional gradient in OSA than that used in LMDZ. Interestingly, the new scheme produces OSA values bracketed by those of old algorithms, except for small solar zenith angle. Under this condition, the effect of winds is to increase OSA up to 0.072. It also displays a greater response to variations in solar zenith angle which differs substantially from those given by the old schemes.

Compared to the Payne (1972) OSA scheme, old and new schemes used in the two atmosphere models exhibit a weaker meridional gradient in OSA (Fig. 4a). However, the meridional gradient as estimated by Payne (1972) is similar to that produced by Taylor et al. (1996) because their formulations solely differ by a coefficient; that is, 0.037 for Taylor et al. (1996) and 0.05 for Payne (1972).

Differences in OSA for diffuse radiation presented in Fig. 4b are noticeable. They clearly illustrate modeling assumptions in the old schemes. Indeed, old schemes have been built on ad hoc formulations. Neither Taylor et al. (1996) nor Larsen and Barkstrom (1977) have provided a differentiated OSA for direct and diffuse radiation. This is why OSA for diffuse radiation is set to 0.06 (corresponding to the angular average of the OSA for direct radiation) in ARPEGE-Climat, whereas that of LMDZ is equal to the OSA for direct radiation from Larsen and Barkstrom (1977).

Figure 5 Ocean surface albedo at COVE station (36.905^∘ N, 75.713^∘ W) from 2005 to 2009. Panel (a) compares daily mean time series of ocean surface albedo as derived from ground-based observations (in red) and as reconstructed with ARPEGE-Climat nudged toward ERA-Interim (dark blue). Panel (b) displays, for the sake of clarity, time series of daily mean ocean surface albedo smoothed using a 5-day moving average for both observations and model results. All daily mean time series from 2001 to 2015 are displayed in Fig. S1. Ocean surface albedo simulated by ARPEGE-Climat (in blue) and LMDZ (in green) using the old and new interactive scheme are indicated with dashed and solid lines, respectively.

Figure 4b shows that the new interactive scheme displays features similar to the diffuse OSA used in ARPEGE-Climat or that estimated from Payne (1972). This scheme produces nonetheless slightly larger values which can fluctuate in response to other drivers. The old OSA for diffuse radiation employed in LMDZ responds to variations in solar zenith angle although it should not. Errors related to this erroneous representation of OSA for diffuse radiation are also modulated by the partitioning between direct and diffuse radiation estimated by the atmospheric model.

Figure 6 Probability density function of daily mean ocean surface albedo at COVE station (36.905^∘ N, 75.713^∘ W) derived from daily mean time series over years 2001 to 2013. Ocean surface albedo derived from ground-based observations and as reconstructed with ARPEGE-Climat nudged toward ERA-Interim are indicated in red and dark blue, respectively. Ocean surface albedo simulated by ARPEGE-Climat (in blue) and LMDZ (in green) using the old and new interactive scheme are indicated with dashed and solid lines, respectively.

6 Evaluation at COVE station (36.905^∘ N, 75.713^∘ W)

In this section, we employ COVE daily data to assess the simulated OSA by both atmospheric models at local scale. OSA is computed here as the ratio of up-to-down radiation fluxes averaged at daily resolution for both ground-based observations and models. Such an evaluation is fundamental because it relies on direct ground-truth observations over the ocean surface and hence provides a more accurate assessment of the OSA scheme as compared to the global-scale satellite-derived estimates developed in the following sections.

Figure 5 shows how well the model using old and new interactive OSA scheme behaves at daily frequency compared to the ground-based observations at COVE station from 2001 to 2009. Figures 5 and S1a clearly show that both old OSA schemes of ARPEGE-Climat and LMDZ fail at replicating day-to-day OSA variations at the COVE station. Comparatively, Figs. 5 and S1b emphasize how much the new interactive scheme improves OSA as simulated by both atmospheric models. Indeed, the simulated OSA values are now consistent with observation at COVE station, with temporal correlation greater than 0.3. However, the models fail at replicating the large OSA values occurring during the winter in ground-based observations.

Figure 7 Mean seasonal cycle of ocean surface albedo at COVE station (36.905^∘ N, 75.713^∘ W) derived from daily mean time series over years 2001 to 2013. Ocean surface albedo derived from ground-based observations and as reconstructed with ARPEGE-Climat nudged toward ERA-Interim are indicated in red and dark blue, respectively. Ocean surface albedo simulated by ARPEGE-Climat (in blue) and LMDZ (in green) using the old and new interactive scheme are indicated with dashed and solid lines, respectively. For comparison, the mean seasonal cycle of ocean surface albedo at COVE as simulated by available CMIP5 models is represented by thin grey lines.

Those findings are reinforced when we compare the probability density function (pdf) estimated from daily mean OSA as simulated by models against that derived from ground-based observations (Fig. 6). This analysis provides further insight into how old and new interactive OSA schemes behave at COVE station. Figure 6 confirms that old schemes fail at capturing the day-to-day variations in OSA. Indeed, day-to-day variations in OSA estimated from old schemes arise from day-to-day variations in SZA and to a lesser extent to variations in direct-to-diffuse ratio of incident radiation which is related to the cloud cover. As shown in Fig. 4, old OSA schemes crudely represent diffuse albedo. Therefore, errors in direct-to-diffuse ratio of incident radiation incurs errors in the simulated OSA. Consequently, day-to-day variations are better reproduced when the albedo for diffuse radiation is realistically simulated (Fig. 6). In particular, the new interactive scheme captures the minimum OSA values occurring during the summer which are lower than 0.06.

At seasonal scale, OSA estimated from averaged radiative fluxes agrees with the above-mentioned findings for ground-based observations and models. Figure 7 clearly shows that old OSA schemes do not capture seasonal variations of observed OSA. Correlation between observation-derived OSA and that simulated by both models is 0.32 for ARPEGE-Climat and 0.28 for LMDZ, which is very low, indicating an unrealistic representation of OSA. Comparison with CMIP5 atmosphere models shows that OSA as simulated by ARPEGE-Climat or LMDZ is in the range of CMIP5 models (0.04–0.17), confirming the large uncertainties related to simulated OSA in a state-of-the-art climate model. While several CMIP5 models replicate seasonal variation in OSA, most of them exhibit large biases in simulated OSA compared to the observation-based estimate. Only ACCESS1-3, BNU-ESM, HadCM3, MIROC-ESM and MIROC-ESM-CHEM display a mean seasonal cycle of OSA comparable to the observation-based estimate at COVE station. For ARPEGE-Climat, this erroneous representation of OSA at seasonal-scale leads, at least for this location, to a systematic bias in the surface energy budget of +3 W m⁻² in winter and −1.5 W m⁻² in summer. It is thus likely that large deviation in OSA as simulated by CMIP5 leads to substantial errors in energy flow at the air–sea interface.

Figure 7 shows significant improvements in the simulated OSA in both models using the new interactive scheme. In both models, the simulated seasonal cycle of OSA replicates the minimum observed during the summer, although with the new interactive OSA scheme, both models do not capture large values of OSA of 0.10 occurring during the winter. That said, model–data comparison shows that correlation with observations has been improved. Indeed, correlation between observed and simulated daily values over a mean yearly cycle has increased from 0.23 to 0.84 in LMDZ to 0.32 to 0.86 in ARPEGE-Climat.

Figure 8 Decadal-mean climatology ocean surface albedo as (a) estimated from CERES satellite observations (Wielicki et al., 1996) and as simulated by LMDZ (b, d) and ARPEGE-Climat (c, e). In panels (b) and (c), LMDZ and ARPEGE-Climat use old ocean albedo schemes, that is, Taylor et al. (1996) and Larsen and Barkstrom (1977), respectively. In panels (d) and (e), LMDZ and ARPEGE-Climat use employ the new interactive ocean surface albedo scheme. Decadal-mean climatology is derived from radiative fluxes averaged over years 2001 to 2014 for CERES estimates and 2000 to 2012 for both climates models.

Although improved, Figs. 5, 6 and 7 show that new interactive OSA scheme seems to suffer from a systematic bias in winter and miss OSA values greater than 0.10. This is supported by the fact that this systematic bias is displayed for all model estimates independently from the atmospheric physics and dynamics (i.e., LMDZ, ARPEGE-Climat and nudged OSA). That being said, some other possible reasons can explain such deviations between models and ground-based data. First, current atmosphere models suffer from systematic errors in the ratio of direct-to-diffuse radiation which could be related to bias in cloud cover or aerosol optical thickness (as shown in Fig. S2). A greater-than-observed atmospheric optical depth in winter may favor the diffuse path over the direct path, resulting in a lower-than-observed OSA. Second, coarse-resolution atmospheric models are not able to replicate the mesoscale meteorological and oceanic conditions at this very location. Differences in surface wind between the models and field conditions can increase the contribution of whitecap albedo with the respect to that of the Fresnel reflectance. Third, local ocean conditions and the presence of ocean waves resulting from remote wind (i.e., swells) are not simulated by the atmospheric models. This would lead to an underestimation of the contributions of both whitecaps and Fresnel reflectance.

7 Global-scale evaluation

7.1 Climatological mean

This section details the evaluation of OSA at the global scale using global satellite product that are routinely used in Earth system model evaluation (Eyring et al., 2016b; Gleckler et al., 2008). We thus use OSA retrieved from CERES surface product to assess the simulated OSA by ARPEGE-Climat and LMDZ.

Figure 8 presents geographical pattern of OSA as simulated by ARPEGE-Climat and LMDZ using Taylor et al. (1996) and Larsen and Barkstrom (1977) schemes, respectively. These old OSA schemes were used during CMIP5 and thus give an idea of the errors in the models' radiative budgets.

Globally, the simulated OSA overestimates the CERES-derived estimate (∼ 0.058) by about 0.007. However, the most striking feature is the substantial differences in the meridional structure. CERES-derived OSA shows maximum values over the high-latitude oceans and minimum values over the tropical oceans. None of the models using old schemes are able to capture this meridional structure. Deviations are particularly high for LMDZ, which only barely replicate maximum OSA over high-latitude oceans and minimum OSA over tropical oceans. Both models exhibit poor spatial correlation with −0.03 for LMDZ and 0.40 for ARPEGE-Climat. Model–data errors in OSA mirror model bias in surface upwelling shortwave radiation, which amounts to ∼ 7 W m⁻² over the tropical oceans compared to CERES.

The new interactive scheme improves favorably the comparison with observations (Fig. 8). Indeed global mean OSA is equal to 0.062 for LMDZ and 0.057 for ARPEGE-Climat, which better matches the value derived from CERES data. As such, the model bias in surface upwelling shortwave radiation has been reduced by ∼ 1 W m⁻² on average over the ocean and by up to ∼ 5 W m⁻² over the tropical oceans.

Both models capture the meridional structure of the OSA with spatial correlations of about 0.82 for LMDZ and 0.86 for ARPEGE-Climat. Nonetheless, the simulated OSA displays some biases. In LMDZ and ARPEGE-Climat, the modeled OSA over the North Atlantic is slightly overestimated and shifted to the south. Major differences between simulated OSA are noticeable over the tropical oceans, where models differ in terms of zonal structure. LMDZ displays OSA of ∼ 0.06 over eastern boundary upwelling systems, which is slightly too high compared to CERES. Differences in the OSA geographical structure between ARPEGE-Climat and LMDZ arise from differences in 10 m wind speed (Fig. S3) and direct-to-diffuse incident radiation as determined from the simulated cloud cover (Fig. S4). Large-scale deviations between models and observations seem to be related to differences in 10 m wind fields (Fig. S3). Model–data deviations in OSA at the regional scale rather mirror biases in total cloud cover (Fig. S4). This is especially clear over low-latitude oceans where LMDZ overestimates OSA over the eastern boundary upwelling systems where LMDZ overestimates the cloud cover (Fig. S4). This result is expected since over the low-latitude oceans the contribution of diffuse OSA is stronger than that of direct OSA (Figs. 3 and 4).

7.2 Seasonal variability

Figure 9 compares the simulated and CERES-derived OSA on the seasonal scale. This timescale matters for accurately modeling the Earth's climate because the flow of incoming radiation fluctuates up to 1 order of magnitude between winter and summer at high latitudes.

Figure 9a, b and c show that both models using old OSA schemes only just reproduce the seasonal cycle of OSA derived from CERES. This is particularly the case for LMDZ, which produced an unrealistic seasonal cycle for OSA. LMDZ fails at simulating maximum OSA during the winter of both hemispheres. Instead, extreme values of simulated OSA occur at 50^∘ N and 50^∘ S during the summer. Simulated OSA in ARPEGE-Climat does not present these features but is biased high during all seasons.

Figure 9 Hovmöller diagram representing the zonally averaged ocean surface albedo as a function of month. The various panels display the ocean surface albedo as (a) estimated from CERES satellite observations (Wielicki et al., 1996) and as simulated by (b) LMDZ and (c) ARPEGE-Climat using old ocean albedo schemes, that is, Taylor et al. (1996) and Larsen and Barkstrom (1977), respectively. Panels (d) and (e) show OSA as simulated by the new interactive ocean surface albedo scheme for LMDZ and ARPEGE-Climat, respectively. Monthly means are derived from radiative fluxes averaged over years 2001 to 2014 for CERES estimates and from years 2000 to 2012 for both climates models.

With the new interactive scheme, the seasonal OSA is improved in both models (Fig. 9d, e). The simulated OSA matches that derived from CERES at seasonal scale, with high values during the winter and low values between 30^∘ S and 30^∘ N. Improvement is especially noticeable for LMDZ which captures the observed seasonal cycle of OSA.

However, a few errors remain in the simulated OSA. In LMDZ, OSA is slightly too high compared to CERES (∼ 0.002) in boreal and austral summer. Nonetheless, simulated OSA reproduces realistic OSA values in the tropics (∼ 0.052). In ARPEGE-Climat, instead, the simulated OSA seems slightly too low compared to CERES (∼ −0.002). This leads ARPEGE-Climat to overestimate the fraction of low-OSA ocean. Interestingly this bias solely concerns the tropical oceans. Indeed, simulated OSA over high-latitude oceans displays realistic features at the seasonal scale. The fact that errors in ARPEGE-Climat and LMDZ are of different signs tends to suggest that the new interactive scheme is not intrinsically biased. It rather points to biases in driving fields such as the surface wind speeds or the ratio between direct and diffuse shortwave radiation simulated by either ARPEGE-Climat or LMDZ (Fig. S2).

8 Conclusions

In this paper, we have detailed a new interactive scheme for ocean surface albedo suited for Earth system models. This scheme computes the ocean surface albedo accounting for the spectral dependence (across a range of wavelengths between 200 and 4000 nm), the characteristics of incident solar radiation (direct of diffuse), the effects of surface winds, chlorophyll content and whitecaps and the canonical solar zenith angle dependence. This scheme enables improved air–sea exchange of solar radiation. It thus provides a much more physical basis to resolve the radiative transfer at the interface between the atmosphere and the upper ocean and offer a suite of processes that are included in complex stand-alone ocean radiative transfer software such as HYDROLIGHT (Ohlman et al., 2000). This work can be extended to include coupling to an ocean wave model that would provide a more realistic distribution of ocean surface state.

Although direct and diffuse albedos were included in the old ocean albedo schemes of the two atmospheric models used here, our results demonstrate that their assumptions employed for diffuse albedo (i.e., fixed values or equal to the direct albedo) are not realistic. The new interactive scheme improves its representation which leads to substantially reduced model–data error in ocean surface albedo over the low-latitude oceans.

Comparison to available data set shows, for at least two state-of-the-art climate models, a noticeable improvement in terms of simulated ocean surface albedo compared to their old ocean surface albedo schemes. At the global scale, the geographical pattern of simulated ocean surface albedo has been improved in both models. The simulated seasonal cycle also shows a noticeable improvement, especially in LMDZ, with a better correlation to CERES data (up to 0.8). At the local scale, simulated ocean surface albedo also fits ocean surface albedo derived from ground-based radiative measurements at daily resolution with an improved correlation up to 0.8.

Compared to old schemes, the new interactive scheme is more complex and induces a small increase in elapsed model time of about 0.2 %. Although noticeable, this increase does not preclude centennial-long simulation or high-resolution model simulations.

Improved ocean surface albedo might lead to differences in the simulated climate or marine biogeochemistry dynamics which will be assessed in future work. Indeed, a difference of about 1 % of simulated ocean surface albedo for a global mean irradiance of ∼ 180 W m⁻² can induce a deviation in the energy flow of the Earth system comparable to the impact of land-cover changes over land (Myhre et al., 2013).