© ESO 2019

1 Introduction

Liquid water is considered to be a prerequisite for life as we know it. The presence of liquid water on the surface of an exoplanet will depend on the planet’s surface temperature and therefore on the separation between the planet and its host star, the stellar luminosity, the surface albedo, the composition and structure of the planetary atmosphere (which also strongly influence the planet’s albedo) and the atmospheric surface pressure. While the amount of stellar flux that is incident on the planet can be measured, the surface albedo, and the atmospheric properties are more difficult to determine. Yet they are crucial factors in determining a planet’s habitability.

Some information about the composition and structure of the upper atmosphere of exoplanets can be derived using transit spectroscopy, where the spectral variation of the stellar flux that is transmitted through the atmosphere of a planet during its transit, is measured. Indeed, detections have been reported of water vapor in the atmospheres of hot Jupiter-like planets (e.g., Tinetti et al. 2007; Deming et al. 2013; Kreidberg et al. 2014) and a low-mass Neptune sized planet (Fraine et al. 2014). A large number of transiting exoplanets is already known, mostly thanks to the space observatories Convection, Rotation and planetary Transits (CoRoT) by the European Space Agency (ESA) and the National Centre for Space Studies (CNES), and, of course, especially Kepler by the National Aeronautics and Space Administration (NASA). New and upcoming space observatories, such as the Transiting Exoplanet Survey Satellite (TESS), NASA, launched in 2018, Ricker et al. (2015) and PLAnetary Transits and Oscillations of stars (PLATO), ESA, planned for launch in 2026, Rauer et al. (2014), and ground-based telescope networks such as the Wide Angle Search for Planets (WASP; Ricker et al. 2015) and the Multi-site All-Sky CAmeRA (MASCARA; Talens et al. 2017) are designed to search for transiting exoplanets around bright, nearby stars. Such transiting exoplanets will be excellent targets for follow-up transit spectroscopy measurements with, for example, the James Webb Space Telescope (JWST), NASA, and the Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL), ESA, planned for launch after PLATO.

Transit spectroscopy, however, does not provide information about the lower atmospheric layers because the optical path-length through these layers is too large to allow the transmission of starlight, in particular in the presence of clouds. And, rocky exoplanets are so small compared to their star, and their atmospheresare so thin compared to the planetary radius, that obtaining a reliable signal-to-noise ratio (S/N) could take several tens to hundreds of planetary transits. For rocky exoplanets in the habitable zones of solar-type stars, that would thus take several tens to hundreds of Earth-years (Kaltenegger & Traub 2009).

A secondary eclipse spectrum, derived from the change in the total flux spectrum when a planet starts to disappear behind and/or reappears from behind its star as seen from the observer, could provide information about the planetary albedo, and thus potentially about the planet’s surface albedo, the atmospheric gaseous scattering, and the presence of clouds (not their actual composition), although all blended together due to the lack of spatial resolution. Unless the atmospheric gaseous optical thickness, the cloud coverage fraction and the cloud physical properties (especially cloud optical thickness) are known, a secondary eclipse spectrum will thus not provide information about the planet’s surface albedo (thus, without the atmosphere) and the probability of liquid surface water (which is expected to have a very small albedo) will thus remain elusive in secondary eclipse measurements.

For identifying liquid surface water, measuring the light of the parent star that is reflected by an orbiting exoplanet, appears to be the most straightforward tool. The spectral and temporal variations of this reflected light as the planet orbits its parent star, holds information about the composition and structure of the planetary atmosphere and/or the surface (if present). This is particularly true for the state of polarization of the reflected light (see Wiktorowicz & Stam 2015, and references therein). The angular distribution of the degree of polarization of singly scattered light is much more sensitive to the size, composition and shape of the scattering particles in the exoplanet’s atmosphere than the total flux (unpolarized + polarized fluxes) is (Hansen & Travis 1974), and in addition, the angular distribution of the polarized signal is preserved when multiple light is added. In particular, the phase angles¹ where the degree of polarization of the reflected starlight attains a local maximum value or where it equals zero (the so-called neutral points) can be directly translated into the scattering angles where the degree of polarization of the singly scattered light attains a local maximum value or zero². And, because the degree of polarization is a relative measure (i.e., the ratio of the polarized flux to the total flux), it is insensitive to parameters such as distance, planetary radii, and the incoming flux. The degree of polarization is also relatively insensitive to the Earth’s atmosphere for ground-based observations, as the transmission of light through the Earth’s atmosphere is virtually independent of the state of polarization of the light (F, Q, and U are equally affected by scattering and/or absorption) and it is insensitive to several instrumental parameters (depending on the instrument design), such as the degradation of filters or lenses in time.

A famous application of these advantages of polarimetry is the characterisation of the particles composing the main cloud deck on Venus by Hansen & Hovenier (1974). They used disk-integrated polarization measurements at wavelengths across the 0.34–1.00 μm region of Venus with Earth-based telescopes at a wide range of phase angles to determine that the clouds particles consist of 75% sulfuric acid solution, and have an effective radius of 1.05 μm with an effective variance of 0.07. Because, as seen from Earth, Venus is an inner planet, it can be observed at phase angles ranging from almost 0° to almost 180°, while outer Solar System planets can only be observed at small phase angles (except of course from a spacecraft in the vicinity of an outer planet). The phase angle range at which an exoplanet can be observed depends on the orbital inclination angle i (defined such that i = 0° for a face-on orbit), as seen fromEarth: 90° − i ≤ α ≤ 90° + i. In particular, transiting exoplanets, which TESS (Ricker et al. 2015) and PLATO (Rauer et al. 2014) are and will be searching for, will be observable across a broad phase angle range with direct detection techniques.

Simulations of orbital variations of the (polarized) flux and degree of polarization of starlight that is reflected by habitable exoplanets explore the potential observables of these planets, are crucial for the design and optimization of future instruments and telescopes aimed at the direct detection of planetary signals, and can be used to devise observational strategies (such as integration times and temporal coverage) and (future) data analysis algorithms. An exampleof a direct detection instrument that is designed to measure polarization of starlight reflected by exoplanets is the Exoplanet Imaging Camera and Spectrograph (EPICS) for the European Extremely Large Telescope (E-ELT; Kasper et al. 2010). EPOL is the imaging polarimeter for EPICS (Keller et al. 2010). NASA’s Large UV Optical Infrared Surveyor (LUVOIR) space telescope concept (LUVOIR Team 2018) includes POLLUX, a polarimeter for the ultraviolet wavelength region (Bouret et al. 2018). Another concept study by NASA is the Habitable Exoplanet Observatory (HabEx; Gaudi et al. 2018), that, as far as we know, does not (yet) include a polarimeter.

In this paper, we focus on the numerical simulations of the total and polarized fluxes and the degree of polarization of starlight that is reflected by Earth-like exoplanets that are covered by a liquid ocean with surfaces that are rough due to the presence of waves. Fresnel reflection of light by an ocean on an exoplanet results in two main observable phenomena: a glint on the ocean surface, whose size increases with increasing surface roughness (i.e., increasing wind speed), and a maximum degree of polarization at the Brewster angle. Indeed, as the strength of both phenomena depends on the local reflection angles, they are strongly phase angle dependent.

The signatures of an ocean in the orbital variations of the total flux F of starlight that is reflected by an exoplanet have been studied before. Robinson et al. (2010) modeled the phase curves of the total flux of the Earth as an exoplanet, and Oakley & Cash (2009), Cowan et al. (2009) and Lustig-Yaeger et al. (2018) aim to reconstruct a surface type map from modeled total flux phase curves at a single or multiple wavelength(s). Those papers, however, did not include the polarization of the reflected starlight. Stam (2008) computed the total flux F and degree of polarization P as functions of the wavelength and planetary phase angle for Earth-like exoplanets some of which with a Fresnel reflecting ocean, but did not include wind-driven waves on the ocean surface, thus yielding an infinitely narrow beam of reflected light representing the stellar glint (the spatial distribution of this glint was lost upon integration of the reflected starlight across the illuminated and visible part of the planetary disk). Williams & Gaidos (2008) presented the phase curves of an ocean planet with a rough surface and found strong signatures of Fresnel reflection in the phase angle dependence of the degree of polarization P. They, however, neglected atmospheric contributions such as Rayleigh scattering and scattering by clouds. Zugger et al. (2010, 2011a,b) did include atmospheric Rayleigh scattering, and the reflection by clouds modeled as horizontal, Lambertian (i.e., isotropic and nonpolarizing) reflecting surfaces³ (not only ignoring gas below and in the clouds, but also above the clouds), and scattering by maritime aerosols (modeled as spherical particles) in the cloud-free atmosphere above a rough ocean surface. They found that the ocean signature in P is limited by the scattering by maritime aerosols and clouds.

The model planets of Zugger et al. (2010, 2011a,b) are horizontally homogeneous, that is, the atmosphere and surface are invariant with latitude and longitude. Therefore, they could, for example, not study the signatures of the glint appearing and hiding behind patchy clouds, which are clouds of varying horizontal shapes distributed across the planetary disk. Rossi & Stam (2017) showed the features that could potentially be observed in the flux and polarization phase curves of horizontally inhomogeneous ocean planets using various cloud cover types. Their ocean surface was, however, still a flat air-water interface, without waves, and thus with an infinitely narrow glint. Also, Karalidi et al. (2012a) and Karalidi & Stam (2012) presented computed flux and polarization phase curves for horizontally inhomogeneous, cloudy Earth-like exoplanets, focusing on the signatures of the clouds, and representing the ocean by a flat, zero albedo (i.e., black) surface (without the glint).

In this paper, we present the computed planetary phase curves of the total and polarized fluxes and the degree of polarization of starlight that is reflected by oceanic exoplanets with rough, wavy surfaces beneath terrestrial-type atmospheres. We fully include scattering of (unpolarized and polarized) light within the water body (i.e., below the rough ocean surface) and the reflection by wind-generated foam on the ocean surface. We compute the phase curves for various wind speeds and surface pressures, at wavelengths ranging from 350 to 865 nm. We use two types of model ocean planets: cloud-free planets, and cloudy planets. The clouds are patchy and embedded within the gaseous atmosphere (there is thus gas below, within and above the clouds), just like those on Earth, and our model clouds consist of liquid water droplets.

This paper is structured as follows. In Sects. 2.1 and 2.2, we provide our definitions and describe our numerical method to compute both the total and polarized fluxes, and the degree of polarization P of starlight that is reflected by our model ocean planets. In Sect. 2.3, we describe the physical and optical properties of the atmosphere and the ocean surface and water body below. In Sect. 3, we present the computed reflected total and polarized fluxes and P for our model planets for various atmospheric and surface parameters, and we describe the influence of these parameters on the planet’s color change in total flux F, polarized flux Q, and P. In Sect. 4, we discuss the color change and with that the detectability of an ocean beneath clouds. In Sect. 5, finally, we summarize our results.

2 Numerical method

2.1 Definitions of fluxes and polarization

We describe light as the following column vector (see, e.g., Hansen & Travis 1974; Hovenier et al. 2004) $$F=\left[\begin{array}{c}F\\ Q\\ U\\ V\end{array}\right],$$(1)

where F describes the total flux, Q and U the linearly polarized fluxes, and V the circularly polarized flux. The dimension of these fluxes is Wm⁻² or in Wm⁻³, when accounting for the dependence on the wavelength λ, but because of our normalization (see below), we do not use these dimensions explicitly. Linearly polarized fluxes Q and U are defined with respect to a reference plane. For our computations of starlight that is reflected by a planet as a whole, we use the planetary scattering plane, that is, the plane that contains the centers of the star, the planet, and the observer, as the reference plane. Fluxes Q and U can be redefined with respect to another reference plane, such as an instrument’s optical plane, using a so-called rotation matrix (see Hovenier & van der Mee 1983).

We assume that the starlight that is incident on the planet is unidirectional and unpolarized. The latter is based on the very small disk-integrated polarized fluxes of solar-type stars, as measured for active and inactive FGK–stars by Cotton et al. (2017) (with, respectively, 23 ± 2.2 ppm and 2.9 ± 1.9 ppm polarized flux) and Kemp et al. (1987) for the Sun itself, and the very small degree of stellar linear polarization (on the order of 10⁻⁶) that is expected from symmetry-breaking disturbances such as spots and/or flares (Berdyugina et al. 2011; Kostogryz et al. 2015). The incident starlight is thus described as F₀ = F₀1, with πF₀ the stellar flux measured perpendicular to the direction of propagation, and 1 the unit column vector.

Starlight that is reflected by a planet will become polarized by scattering by gases and aerosol or cloud particles in the planetary atmosphere and/or by reflection off the surface. We define the degree of polarization of the reflected starlight as $$P=\frac{\sqrt{Q^2+U^2}}{F}.$$(2)

Unlike Q and U, P is independent of a reference plane. We ignore the circularly polarized flux V of the reflected starlight. This flux is expected to be very small (for examples for rocky exoplanets, see Rossi & Stam 2018; García Muñoz 2015, and references therein), and ignoring circular polarization in the radiative transfer computations yields only very small errors in the total and linearly polarized fluxes (Stam & Hovenier 2005).

2.2 Calculating reflected starlight

We compute starlight that is reflected by different types of exoplanets, both the total flux and polarization, and for a range of wavelengths λ. Given a planet orbiting a star, with πF₀1 the flux vector that is incident on the planet, we compute the flux vector πF that is reflected by the planet as a whole and that arrives at an observer situated at a large distance d of the planet and star, by integrating the locally reflected light across the illuminated and visible part of the planetary disk. We perform this integration using a grid of pixels that cover the planetary disk (see Rossi et al. 2018), and by summing the fluxes reflected by each of the N pixels that are both illuminated by the parent star and visible to the observer, as follows: $$\pi F(\lambda ,\alpha )=\frac{F_0(\lambda )}{d^2}{\sum }_{i=1}^N{\mu }_i{\mu }_{0i}L({\beta }_i)R_i(\lambda ,{\mu }_i,{\mu }_{0i},{\varphi }_i-{\varphi }_{0i})1d{O}_i.$$(3)

Here α is the planetary phase angle, which is defined as the angle between the star and the observer measured from the center of the planet (at α = 0°, the planet is thus behind the star as seen from the observer, and at α = 180°, the planet is precisely in front of the star). Angle α and the location on the planetary disk determine the local illumination and viewing angles across the planet: μ = cos θ, with θ the local viewing zenith angle, μ₀ = cos θ₀, with θ₀ the local illumination angle, and ϕ − ϕ₀, the local azimuthal angle (see de Haan et al. 1987, for the precise definitions of these angles).

Matrix R_i in Eq. (3), is a local reflection matrix, that is, it describes how the area that pixel i covers on the planet, reflects the incident light. This matrix depends on the composition and structure of the local atmosphere and the underlying surface (see Sect. 2.3), on the wavelength, as the optical properties of the atmosphere and surface are wavelength dependent, and on the local illumination and viewing angles (hence, on phase angle α and/or the location of the pixel on the planet). Matrix L is a rotationmatrix (see Hovenier & van der Mee 1983) that rotates each locally reflected flux vector, that is defined with respect to the local meridian plane through the local zenith and the local direction toward the observer, to the planetary scattering plane. Angle β_i is the local rotation angle. Furthermore, dO_i is the area that grid pixel i covers on the 3D-planet. The area dO_i a pixel covers on the 3D-planet depends on the local viewing angle θ_i. Indeed, μ_idO_i equals the pixel size, which we choose the same for all pixels.

With Eq. (3), we can compute the flux vectors reflected by locally horizontally homogeneous planets, that is, where the surface and atmosphere are the same for each pixel across the planetary disk, and by horizontally inhomogeneous planets, that is, where the surface and atmosphere can be different for different pixels. The error that results from using the summation of Eq. (3) for the integration across the planetary disk depends on the size of the pixels, on the horizontal variations of the surface and atmosphere, and on the phase angle α (see Rossi et al. 2018). For our computations, we use a minimum of 40 pixels across the planetary equator for every α, which provides a good relation between the integration error and the computing time.

To compute πF, we first assign an atmosphere and surface model to each pixel on the planetary disk. Then, given the phase angle α, and thus the local illumination and viewing angles, the local matrices R_i are computed and rotated using the local rotation angles β_i. Then, the summation in Eq. (3) is evaluated. The computation of the matrices R_i is performed with an efficient adding–doubling algorithm (de Haan et al. 1987; Stam 2008), extended with reflection by and transmission through a rough ocean surface and scattering in the ocean body (see Sect. 2.3), fully including polarization for all orders of scattering. Our adding–doubling algorithm is a multistream radiative transfer code. To accurately capture the angular variations in the reflected light signals, we use 130 streams for our simulations and a Fourier-series representation for the azimuthal direction (for details, see de Haan et al. 1987; Rossi et al. 2018).

We normalize computed reflected flux vectors such that at α = 0°, the total reflected flux equals the planet’s geometric albedo A_G (see, e.g., Stam et al. 2006; Rossi et al. 2018). The total and polarized fluxes presented in this article can thus straightforwardly be scaled to the relevant parameters of a given star-planet system. The degree of polarization P is a relative measure and thus independent of the stellar flux, planetary radius, and/or the distance to the planet.

2.3 The planetary surface–atmosphere models

Each pixel on the planetary disk has an atmosphere bounded below by a surface. The local atmosphere and surface are horizontally homogeneous, and rotationally symmetric around the local vertical. In principle, every pixel can have a unique atmosphere and surface model. In practice, we limit the number of different models across the disk to two. Detailed descriptions of the atmosphere and surface models are given below.

2.3.1 The atmosphere

Our model atmosphere is Earth-like. Locally, the planetary atmosphere is composed of a stack of horizontally homogeneous layers, each containing either only gas or gas with a cloud added. The gas molecules scatter as anisotropic Rayleigh scatterers with a wavelength independent depolarization factor equal to 0.03 (see Stam 2008, for details regarding their scattering properties). We perform our computations at wavelengths where absorption by the gas is small, and can safely be ignored. The standard surface pressure is 1 bar, the acceleration of gravity 9.81 m s⁻¹², the mean molecular mass 29 g mole⁻¹, and we assume hydrostatic equilibrium, which yields a total atmospheric gas optical thickness of 0.096 at λ = 550 nm (see Stam 2008, and references therein for the relevant equations).

Clouds strongly affect the signal of the glint on a planet, because they prevent light from reaching the ocean and they will block the light reflected by the ocean from reaching the observer. Our model clouds are composed of spherical particles, which we choose to bemade of liquid water. The particles are distributed in size according to a two-parameter gamma distribution (see Hansen & Travis 1974; De Rooij & van der Stap 1984), with the cloud particle effective radius equal to 10 μm and the effective variance equal to 0.1. The optical properties of the cloud particles are computed using Mie-theory as described by De Rooij & van der Stap (1984) and a wavelength independent refractive index of 1.33. Table 1 lists the parameters describing the atmosphere model including the clouds. We would like to point out that the reflection by cloudy regions on the planet models used by Zugger et al. (2010, 2011a,b) is described by the reflection bywhite, Lambertian surfaces, and scattering by gas below, in and above the clouds is ignored by Zugger et al. (2010, 2011a,b).

The influence of clouds on the planet signals will not only depend on the cloud physical properties, but also on the location of the cloudy pixels on the planetary disk. We will investigate the effect of the cloud cover on the planet’s fluxand polarization signals for patchy clouds, modeling them as horizontal slabs with optical thickness b^a that are distributed across the planetary disk. The location of the clouds is determined by the cloud fraction f_c (defined on the full disk, not only on the illuminated part) and asymmetry parameters that describe a zonal-oriented pattern similar to that observed on Earth (for the algorithm, see Rossi & Stam 2017). Examples for different cloud coverage fractions f_c can be seen in Fig. 7.

2.3.2 The surface

The reflection by the surface is either Lambertian, that is, isotropic and completely depolarizing, or oceanic, which means that the reflection is described by a Fresnel reflecting and transmitting air–water interface and subinterface water layers.

The air–water interface can be rough due to wind-driven waves. The roughness is modeled with a Gaussian distribution for the wave facet inclinations, with a standard deviation determined by the horizontally isotropic wind speed (Cox & Munk 1954). Unless stated otherwise, we use a wind speed of 7 m s⁻¹ which corresponds to a moderate breeze⁴ and is a reasonable value for above the Earth’s oceans where the long-term annual mean wind speed varies approximately between 5 and 9 m s⁻¹ (Mishchenko & Travis 1997; Hsiung 1986). The reflection and transmission matrices for polarized light for a flat wave facet follow from the Fresnel equations and a proper treatment of the changing solid angle of the light beam when traveling from one medium to another (Garcia 2012a,b; Zhai et al. 2012, 2010). We use the shadowing function of Smith (1967) and Sancer (1969) to account for the energy budget across the rough air-water interface at grazing angles caused by neglecting wave shadowing (see also Tsang et al. 1985; Mishchenko & Travis 1997). The interface reflection and transmission matrices are corrected for the energy deficiency caused by neglecting reflections of light between different wave facets (Nakajima 1983). The elements of our matrix describing light reflected by the air–water interface have been verified against the code of Mishchenko & Travis (1997), and the energy balances across the interface (i.e., the sum of the reflection and transmission, for both illumination from below and above) are compared to Fig. 4 of Nakajima (1983). The elements of our reflection and transmission matrices have been verified against those of Zhai et al. (2010).

The subinterface ocean is modeled as a stack of horizontally homogeneous layers bounded below by a black surface. The reflection matrix of the subinterface ocean is computed with the adding–doubling algorithm (de Haan et al. 1987). Scattering of light in pure seawater is caused by small density fluctuations. However, the matrix describing the single scattering in this water may be approximated by that for anisotropic Rayleigh scattering (Hansen & Travis 1974), using a depolarization factor of 0.09 (Morel 1974; Chowdhary et al. 2006). The single scattering albedo of the water is computed according to $${\omega }_{\text{w}}(\lambda )=\frac{b_{\text{w}}(\lambda )}{b_{\text{w}}(\lambda )+a_{\text{w}}(\lambda )},$$(4)

where b_w and a_w are the wavelength dependent scattering and absorption coefficients for pure seawater, respectively, which are both expressed in m⁻¹. For b_w, we use the values tabulated by Smith & Baker (1981) and for a_w the values ofPope & Fry (1997)⁵. The ocean optical thickness is computed from the beam attenuation coefficient c(λ) = b_w(λ) + a_w(λ) and the ocean’s geometric depth. The ocean water does not include hydrosols (such as phytoplankton, detritus and bubbles). This would require a determination of the hydrosol scattering matrix elements (see, e.g., Chowdhary et al. 2006), which is beyond the scope of this study. We compared our wavelength dependent subinterface ocean albedo at a local solar zenith angle θ₀ of 30° with the bio-optical model for the so-called Case 1 Waters of Morel & Maritorena (2001) and find our ocean albedos, at the wavelengths employed in this paper, correspond to those of water with phytoplankton chlorophyll concentrations between 0 and 0.1 mg m⁻³.

The reflection by the ocean system, that is, by both the air–water interface and the subinterface ocean, is computed with the adding equations described by de Haan et al. (1987) and by assuming an infinitely thin interface (see also Xu et al. 2016). Light traveling from the atmosphere into the water is refracted into a sharp cone, resulting in a lower angular resolution in the cone compared to the resolution in the atmosphere. In order to accurately cover the sharp cone, we use rectangular super-matrices (see de Haan et al. 1987, for the use of the super-matrices) for the interface in the adding-doubling algorithm, allowing us to employ more and different Gaussian quadrature points inside the ocean layers than inside the atmosphere layers (see also Sect. 3.2.4 by Chowdhary 1999).

Finally, we adapt the hence obtained “clean” ocean reflection for the reflection by wind-generated foam, or white caps, using a weighted sum of the clean ocean reflection matrix and the Lambertian (i.e., isotopic and nonpolarizing) reflection matrix by foam. As a baseline, we assume an effective foam albedo of 0.22 (Koepke 1984) and the empirical relation of Monahan & Ó Muircheartaigh (1980) for the wind speed dependent fraction of white caps.

In Appendix A, we provide the governing equations that we use to compute the ocean reflection of polarized light as described in this section. The ocean parameters are listed in Table 2. The real refractive index of air is set equal to 1.0 in the air–water interface computations, while in reality, it will vary slightly with the wavelength (this variation is indeed taken into account in the computation of the scattering optical thickness of the atmospheric layers, see Sect. 2.3.1). The refractive index of water is assumed to be wavelength independent and equal to 1.33 (Hale & Querry 1973). Varying this refractive index between 1.34 (representative for liquid water and short wavelengths) and 1.31 (representative for water ice) yielded negligible differences in the phase angle dependent, disk-integrated, reflected fluxes and degree of polarization.

Fig. 1 Disk-integrated reflected total flux F (left column), linearly polarized flux Q (middle column) and degree of polarization P (right column) of an ocean planet with a rough, Fresnel reflecting surface (top row) and a planet with a black surface (bottom row) as functions of the phase angle, for various wavelengths. Both planets have a pure gaseous atmosphere and are cloud-free. The wind speed above the ocean is 7 m s−1, and the surface coverage of white caps 0.28%. The black dashed line in the ocean planet plots (top row) pertains to a planet covered by the rough Fresnel surface, with a black surface below, thus without sub-surface ocean, and without atmosphere above.

3 Results

Here, we present the results of our computations of the flux and polarization phase curves of ocean planets, starting with horizontally homogeneous, cloud-free planets (Sect. 3.1), followed by planets with different cloud coverage fractions (Sect. 3.2).

3.1 Horizontally homogeneous, cloud-free planets

Figure 1 shows the phase curves of an “ocean planet”, that is, a planet with a rough, Fresnel reflecting interface on top of a 100 m deep ocean that is bounded below by a black surface, and a “black surface planet”, that is, a planet with a black surface. The wind speed is 7 ms⁻¹, and the percentage of white caps thus 0.28% (see Monahan & Ó Muircheartaigh 1980). The planets have the same gaseous atmosphere, with a 1 bar surface pressure. The curves are shown for λ ranging from 350 to 865 nm. Because these model planets are symmetric with respect to the reference plane, their polarized flux U (see Eq. (1)) is zero at all phase angles. Also included in Fig. 1, are the phase curves of a planet with the same rough Fresnel reflecting interface, but without an atmosphere and without an ocean (on this planet, the Fresnel interface is thus bounded below by a black surface). While such a planet is physically impossible (without atmosphere, there would be noliquid water), its curves are included for comparison. Because we assume a wavelength independent refractive index for the Fresnel interface and because this planet has no atmosphere, its phase curves are wavelength independent.

The phase curves for the black surface planet (with atmosphere) are as expected (see, e.g., Buenzli & Schmid 2009; Stam 2008): with increasing λ, the contribution of the gaseous atmosphere decreases because of the ~ λ⁻⁴ dependence of the Rayleigh scattering cross-section, hence fluxes F and |Q| decrease with λ at all phase angles. Degree of polarization P increases with λ, because of the decreasing amount of multiple scattering, which usually lowers P, and the nonreflecting, black surface. Themaximum P occurs around α = 90° for all λ (α where P is maximum increases slightly withλ). With increasing surface albedo, the location of this maximum would shift to much larger phase angles while the maximum P would decrease (see, e.g., Fig. 4 of Stam 2008).

Up to α ~ 80°, the flux phase curves for the ocean planet are very similar to those of the black surface planet, although the ocean planet is somewhat brighter, especially at shorter wavelengths, where the ocean is a more efficient reflector. In particular, at 443 nm, the ocean planet’s geometric albedo is more than 20% larger than that of the black surface planet. Retrieving such small geometric albedo differences from observational data would require knowing the planetary radius to within a few percent.

With increasing α, the differences between the black surface and the ocean planet increase. To start with, at the longer wavelengths, where the atmosphere scatters less, the Fresnel reflection significantly brightens up the planet. This can also be seen from the phase curves of the planet with the Fresnel interface but without the atmosphere. Also, the ocean planet’s flux curves for the different wavelengths cross each other at α ~ 134°. Indeed, while at smaller phase angles, the ocean planet is blue, around 134° it would appear white. With further increasing α, the ocean planet would turn red, before darkening completely when α > 160°. The ocean planet’s polarized flux Q shows a similar crossing of the phase curves, and hence a color change from blue, through white, to red in polarization, except the crossing happens at a smaller phase angle: at about 123°.

This color change of the total and polarized fluxes with phase angle could be used to identify an ocean on an exoplanet, as the black surface planet does not show this color change with phase angle. It should be noted that a color change in the total reflected flux F only, thus without a detection of a color change in the polarized flux, would be ambiguous because a color change in the total flux can also occur in the presence of clouds (see Sect. 3.2) or for planets without an ocean and with a nonzero surface albedo⁶. Because the black surface planet does not show a crossing in Q in the presence of clouds (Fig. 6), a color change in Q would indicate an exo-ocean. Zugger et al. (2010) emphasize that it will be difficult to distinguish between a liquid and a frozen water surface when measuring total flux F, because of the small difference in refractive index (see Sect. 2.3.2). Measuring a crossing in Q, however, would allow such a distinction to be made, because Fresnel reflection by an icy surface strong enough to allow confusion with a liquid watersurface will only occur when the ice is clean, that is, not covered by rocks, snow or other rough materials.

The color change can be attributed to the wavelength dependence of the molecular scattering cross-section: with increasing α, the average atmospheric gaseous optical path-length encountered by the light that is received by the distant observer increases. While the light with short wavelengths (blue) is scattered in the higher atmospheric layers, the light with the long wavelengths (red) will reach the surface and is reflected by the ocean, back toward space and the observer. Without the gaseous atmosphere, all colors of incident light would be reflected by the ocean at the larger phase angles, as can be seen from the curves of the planet without atmosphere but with a Fresnel reflecting interface. Such an atmosphere-free planet would thus appear white at the largest phase angles. The influence of the atmospheric optical thickness (i.e., the surface pressure) on the color change of a planet will be discussed in Sect. 3.1.1.

We assume the incident starlight to be white: a strong wavelength dependent stellar spectrum would influence the planet’s actual color in total and polarized flux. On the other hand, in the case of real observations, the planet’s fluxes can be normalized to the spectrum of the incident starlight to retrieve the normalized planetary color shown in our figures.

Up to α ~ 80°, the phase curves of the degree of polarization P for the ocean planet are very similar to those of the black surface planet, although the ocean planet’s P is slightly lower at the same wavelength and phase angle. At the shorter wavelengths, the maximum P of the ocean planet is smaller than that of the black surface planet due to the increased amount of multiple scattered, thus low polarized light in the ocean planet’s reflected signal. With increasing λ, the influence of the atmosphere on the planet signal decreases and that of the Fresnel reflecting surface increases, shifting the maximum P to larger phase angles. The maximum P shifts toward twice the Brewster angle, thus toward 2arctann₂∕n₁ = 106° for our water ocean. This is indeed thephase angle of maximum P for the ocean planets without an atmosphere, as can be seen in Fig. 1, and as was also shown by Zugger et al. (2010).

The degree of polarization P is of course a relative measure and usually not associated with colors, although the color of polarization has also been used by, for example, Sterzik et al. (2019) in their analysis of polarized Earth-shine observations. Translating the phase curves of P as shown in Fig. 1 into colors, both the black surface and the ocean planet would be white at smaller phase angles, and red at intermediate angles. With increasing phase angle, the ocean planet would remain red in P, while the black surface planet would return to its white color.

3.1.1 Influence of the surface pressure p_s

To have a better look at the color changing ocean planets, Fig. 2 shows the phase curves of F and Q reflected by the ocean planet and the black surface planet for surface pressures ranging from 0.5 to 10 bar. The curves are shown in a Red-Green-Blue (RGB) color-scheme, in which the fluxes at various wavelengths are combined to mimic the actual perceived color. Figure 3 shows a few cross-sections of Fig. 2 for F, Q, and P (the latter is not shown in Fig. 2) for surface pressures of 0.5, 5.0, and 10.0 bar. We note that the boiling temperatures of water at these pressures are, respectively, 354 K (81°C), 425 K (152°C), and 452 K (179°C).

Figure 2 clearly shows that a cloud-free, black surface planet will appear mostly blue in F and Q across this surface pressure range. Only at the largest phase angles and surface pressures, it will appear white. The figure also shows the color change of an ocean planet from blue, through white, to red, and back to white, with increasing α, in both F and Q.

The difference between the total fluxes F reflected by the black surface planet and the ocean planet starts to disappear with increasing surface pressure, because with increasing pressure, the contribution of the surface reflection to the planetary signal decreases, as the scattering within the atmosphere increases and less light reaches the surface and subsequently the top of the atmosphere and the observer. The color change of the ocean planet thus also starts to disappear: for surface pressure higher than about 8 bar, the ocean planet will no longer be obviously red in any phase angle range and the presence of liquid surface water would not be recognizable in F. For surface pressures decreasing below 0.5 bar (not shown), the planet will become whitish at all phase angles because the influence of the atmosphere will disappear completely, while the Fresnel reflection is wavelength independent.

Interestingly, in polarized flux Q, the ocean planet’s color change from blue, through white, to red, and back to white with increasing α, remains strong, and is even still present at a surface pressure p_s of 10 bar. At such high surface pressures, the best approach to detecting an ocean appears to be using a combination of a short wavelength (λ ~ 550 nm) and a long wavelength (λ ~ 865 nm) filter, because the strength of the color change (i.e., not the relative but the absolute polarized fluxes) decreases with increasing p_s. This can be seen in Fig. 3: with increasing pressure, Q at intermediate phase angles increases for all λ. At phase angles beyond the crossing, Q decreases because of the increase of the atmospheric path length, with Q at longest (reddest) wavelengths remaining the largest, because the red light is still more likely to reach the surface than the blue light. The phase angle where the crossing takes place is approximately conserved for surface pressures above 2 bar (see upper right plot in Fig. 2). In addition, the phase angle range across which the ocean planet would be red in Q is much wider than the range across which the planet would be red in F. Indeed, with the color transition in Q occurring around α = 123°, the angular distance between the planet and the star would be larger than at α = 134°, where the color transition in F occurs, that is, 0.83″ versus 0.72″ for a planet orbiting its star at 1 AU, with a 1 pc distance between the observer and the star. We would like to point out that Q of the black surface planet also shows a color change from blue to red and back to blue/white at large phase angles, but the actual values of Q are very small here, for all wavelengths, as can also be seen in Fig. 1.

Zugger et al. (2010) suggest that the shift of the maximum P toward larger α due to the Fresnel reflection (see Fig. 1) could be used for identifying an exo-ocean. However, as can be seen in Fig. 3, this shift of the maximum P will depend on the surface pressure. Indeed, with increasing pressure, the maximum P of an ocean planet shifts toward α = 90°. Thus, a shift of the maximum P toward larger α’s would reveal the presence of an exo-ocean and indicate, necessarily, a small surface pressure p_s, as with a high surface pressure, there would be no observable phase angle shift of the maximum P.

Fig. 2 RGB colors of the (disk-integrated) reflected F (left column) and |Q| (right column) of the cloud-free ocean planet (top row) and the cloud-free black surface planet (bottom row) for surface pressures ranging from 0.5 to 10 bar in steps of 0.5 bar. Like in Fig. 1, the wind speed is 7 m s−1. The RGB colors are computed using weighted additive color mixing of the fluxes at λ = 443 nm (blue), 550 nm (green), and 670 nm (red), such that when these fluxes are equal, the color is white.

3.1.2 Influence of the wind speed v

Figure 4 shows the phase curves of the ocean planet for various wavelengths and three wind speeds v: 1, 7 (our baseline wind speed), and 13 m s⁻¹. Larger windspeeds create higher waves and thus a wider glint pattern on the planet (see also the disk-resolved planets in Fig. 7) and more white caps (see Monahan & Ó Muircheartaigh 1980). Higher waves, however, also increase extent of shadows.

First, as can be seen in the figure, for λ = 350 nm, the curves for F and Q do not show the glint for any v because of the large atmospheric gas optical thickness. For λ = 865 nm, the glint feature is the sharpest because at long wavelengths, this optical thickness is very small and thus most of the incident sunlight will reach the surface and, after being reflected, space without being scattered in the atmosphere. The wind speed does not affect the phase curves of F and Q, except slightly at the largest wavelengths, up to at least the crossing phase angles, where the curves for the different λ’s coincide, that is, α = 134° for F and α = 123° for Q. Indeed, the crossing phase angles themselves are virtually independent of v. At phase angles beyond the crossing, F and Q do depend on v: the strength (amplitude) of the glint in F and Q increases with decreasing v at every λ. The main reason is that the smaller v, the smaller the variation in wave angles, and hence the less diffuse the reflection, and the contribution of shadows will be smaller.

Figure 4 also shows P of the ocean planet. For λ > 350 nm, the maximum P decreases slightly with increasing wind speed v. With increasing λ, the locationof the maximum also shifts to slightly larger phase angles (i.e., from 100° at 865 nm and 1 m s⁻¹, to 108° at 865 nm and 13 m s⁻¹). The decrease and shift of the maximum P are not directly due to the increase of the surface roughness with increasing v. Indeed, in the absence of foam and without an atmosphere (or with an optically thin atmosphere), the phase curve of P is actually independentof v (see Appendix A) because the probability distribution function describing the wave facet inclinations (Cox & Munk 1954, 1956) influences F and Q in the same way, thus leaving P undisturbed. In the presence of an atmosphere, however, light that has been reflected by the Fresnel interface can subsequently be scattered in the atmosphere, and such scattering processes can influence F and Q differently, and thereby influence P.

Another factor that influences P and in particular the decrease and shift of the maximum P with increasing v is the increasing surface coverage by reflecting foam. Indeed, the white cap percentages are 0.0003, 0.28 and 2.46% for 1, 7, and 13 m s⁻¹, respectively. An increase in white caps slightly increases F (although in Fig. 4, the influence of the ocean roughness on F is much larger at the large phase angles), and only slightly decreases Q because the white cap reflection is assumed to be Lambertian. In Appendix B, we provide some more results, such as for different white cap albedos.

Fig. 3 Similar to Fig. 1, but for the ocean planet, and for surface pressures equal to 0.5 bar (top row), 5 bar (middle row), and 10 bar (bottom row).

Fig. 4 Similar to Fig. 1, but for the ocean planet and different wind speeds. The wind speeds and associated white cap percentages are: 1 m s−1 and 0.0003%(solid lines), 7 m s−1 and 0.28% (dashed lines) (cf. Fig. 1), 13 m s−1 and 2.46% (dotted lines).

Fig. 5 Similar to Fig. 1, but for the ocean planet (solid lines) and a planet with the same gaseous atmosphere and the same rough Fresnel interface, but without the water body underneath (dashed lines). The wind speed v is 7 m s−1.

3.1.3 Influence of the ocean color

In Fig. 1, we showed that the flux F that is reflected by the ocean planet is larger than that reflected by the black surface planet at the same phase angle, in particular at small phase angles. As described in Sect. 2.3, the water body below the Fresnel interface is not black (the surface below the water body is). Here, we will show the influence of the water body’s color on the reflected signals.

Figure 5 shows the phase curves of the ocean planet with the reflection of the water body below the Fresnel interface described by Rayleigh scattering, together with the curves for a planet with the same, cloud-free atmosphere and the same Fresnel interface, but with a black surface below the interface (thus without the water body). The difference between the curves shows how the scattering within the water increases F and |Q| at the shorter wavelengths: the largest effect is seen for λ = 443 nm, while there is virtually no effect for λ > 550 nm. Concluding, the blueness of the ocean is a main contributor to the increase in F at small phase angles and for λ < 550, as shown in Fig. 1.

In F, the influence of the water body decreases with α: a global ocean makes a planet slightly more blue for α < 90°. In polarized flux |Q|, the influence is largest for 60° < α < 80°. The scattering in the water decreases P somewhat at the shorter wavelengths and intermediate phase angles. Zugger et al. (2010) mention that the maximum P is limited by scattering within the water. They, however, model the reflection by the water body below the interface as a Lambertian surface (which scatters isotropically and unpolarized), and their results cannot be directly compared against ours that do include (polarized) scattering within the water body.

3.2 Cloudy planets

The planets in the previous sections were cloud-free. A planet with an ocean on its surface, will, however, very likely also have clouds. Here, we will discuss the influence of the clouds on the detectability of the ocean, starting with the influence of a horizontally homogeneous cloud deck (Sect. 3.2.1), followed by the influence of broken clouds (Sect. 3.2.2).

3.2.1 Influence of horizontally homogeneous clouds

Figure 6 shows the phase curves of an ocean planet that is completely covered by a horizontally homogeneous cloud deck for different wavelengths and wind speeds v (white caps are included). The details on the cloud deck parameters can be found in Table 1. The curves for a completely cloudy planet with a black surface (a_s = 0) are also shown.

At all wavelengths, F of the cloudy ocean planet is somewhat higher than that of the cloudy black surface planet, thanks to the added reflection by the Fresnel interface and the ocean body. As can be seen, F of the cloudy ocean planet is virtually independent of v at all phase angles. For α up to about 134°, this was to be expected, as F of the cloud-free ocean planet was also virtually independent of v (see Fig. 4). For the larger phase angles, clouds apparently smooth out the flux differences that were due to the wind on the cloud-free ocean planet.

While for the cloud-free ocean planets, the phase curves of F crossed at α = 134° (see Fig. 4), those for the cloudy ocean and the cloudy black surface planet cross at α = 108°. Analogous to the reflection by the cloud-free ocean planets (Sect. 3.1), the crossing is due to the combination of Rayleigh scattering above the clouds and the reflection by the clouds, with more light with longer wavelengths reaching the clouds and subsequently, after reflection by the clouds, the top of the atmosphere than light with shorter wavelengths. A completely cloudy ocean or black surface planet will thus change color from blue, through white, to red with increasing α. With increasing cloud optical thickness and at α’s below the crossing, a planet’s blueness decreases, as the contribution of light that has been scattered by cloud particles increases. With increasing cloud top altitude (i.e., decreasing cloud top pressure) and at the largest α’s, a planet’s redness decreases, as more blue light will reach the cloud and the top of the atmosphere after reflection.

The bumps in the flux phase curves around α = 40° are due to the primary rainbow: light that has been reflected once within the spherical water cloud droplets (Hansen & Travis 1974; Karalidi et al. 2012a; Bailey 2007). There is only a very slight wavelength dependence in the position of this rainbow, it will thus be mostly white. Indeed, as pointed out earlier by Karalidi et al. (2011), the colors of this rainbow due to scattering by the small cloud droplets (the blue rainbow occurs at a larger α than the red rainbow) are reverse from those of the rainbow in rain droplets (red at a larger α than blue). The phase curves of P and Q also show the primary rainbow feature and indeed, at α ~ 56°, also a small bump for the secondary rainbow. The P and Q phase curves for the completely cloudy ocean planet are very similar to those of the cloudy black surface planet. Also, Q shows no color crossing point, and there is no shift of the maximum value of P toward twice the Brewster angle, such as for a cloud-free ocean planet (see Figs. 1–4). Our model clouds have an optical thickness of about 5.0 at 550 nm (see Table 1), which is a normal value for Earth-like clouds. With optically thin clouds, the ocean signal will be stronger. However, we can conclude that it will be impossible to detect an exo-ocean when it is covered by a homogeneous, optically thick cloud deck.

Fig. 6 Similar to Fig. 1, but for the ocean planet with a horizontally homogeneous cloud deck. The wind speed is 1 m s−1 (solid line), 7 m s−1 (dashed line), or 13 m s−1 (dashed-dotted line). The dotted lines are the phase curves of a planet with a black surface and the same cloud deck.

3.2.2 Influence of an inhomogeneous cloud deck

With patchy clouds, both the location and the cloud coverage fraction influences the reflected light signal (assuming for simplicity that the clouds are all the same and at the same altitude). Here, we show the influence of a horizontally inhomogeneous cloud deck on F, Q, and P of light reflected by an ocean planet. Like on Earth, our model cloud deck consists of patches of clouds, as described in Sect. 2.3.1 with a cloud coverage fraction f_c. Figure 7 shows the disk-resolved RGB colors of the reflected F and |Q| for model planets with different cloud coverage fractions f_c and wind speeds v (for every combinationof f_c and v, the model planet has a different cloud pattern) at α = 80°. The glint can be seen to appear between the clouds, and with increasing v, the width of theglint pattern increases while its peak value decreases.

Figures 8–10 show F, Q, and P for the ocean planet for different values of f_c and v. Also included arethe curves for a black surface planet beneath the clouds. We used 160 pixels across the equator, just like in Fig. 7. For every f_c, we computed F, Q, and U (since planets with patchy clouds are usually not symmetric with respect to the reference plane) for 300 different patchy cloud patterns. The figures show the averages and the 1-σ standard deviation (the variability). Despite the nonsymmetric cloud patterns, the disk-integrated linearly polarized flux U appears to be smaller than 0.00003 for all model planets, across all phase angles. We thus do not show U in the figures (even though it is very small, U is included in the computation of P).

As can be seen in Fig. 8, total flux F generally increases with increasing f_c, and the curves for the ocean planet are somewhat higher than for the black surface planet. The influence of the variability of the cloud coverage for a given value of f_c is small and only apparent at the longest wavelengths and the largest α’s. At the smallest wavelengths, the significant Rayleigh scattering above and between the clouds subdues the contrast between the bright clouds and the dark background, thus reducing the variability. At the largest phase angles, the surface area on the disk that provides the signal is smallest, and there is more variability in the coverage of cloudy pixels across that area (f_c is defined on the whole disk). A more in-depth discussion of the relation between α and variability can be found in Rossi & Stam (2017).

Polarized flux |Q| generally decreases with increasing f_c (Fig. 9), and |Q| is slightly larger for the ocean planet than for the black surface planet. The curves for the black surface planet show virtually no variability, while the variability in Q of the ocean planet is significant for the longer wavelengths and α’s larger than 90°, in particular for small wind speeds.

Like the F and Q phase curves of the cloud-free and completely cloudy ocean planets, the F and Q phase curves of the ocean planets with patchy clouds for the different λ’s cross each other at a given phase angle α (see Figs. 8 and 9). Around that α (which is different for F and Q), the planets thus appear white in total and polarized fluxes. The crossing is also present in F for large cloud fractions (f_c ≥ 0.5) over a black surface planet, in agreement with the earlier results for completely cloudy black surface planets (Fig. 6). A black (or dark) surface planet with patchy clouds will, however, appear more whitish than red at the largest phase angles, because the spectral dispersion of F is small.

Figure 10 shows the P phase curves of the ocean and black surface planets with patchy clouds. The largest variation in P is seen at the longest wavelengths, where the contribution of the gaseous atmosphere is smallest and the contribution of the surface is largest, and at intermediate to large α’s, where the contribution of the glint is strongest, and thus where the location of a cloud is relatively more important than at the smaller α’s (cf. Fig. 7). There appears to be little dependence of both the average curves and the variability in P on wind speed v, except for α > 100°.

For ocean planets with patchy clouds, the curves of P show crossings, that is, the planet will show a color change, just like the curves for F and Q. While this crossing in P does not occur for the completely cloudy planet (f_c = 1.0), as was shown in Fig. 6, it does show up for f_c = 0.75 (see Fig. 10). The black surface planet with patchy clouds does not show a color change. Measuring a color change in P in the relevant phase angle range would thus indicate the presence of an ocean, even in the presence of patchy clouds.

Apart from the crossing in the P phase curvefor a cloudy ocean planet, there are more differences from the phase curve of the cloudy black surface planet. For example, Fig. 10 shows that for λ = 865 nm and both f_c up to (at least) 0.50, the slope of P for the ocean planet is positive between 80° < α < 100°, while that for the black surface planet is negative. Additionally, the variability in P is much larger for the ocean planet than for the black surface planet because of the appearance and disappearance of the polarized glint between the clouds.

Figure 11 shows the disk-resolved F and |Q| of the ocean planet and the black surface planet with patchy clouds (f_c = 0.50) at a range of phase angles α to illustrate the reddening of both planets in F and the reddening of the ocean planet in |Q|. It can also be seen that the glint makes the ocean planet somewhat lighter in blue than the black surface planet. In the images of the black surface planet at the largest phase angles the reddening of the planet in F due to the clouds is visible. This cloud reddening also happens with the ocean planet, but there the reddening of the glint is more significant.

Fig. 7 Disk-resolved (160 × 160 pixels) model output RGB color images in F (Cols. 1–3) and |Q| (Cols. 4–6)of ocean planets with patchy clouds. The brightness of a pixel is scaled with its average spectral F. Phase angle α = 80°. The cloud coverage fraction fc is 0.25 (Cols. 1 and 4), 0.5 (Cols. 2 and 5), or 0.75 (Cols. 3 and 6), and the wind speed v is 1 m s−1 (top row), 7 m s−1 (middle row), or 13 m s−1 (bottom row). For every combination of fc and v, the cloud pattern across the disk is different.

4 Discussion of the color change

Our results show that the presence of an exo-ocean could be revealed through a color change in the polarized flux Q of starlight that is reflected by a planet, except when the planet is completely covered by a thick cloud deck. The total reflected flux F and degree of polarization P can also show color changes. The crossing α’s, or α_cc ’s, where the color change takes place, are generally different for F, Q, and P, and depend on the cloud fraction f_c.

An indication of such a crossing in measured phase curves of P of sunlight that is reflected by the Earth, has recently been presented by Sterzik et al. (2019). They measured P for α ranging from ~50° to ~ 140° as derived from Earth-shine data obtained with the Very Large Telescope (VLT) at the European Southern Observatory (ESO) in Chile. The phase angle range of these Earth-shine data is limited by the lunar phase during the observations⁷, so unfortunately, the phase curves do not include the rainbow angle, which would have provided an interesting comparison with the strength of the rainbow in our simulations. In the Earth-shine measurements, α_cc appears to vary from about 120° to 135°. The observations span a period of about two years and during the observations either South America and the Atlantic Ocean (east of Chile) were turned toward the Moon or the Pacific Ocean (west of Chile) (Sterzik et al. 2019). The α_cc of about 135° was measured over the Pacific Ocean.

Figure 12 shows RGB color phase curves of F and Q reflected bythe ocean planet and the black surface planet for cloud fractions f_c ranging from 0.0 to 1.0 in steps of 0.01. The phase curves for each wavelength are computed using so-called quasi horizontally inhomogeneous ocean planets, in the weighted sum approach as described by Stam (2008). In this approach, F and Q of a horizontally inhomogeneous planet are computed by taking a weighted sum of the F’s and Q’s computed for horizontally homogeneous planets with the atmosphere-surface combinations that occur on the inhomogeneous planet. Each weighting factor equals the fraction of the planetary disk that is covered by a given atmosphere-surface combination. A careful comparison (see Appendix C) shows that while the weighted sum approach does not directly provide variances due to the spatial distribution of cloud patches across the disk, it does indeed allow us to rapidly and accurately compute the average phase curves of F, Q, and hence the respective values of α_cc’s, using a small step size for f_c. Figure 12 illustrates that the color change in Q (i.e., where the RGB color changes from blue, through white, to red) is unique for ocean planets for f_c up to 0.98. For f_c > 0.98, the crossing of the phase curves at the RGB wavelengths is too small to detect.

Figure 13 shows α_cc for F, Q, and P as functions of the cloud fraction f_c, using the weighted sum approach. The P of the inhomogeneous planet is computed from the summed total F and Q (U of such a quasi horizontally inhomogeneous planet is always zero). The figure also includes α_cc ’s for the ocean planets as shown in Figs. 8–10, and thus computed for actual horizontally inhomogeneous planets. A comparison between these data points (the symbols) and the data computed assuming quasi horizontally inhomogeneous planets (the lines) in Fig. 13, shows that the weighted sum approach predicts α_cc with a maximum deviation of 0.35° (in P for f_c = 0.50). Also shown in Fig. 13, is α_cc computed for higher clouds (with a cloud top pressure, p_ct, of 0.5 bar) and optically thicker clouds (with b^a = 10). We did not include curves for planets with black surfaces in Fig. 13, because there are no color changes in Q and P for such planets (cf. the bottom rows of Figs. 9 and 10), while the color change in F is not very strong due to the small dispersion. The horizontal axis at the top indicates the minimum inclination angle orbital i^* that is required to be able to observe an exoplanet at that particular phase angle; that is, a color change expected to take place at α_cc = 140° can be observed if the inclination angle of the planetary orbit is at least 50°.

Increasing the cloud coverage f_c on an ocean planet decreases α_cc for the total flux F, as can be seen in the figure, because more long wavelength light will be reflected at a given α: α_cc ranges from ~ 134° for cloud-free planets (f_c = 0.0) to ~ 108° for a fully cloudy planet (f_c = 1.0). Increasing the cloud top altitude for a given f_c, increases α_cc, because the higher the clouds, the smaller the contribution of short wavelength (Rayleigh scattered) light at a given α. An optically thicker cloud reflects more light, and will thus decrease α_cc in F.

Increasing f_c increases α_cc for the polarized flux Q, because more light will be reflected by the clouds, and less by the ocean at a given α: α_cc ranges from ~ 123° for cloud-free planets (f_c = 0.0) to ~157°, where the color change starts to disappear for an almost fully cloudy planet (f_c = 0.98). Increasing the cloud top altitude (i.e., decreasing p_ct) at a given f_c, increases α_cc because with higher clouds, a longer atmospheric path length, thus a more slanting incoming direction, is required for the short wavelength light to be sufficiently scattered out of the incident beam. Increasing the cloud optical thickness, b^a, has virtually no effect on α_cc of Q, because a thicker cloud will add mostly unpolarized light to the reflected signal.

For the degree of polarization P, Fig. 13 shows that increasing f_c, increases α_cc. Indeed, the range of α_cc appears to be largest for P: for f_c as small as 0.03, α_cc ~ 88°, while for f_c = 0.98, the color change starts to disappear at α_cc ~ 161° (cf. Fig. 6). There is no crossing in P for f_c < 0.03, which is consistent with our conclusion that a crossing in P indicates an ocean that is partly covered by clouds: without clouds, there is no color change. Increasing the cloud top altitude, decreases α_cc for P, while increasing the cloud optical thickness b^a, increases α_cc slightly and only for f_c < 0.3.

The crossings in the phase curves of P derived from Earth-shine measurements by Sterzik et al. (2019) occurred between about 120° and 135°. Taking into account the order of magnitude of the variation in α_cc with the ± 1-σ standard deviation of the phase curve, the cloud coverage fraction would have ranged between 0.2 and 0.6. The later value pertains specifically to observations with the Pacific Ocean (west of Chile) turned toward the Moon. Observations taken when South America and the Atlantic Ocean were turned toward the Moon (east of Chile) would be influenced not only by the reflection by the ocean but also by South America. The influence of land surface reflection (albedo, reflection properties and coverage fraction) on α_cc will be part of a future study.

Finally, our model clouds consist of liquid water particles only, but on Earth, high altitude clouds can of course also contain ice particles. If these particles are flat and horizontally oriented, they will produce their own glint feature (examples from Earth-observations, are described in Marshak et al. 2017; Bréon & Dubrulle 2004; Chepfer 1999; Noel & Chepfer 2004). Whether or not the polarization signature of the glint due to such (usually optically thin) ice clouds will mimic and/or influence that of a surface ocean should be part of a future study. In particular, the influence of bright surface areas, such as high planetary latitudes that could be covered by ice and/or snow, on the disk–integrated fluxes and polarization would be interesting to investigate further. Cowan et al. (2012) postulate that the larger relative contributionof light reflected by polar regions at a planet’s crescent phase (assuming the light equator coincides more or less with the planet’s true equator), could yield a false positive for the glint of light reflected by an ocean. However, unless such snow and ice covered surfaces are very clean, we would mostly expect an increase in total flux F, not in Q and/or U, and thus a decrease of P. Such an increase in F and decrease in P would mainly happen at the longest wavelengths, and would thus not result in a color change in Q and/or P and a false-positive exo-ocean detection.

Fig. 8 Disk-integrated reflected F of ocean planets with patchy clouds for 5 wavelengths: 350 nm (pink), 443 nm (blue), 550 nm (green), 670 nm (red), and 865 nm (brown), and for wind speeds of 1 m s−1 (first row), 7 m s−1 (second row), and 13 m s−1 (third row), and cloud fractions fc of 0.25 (left column), 0.50 (middle column), and 0.75 (right column). Bottom row: corresponding curves forplanets with a black surface instead of an ocean. Each curve is the average of 300 curves computed for 300 different cloud patches configurations, which are different for each phase angle α. The 1-σ standard deviation is represented by the shaded areas around each curve (except for the longest wavelengths and largest phase angles, this deviation is very small).

Fig. 9 Similar to Fig. 8, but for the linearly polarized flux Q.

Fig. 10 Similar to Fig. 8, but for the degree of linear polarization P.

Fig. 11 Similar to Fig. 7 for ocean planets (top) and black surface planets (bottom) covered by patchy clouds. The cloud coverage fraction fc is 0.50, and the wind speed v on the ocean planet is 7 m s−1. Phase anglesα range from 0° (left), through 30°, 60°, 80°, 100°, 120°, and 140°, to 160° (right).

5 Summary and outlook

We have presented computed total and (linearly) polarized fluxes, F and Q, and the degree of (linear) polarization, P, of starlight that is reflected by ocean planets, that is, planets with a liquid ocean ruffled by the wind on their surfaces, with the aim of identifying observables for the presence of an ocean, by comparing their signals against those of so-called black surface planets, with the same atmosphere but with a black surface⁸. Below, we summarize and reflect on our results.

The total flux F of light reflected by an ocean planet will be slightly larger than that of a black surface planet, up till phase angles α of about 100°. Across these phase angles, both ocean and black surface planets will appear blue, for every cloud coverage fraction f_c. With α increasing from about 100°, the color of a black surface planet will change from blue to white, while that of an ocean planet will change from blue to white, then to red, and then back to white (but extremely dark) at the largest phase angles.

The crossing phase angle α_cc, that is, the α where the planet’s color changes from blue, through white, to red, increases with increasing surface pressure: at 1 bar, α_cc ~ 134°, while at 6 bar, α_cc ~ 140°. The phase angle where the cloud-free ocean planet changes back from red to white, decreases with increasing surface pressure p_s. For p_s increasing up from about 6 bars, the planet’s reddish phase fades and with that the difference between the ocean planet and the black surface planet starts to disappear.

Angle α_cc decreases with increasing f_c: while for a cloud–free planet (f_c = 0.0) and p_s ~ 1 bar, α_cc ~ 134°, for a fully cloudy planet (f_c = 1.0), α_cc ~ 108°. Angle α_cc is insensitive to the wind speed v over the ocean, although the strength of the color change does depend on v: the smaller v, and thus the lower the waves and the less whitecaps, the redder the ocean planet would appear.

The polarized flux Q of light reflected by an ocean planet shows a similar color change from blue, through white, to red with increasing α as the total flux F. Crucially, for Q, the color change does not occur for black surface planets, only for ocean planets, independent of the cloud coverage fraction f_c. Detecting a color change in Q would thus indicate the presence of an ocean on the surface of a planet. Angle α_cc for Q increases with increasing f_c. For f_c = 0.0, α_cc ~ 123°, while for almost completely cloudy planets (f_c = 0.98), the color change starts to disappear at α_cc ~ 157°. Also, the color change in Q remains strong at higher surface pressures, and is even still present (somewhat) at a surface pressure p_s of 10 bar. Furthermore, for p_s ≿ 2 bar, the crossing angle α_cc is more or less independentof both p_s. Angle α_cc also appears to be independent of v, although, like with F, the smaller v, the redder the ocean planet would appear in Q.

For ocean planets with patchy clouds (0.0 < f_c < 1.0), the variability of Q due to variation in the location of the clouds is, at each α, much larger than a black surface planet with a similar value of f_c, in particular at longer wavelengths. This variability is caused by the highly polarized glint appearing and disappearing from behind the patchy clouds.

The degree of polarization P of light reflected by a cloud-free ocean planet is slightly lower at λ = 350, 443, and 550 nm compared to that of a cloud-free black surface planet, mostly due to the larger flux F that is reflected by the ocean. Because the probability distribution of the slopes of the wave facets and the shadow function affect fluxes F and Q in the same way, the dependence of P on the wind speed v is only due to the white cap distribution and scattering within the atmosphere, and thus relatively small. Indeed, for v = 13 m s⁻¹, the maximum P is about 5% smaller than for v = 7 m s⁻¹ at λ = 865 nm and using a whitecap albedo of 0.22. The phase angle where P is maximum, increases from around 90° to almost twice the Brewster angle of 106° with increasing λ, as the optical thickness of the atmosphere decreases and the contribution of the ocean to the reflected signal becomes more dominant. The location of the maximum P, however, depends on the surface pressure p_s. Indeed, an ocean-induced shift of the phase angle where P is maximum toward larger values only occurs if p_s is small (if p_s is large, the presence of an ocean does not lead to a phase angle shift).

Like in Q, the variability in P of an ocean planet, especially at the longer wavelengths and for all wind speeds, due to the glint appearing and disappearing from behind the clouds, would be much larger than the variability in P of a black surface planet. A measurement of such variability could thus be used to identify an exo-ocean.

An ocean planet will change color in P, but only if there are clouds. The color change phase angle α_cc strongly increases with f_c, from about 88° for f_c = 0.03 to about 161° for f_c = 0.98. It should be noted that this range of α_cc is much broader than that for F or Q, and it includes phase angles with a favorite planet-star separation. For f_c < 0.03, an ocean planet does not change color in P. It would be interesting to investigate whether such low cloud fractions were to be expected in the presence of an ocean. The P phase curves as derived from Earth-shine measurements by Sterzik et al. (2019) show crossings that indicate cloud coverage fractions f_c between 0.2 and 0.6.

While these cloud coverage fractions appear to be realistic for the Earth, the Earth-shine measurements covered a time period of several years, and two regions on Earth: the Pacific ocean (west of Chile) and the Atlantic ocean and South America (east of Chile). In addition, although Sterzik et al. (2019) took careful measures to correct their observations for the depolarizing reflection by the lunar surface, the Earth’s polarization phase functions might still include some unknown, and possibly wavelength dependent residual of the reflection by the Moon, that could influence α_cc as derived from the curves. Indeed, flux and polarization measurements of sunlight that is reflected by the whole Earth, across the visible spectrum, covering our planet’s daily rotation and all phase angles, would allow the confirmation of using the color change of, in particular, Q for identifying the presence of an ocean. Such measurements (which would also show the rainbow feature that is not covered in the Earth-shine data) could be performed from the lunar surface or from a lunar orbiter (for measurement concepts, see Karalidi et al. 2012b; Hoeijmakers et al. 2016).

Our simulations only include scattering by the atmospheric gas and by clouds consisting of liquid water droplets. Scattering and absorption by aerosol particles has not been included. Over the ocean, most aerosol particles are expected to be sea salt particles. We note that most of these aerosol particles will reside in the boundary layer, below the clouds. The influence of aerosols on the color change of a planet would depend on the aerosol column number density, and the aerosol microphysical properties, such as size, shape, and composition. The observation of the color change in P in Earth-shine measurements by Sterzik et al. (2019) strongly suggests that the influence of aerosol particles will be minor, but a follow-up study with realistic ocean-type aerosols added to our atmosphere model, will give more insight into their role. In such a follow-up study, also the influence of oriented ice crystals can be investigated.

In our numerical simulations, we assume that the incident starlight is wavelength independent. This assumption is irrelevant for our results, as measured fluxes F and Q could straightforwardly be normalized to the color of the incident stellar light⁹. The degree of polarization P is of course independent of any wavelength dependence of the incident stellar spectrum¹⁰.

Our simulations provide yet another example of the value of polarimetry as a tool for exoplanet characterization: especially the detection of a color change of the polarized flux Q would uniquely identify a liquid surface ocean, and a color change in P would indicate an ocean that is at least partly covered by clouds. The phase angle α_cc at which the color change in Q takes place depends on the cloud coverage fraction f_c: α_cc ~ 123° for a cloud-free planet (f_c = 0.0) and α_cc ~ 157° for an almost fully clouded planet (f_c = 0.98) where the color change starts to disappear. The color change in P takes place at phase angles from about 88° for f_c = 0.03 to ~161° for f_c = 0.98. Whether or not an exoplanet ever attains the phase angle at which a color change takes place, depends on the inclination angle i of the planetary orbit, as can be seen in Fig. 13. A phase angle of 123° is attained along orbits with an inclination angle of at least 33°, while for a phase angle of 161°, the orbital inclination angle has to be at least 71°. For intermediate cloud fractions, such as that on Earth, inclination angles larger than about 40° are required.

Current and near future space missions such as NASA’s TESS (Ricker et al. 2015) and ESA’s PLATO mission (Rauer et al. 2014) will provide transiting exoplanets, thus with i close to 90°, around bright, nearby stars. These planets will be excellent targets for searching for a color change in the polarized light signal. This color change can be measured at smaller phase angles, where the angular separation between the star and the planet is larger, than the glint, that is, the increased reflected flux F of starlight on the ocean surface, that was proposed as exo-ocean identifier by for example, Williams & Gaidos (2008) and Robinson et al. (2010) (our Fig. 8 shows the glint in F around α = 160° for a wind speed of 1 ms⁻¹). Thus, using polarization to identify exo-oceans is not limited to crescent planetary phase angles.

We have shown that the color change in total flux F can also occur with a zero surface albedo planet, with clouds below a Rayleigh scattering gas layer, while color changes in Q and P only occur for the ocean planets. We have also shown that the lack of a color change in F could mistakenly be interpreted as a negative for an exo-ocean in case of high surface pressures. The color change in Q remains present at high surface pressures (see Fig. 2). Furthermore, as shown by Cowan et al. (2012), the relative contribution of high planetary latitudes, which are more likely to be covered by ice and/or snow, to the disk-integrated fluxes increases with increasing phase angle toward the crescent phase, creating a false-positive detection of the glint. However, unless the ice is perfectly clean, we would not expect this increased contribution to yield a color change in P and/or Q. Polarimetry would thus be a more robust tool against the false-positive identification of exo-oceans. Dedicated numerical simulations would be required to confirm this.

The crucial role that liquid water plays for life on Earth should spur feasibility studies into measuring the color change of an exoplanet in polarized flux and/or degree of polarization by future telescopes that are designed to measure starlight that is reflected by exoplanets. Apart from EPOL (Keller et al. 2010), the imaging polarimeter for EPICS, the Exoplanet Imaging Camera and Spectrograph for the European Extremely Large Telescope (Kasper et al. 2010), current space telescope concepts, such as LUVOIR (LUVOIR Team 2018) and HabEx (Gaudi et al. 2018) do not include polarimetry across the visible¹¹, and will thus not be able to measure the color changes.

While our numerical simulations pertain to the planet-signal only, and do not include any background starlight, the color change in polarized flux Q and/or degree of polarization P could also be searched for in the combined light of the star and the planet. Obviously, P of this combined light would be very small because of the added, mostly unpolarized stellar flux, but extracting this very small planet signal from the background stellar signal and noise could possibly be facilitated using the variation of P along the planetary orbit, and the knowledge that the color change would be limited to the planetary signal.

Fig. 12 Similar to Fig. 2, but for cloud fractions ranging from 0.0 to 1.0 in steps of 0.01. The surface pressure is 1 bar. The results are shown for the phase curves generated by the weighted sum of cloud-free (fc = 0.0) and fully cloudy (fc = 1.0) planets (see Appendix C). The purple, yellow and cyan colors at α > 140° result from the high RGB color sensitivity on the very small absolute polarized fluxes |Q| for the different wavelengths at these phase angles (cf. Fig. 9).

Fig. 13 Phase angle αcc where the phase curves of the total flux F (black), the (linearly) polarized flux Q (gray), and the degree of (linear) polarization P (brown) of an ocean planet for the different wavelengths cross each other, as functions of the cloud coverage fraction fc. The wind speed v is 7 m s−1, although αcc does not significantly depend on v. The symbols(squares, circles, diamonds) are the αcc’s shown before in Figs. 8–10 as computed for horizontally inhomogeneous planets, with the error bars representing the 1σ standard deviation of 300 iterations with different locations of the patchy clouds. The solid, dashed, and dotted lines have been computed using the weighted sum approach (see text), with the clouds at different altitudes (pct = 0.5 or 0.7 bar) and for different cloud optical thicknesses (ba = 5 or 10). The horizontal axis at the top indicates the minimum orbital inclination angle i* required to be able to observe the exoplanet at the corresponding αcc (an inclination angle of 90° indicates a planetary orbit that is “seen” edge-on).

Acknowledgements

The authors thank Johan de Haan, Jacek Chowdhary, and Mike Zugger for sharing their knowledge and insights on the topic of this paper, and the authors thank the anonymous referee for the constructive comments that have improved this paper.

Appendix A The reflection matrix of the ocean

Because we fully include the polarization of light, we describe scattering, reflection and transmission processes by 4 × 4 matrices¹² instead of by scalars. The reflection of light by a gas–liquid interface (in our case, an air–liquid water interface), with waves shaped by randomly oriented, flat facets, is described by the following matrix R_I (see Mishchenko & Travis 1997; Zhai et al. 2010) $$\begin{array}{lll}{R}_{\text{I}}(\mu ,{\mu }_0,\varphi -{\varphi }_0)\hfill & =\hfill & S(\mu ,{\mu }_0,\sigma )\frac{\pi }{{\mu }_0}\frac{\text{PDF}({\mu }_n,\sigma )}{4\mu {\mu }_n}\hfill \\ \hfill & \hfill & \times L(-{\xi }_2)R_{\text{F}}({\mu }_i,n_1,n_2)L(-{\xi }_1),\hfill \end{array}$$(A.1)

with R_F the matrix describing the Fresnel reflection by a flat wave facet, with n₁ and n₂ the refractive indices of air and water, respectively (see Zhai et al. 2010, Eq. (47)). In Eq. (A.1), μ = cos θ, μ₀ = cos θ₀, μ_n = cos θ_n, and μ_i = cos θ_i, with θ, θ₀, θ_n, and θ_i the angles between the local zenith direction and the direction of propagation of the reflected light, the incident light (not only the direct beam but also the diffuse incident light), the normal on the wave facet, and the Fresnel reflection angle, respectively.

The matrices L are rotation matrices that are used to rotate the reference plane for the incident light (through the local zenith direction and the direction of incidence) to the reference plane for the reflection by the wave facet (through the local normal on the facet, and the directions of incidence and reflection), and next to the reference plane of the reflected light (through the local zenith direction and the direction of reflection) (see Hovenier & van der Mee 1983), given by $$L(\xi )=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos 2\xi & \sin 2\xi & 0\\ 0& -\sin 2\xi & \cos 2\xi & 0\\ 0& 0& 0& 1\end{array}\right],$$(A.2)

with ξ the angle between the two reference planes, measured rotating in the anticlockwise direction from the old to the new reference plane when looking in the direction of propagation (thus clockwise for an observer that looks toward the planet) (ξ ≥ 0°). The first rotation in Eq. (A.1) is over angle − ξ₁ and the second rotation over angle − ξ₂.

Furthermore, in Eq. (A.1), PDF is the probability density function of the wave facet inclination angles, as follows (see Breón 1993; Zhai et al. 2010) $$\text{PDF}({\mu }_n,\sigma )=\frac{1}{\pi {\sigma }^2{\mu }_n^3}\text{exp}\left(-\frac{1-{\mu }_n^2}{{\sigma }^2{\mu }_n^2}\right).$$(A.3)

The variance of the slopes of the wave facets, σ, depends on the wind speed v (in m s⁻¹) and is given by (see Cox & Munk 1954) $${\sigma }^2=0.003+0.00512v.$$(A.4)

A larger wind speed v thus increases the variance in the slopes.

S in Eq. (A.1) describes the shadowing, that is, the blockage of incident light and reflected light due to inclined wave facets. For this function we use (see Smith 1967; Sancer 1969) $$S(\mu ,{\mu }_0,\sigma )=\frac{1}{1+\Lambda (\mu ,\sigma )+\Lambda ({\mu }_0,\sigma )},$$(A.5)

where $$\begin{array}{ll}\Lambda (\gamma ,\sigma )=\hfill & \frac{1}{2}\{\frac{\sigma }{\gamma }{\left[\frac{1-{\gamma }^2}{\pi }\right]}^{1/2}\hfill \\ \hfill & \text{exp}\left[-\frac{{\gamma }^2}{{\sigma }^2(1-{\gamma }^2)}\right]-\text{erfc}\left[\frac{\gamma }{\sigma \sqrt{1-{\gamma }^2}}\right]\}.\hfill \end{array}$$

Here, erfc is the complementary error function.

Part of the light that is incident on the Fresnel reflecting interface will be transmitted into the liquid ocean. Once in the ocean, this light can be scattered, or absorbed, or it can be absorbed by the black surface that we assume below the water body. The transmission of incident light from above (the star/sun or the sky) through the rough interface into the water, is described by the following transmission matrix: $$\begin{array}{lll}{T}_{\text{I}}(\mu ,{\mu }_0,\varphi -{\varphi }_0)\hfill & =\hfill & S(\mu ,{\mu }_0,\sigma )\frac{\pi }{{\mu }_0}\frac{\text{PDF}({\mu }_n,\sigma )}{\mu {\mu }_n}\hfill \\ \hfill & \hfill & \times L(-{\xi }_2)T_{\text{F}}({\mu }_i,n_1,n_2)L(-{\xi }_1)\hfill \\ \hfill & \hfill & \times \frac{n_2^2{\mu }_t{\mu }_i}{{(n_2{\mu }_t-n_1{\mu }_i)}^2},\hfill \end{array}$$(A.6)

where μ_t = cos θ_t, with θ_t the transmission angle at the wave facet, and T_F is the Fresnel transmission matrix for a flat interface. We note that T_F contains a factor (n₂∕n₁)(μ_t∕μ_i) to account for the changing solid angle of the light beam as it is transmitted through the air-water interface and for the changing flux because of the different permittivity of water compared to that of air (see Eq. (48) of Zhai et al. 2010). The interface reflection and transmission matrices for illumination from below, $R_{\text{I}}^{*}$ and $T_{\text{I}}^{*}$, respectively, may readily be derived following a similar approach (see Zhai et al. 2010).

The matrix describing the reflection by the water below the interface, R_W, is computed using the adding-doubling algorithm (see Sect. 2.3.2). The water body is bounded below by a black surface (the algorithm can also handle nonblack surfaces, for example, when modeling reflection by shallow waters above lightly colored sand).

We compute the reflection matrix, R_CO, of the “clean” ocean body as a whole, that is, of the rough interface and the water body, but without white caps, according to (Xu et al. 2016) $$\begin{array}{lll}\hfill & \hfill & Q_1=R_{\text{I}}^{*}{R}_{\text{W}},\hfill \\ \hfill & \hfill & Q_{p+1}=Q_1{Q}_p,\hfill \\ \hfill & \hfill & Q={\sum }_{p=1}^{\infty }{Q}_p,\hfill \\ \hfill & \hfill & D=T_{\text{I}}+Q{T}_{\text{I}},\hfill \\ \hfill & \hfill & U=R_{\text{W}}D,\hfill \\ \hfill & \hfill & R_{\text{WL}}=T_{\text{I}}^{*}U,\hfill \\ \hfill & \hfill & R_{\text{CO}}=R_{\text{WL}}+R_{\text{I}},\hfill \end{array}$$

where D, T_I, U and $T_{\text{I}}^{*}$ are rectangular super-matrices (see Sect. 2.3.2) and R_WL describes the reflection by the subinterface water body (including the reflections between the water body and the interface).

Finally, we computed R_O, the reflection matrix of the entire ocean (interface plus water body plus black surface below the water), and including white caps as follows (Koepke 1984; Sun & Lukashin 2013) $$R_{\text{O}}=q{R}_{\text{WC}}+(1-q)R_{\text{CO}},$$(A.14)

where R_WC is the matrix describing the Lambertian reflection by the white caps with a foam albedo a_foam and q is the fraction of white caps as determined by the empirical relation of Monahan & Ó Muircheartaigh (1980) $$q=2.95\times {10}^{-6}{v}^{3.52}.$$(A.15)

The influence of the white caps and the white cap albedo will be discussed in Appendix B.

Appendix B The reflection by white caps

Wind blowing across an ocean generates waves on the surface and with that, white caps or foam. The reflection matrix of the ocean as a whole (interface, water body, surface below the water) is computed using a weighted sum of the reflection matrix of the ocean without white caps, with its surface described by wave facets reflecting as Fresnel surfaces, and the reflection matrix of the white caps (see Eq. (A.14)).

Fig. B.1 Similar to Fig. 4, but for 865 nm and different wind speeds and whitecap properties. The speeds and the respective whitecap properties are 7 m s−1 and afoam = 0.22 (solid lines) (cf. Fig. 4), 13 m s−1 and afoam = 0.22 (dotted lines), 13 m s−1 and afoam = 0.33 (dashed lines), and 13 m s−1 but no whitecaps, i.e., q = 0 (dashed-dotted lines).

The reflection by the white caps or foam is assumed to be Lambertian, thus isotropic and nonpolarizing, and the surface that is covered by white caps is flat. As baseline, we assume an effective foam albedo, a_foam, of 0.22, as empirically derived by Koepke (1984). Figure B.1 shows the total and polarized fluxes and the degree of polarization of the ocean planet at 865 nm, for a wind speed v of 7 m s⁻¹ and our standard whitecap albedo a_foam of 0.22, for v = 13 m s⁻¹ and a_foam = 0.22, for v = 13 m s⁻¹ and a_foam = 0.33 (the latter value is based on the uncertainty in the foam albedo value proposed by Koepke 1984), and for v = 13 m s⁻¹ and no whitecaps. The curves clearly show that the effect of the whitecap reflection is very small when v = 7 m s⁻¹, as only 0.0003% of the surface is covered by whitecaps; the Lambertian reflection by the foam only slightly (≲ 1%) decreases the maximum P as compared to the case without whitecaps. For v = 13 m s⁻¹, the whitecaps cover 2.46% of the surface, and decrease the maximum P by about 5 and 10% for a_foam = 0.22 and 0.33, respectively. The whitecaps shift the maximum P from α ~ 100° to about 103° and 108°, respectively. We conclude that the smaller P at v = 13 m s⁻¹ for λ = 865 nm in Fig. 4, compared to P at v = 7 m s⁻¹, is indeed caused by the whitecaps.

If white caps are ignored (q = 0 in Eq. (A.14)), the degree of polarization P of light that is reflected by the ocean is independent of the wind speed v. The reason is that both the shadow function S and the probability density function for the wave facet inclination angles PDF, that depend on v (see Eqs. (A.1)–(A.6)), influence all elements of the ocean’s reflection matrix R_CO (and, if q = 0, matrix R_O) equally. The degree of polarization of an ocean planet as a whole can, however, depend (slightly) on v, because when v is large enough such that white caps cannot be ignored, light that is reflected by the foam and that is subsequently scattered within the atmosphere, will influence the total and polarized fluxes, and hence P of the light that emerges from the top of the atmosphere.

Appendix C Horizontally inhomogeneous planets and the weighted sum approach

A fast way to estimate the average flux and polarization phase curves of a partly cloudy planet, for any cloud coverage fraction f_c, is to compute the weighted sum of horizontally homogeneous planets (Stam 2008). That is, the total and polarized fluxes of a partly cloudy planet are derived from the respective total and polarized fluxes of a completely cloudy planet, F_cloudy, and a cloud-free planet, F_uncloudy, according to $$F(\lambda ,\alpha )=f_{\text{c}}{F}_{\text{cloudy}}(\lambda ,\alpha )+(1-f_{\text{c}})F_{\text{uncloudy}}(\lambda ,\alpha ).$$(C.1)

Because F_cloudy and F_uncloudy have to be computed only once, Eq. (C.1) can be evaluated very fast for arbitrary cloud coverage fractions f_c.

When using a weighted sum, the locations of cloud patches on a planet are not taken into account, while the locations can influence the total and polarized fluxes, and hence P. In this appendix, we compare results of Eq. (C.1) to the average phase curves of planets with various distributions of horizontally inhomogeneous patchy clouds (for which the locations of the patches are thus fully accounted for) (see Sect. 3.2.2).

Figure C.1 shows that the average phase curves of F, Q, and P for partly cloudy ocean planets with randomly distributed cloud patterns follow the phase curves computed using the weighted sum very accurately. For the curves in this figure, the wind speed v is 7 m s⁻¹ and the cloud coverage fraction f_c 0.50. However, we find that for any combination of v (1, 7, or 13 m s⁻¹) and f_c (0.25, 0.50, or 0.75) the agreement between the average phase curves and the weighted sums is excellent.

Fig. C.1 Comparison between F (left), Q (middle), and P (right) of the reflected light as computed using the weighted sum method (dashed lines) and as computed for pixelated, horizontally inhomogeneous ocean planets with patchy clouds (solid lines). The solid lines are the same as those shown in the second rows and second columns of Figs. 8–10. The cloud coverage fraction fc, which isthe weighting factor in the weighted sum approach, is 0.50 and the wind speed v is 7 m s−1. The wavelengths are the same as those used in Figs. 8–10, i.e., 350 nm (pink), 443 nm (blue), 550 nm (green), 670 nm (red), 865 nm (brown).

It should be emphasized that we use random patchy cloud cover patterns. Hence, in the case of a planet with a nonrandom cloud pattern, for example, a cloud pattern that is correlated to underlying continents, horizontally inhomogeneous planet models as presented in this paper may be needed to accurately model measured phase curves. In addition, using a weighted sum of horizontally homogeneous planets to compute phase curves of horizontally inhomogeneous planets does not provide information about the variability in the phase curves that is due to the variability in the cloud pattern, as shown by the shaded areas around the phase curves in Figs. 8–10.

References

All Tables

All Figures

Fig. 1 Disk-integrated reflected total flux F (left column), linearly polarized flux Q (middle column) and degree of polarization P (right column) of an ocean planet with a rough, Fresnel reflecting surface (top row) and a planet with a black surface (bottom row) as functions of the phase angle, for various wavelengths. Both planets have a pure gaseous atmosphere and are cloud-free. The wind speed above the ocean is 7 m s−1, and the surface coverage of white caps 0.28%. The black dashed line in the ocean planet plots (top row) pertains to a planet covered by the rough Fresnel surface, with a black surface below, thus without sub-surface ocean, and without atmosphere above.

In the text

Fig. 2 RGB colors of the (disk-integrated) reflected F (left column) and |Q| (right column) of the cloud-free ocean planet (top row) and the cloud-free black surface planet (bottom row) for surface pressures ranging from 0.5 to 10 bar in steps of 0.5 bar. Like in Fig. 1, the wind speed is 7 m s−1. The RGB colors are computed using weighted additive color mixing of the fluxes at λ = 443 nm (blue), 550 nm (green), and 670 nm (red), such that when these fluxes are equal, the color is white.

In the text

	Fig. 3 Similar to Fig. 1, but for the ocean planet, and for surface pressures equal to 0.5 bar (top row), 5 bar (middle row), and 10 bar (bottom row).
In the text

	Fig. 4 Similar to Fig. 1, but for the ocean planet and different wind speeds. The wind speeds and associated white cap percentages are: 1 m s−1 and 0.0003%(solid lines), 7 m s−1 and 0.28% (dashed lines) (cf. Fig. 1), 13 m s−1 and 2.46% (dotted lines).
In the text

	Fig. 5 Similar to Fig. 1, but for the ocean planet (solid lines) and a planet with the same gaseous atmosphere and the same rough Fresnel interface, but without the water body underneath (dashed lines). The wind speed v is 7 m s−1.
In the text

	Fig. 6 Similar to Fig. 1, but for the ocean planet with a horizontally homogeneous cloud deck. The wind speed is 1 m s−1 (solid line), 7 m s−1 (dashed line), or 13 m s−1 (dashed-dotted line). The dotted lines are the phase curves of a planet with a black surface and the same cloud deck.
In the text

Fig. 7 Disk-resolved (160 × 160 pixels) model output RGB color images in F (Cols. 1–3) and |Q| (Cols. 4–6)of ocean planets with patchy clouds. The brightness of a pixel is scaled with its average spectral F. Phase angle α = 80°. The cloud coverage fraction fc is 0.25 (Cols. 1 and 4), 0.5 (Cols. 2 and 5), or 0.75 (Cols. 3 and 6), and the wind speed v is 1 m s−1 (top row), 7 m s−1 (middle row), or 13 m s−1 (bottom row). For every combination of fc and v, the cloud pattern across the disk is different.

In the text

Fig. 8 Disk-integrated reflected F of ocean planets with patchy clouds for 5 wavelengths: 350 nm (pink), 443 nm (blue), 550 nm (green), 670 nm (red), and 865 nm (brown), and for wind speeds of 1 m s−1 (first row), 7 m s−1 (second row), and 13 m s−1 (third row), and cloud fractions fc of 0.25 (left column), 0.50 (middle column), and 0.75 (right column). Bottom row: corresponding curves forplanets with a black surface instead of an ocean. Each curve is the average of 300 curves computed for 300 different cloud patches configurations, which are different for each phase angle α. The 1-σ standard deviation is represented by the shaded areas around each curve (except for the longest wavelengths and largest phase angles, this deviation is very small).

In the text

	Fig. 9 Similar to Fig. 8, but for the linearly polarized flux Q.
In the text

	Fig. 10 Similar to Fig. 8, but for the degree of linear polarization P.
In the text

	Fig. 11 Similar to Fig. 7 for ocean planets (top) and black surface planets (bottom) covered by patchy clouds. The cloud coverage fraction fc is 0.50, and the wind speed v on the ocean planet is 7 m s−1. Phase anglesα range from 0° (left), through 30°, 60°, 80°, 100°, 120°, and 140°, to 160° (right).
In the text

Fig. 12 Similar to Fig. 2, but for cloud fractions ranging from 0.0 to 1.0 in steps of 0.01. The surface pressure is 1 bar. The results are shown for the phase curves generated by the weighted sum of cloud-free (fc = 0.0) and fully cloudy (fc = 1.0) planets (see Appendix C). The purple, yellow and cyan colors at α > 140° result from the high RGB color sensitivity on the very small absolute polarized fluxes |Q| for the different wavelengths at these phase angles (cf. Fig. 9).

In the text

Fig. 13 Phase angle αcc where the phase curves of the total flux F (black), the (linearly) polarized flux Q (gray), and the degree of (linear) polarization P (brown) of an ocean planet for the different wavelengths cross each other, as functions of the cloud coverage fraction fc. The wind speed v is 7 m s−1, although αcc does not significantly depend on v. The symbols(squares, circles, diamonds) are the αcc’s shown before in Figs. 8–10 as computed for horizontally inhomogeneous planets, with the error bars representing the 1σ standard deviation of 300 iterations with different locations of the patchy clouds. The solid, dashed, and dotted lines have been computed using the weighted sum approach (see text), with the clouds at different altitudes (pct = 0.5 or 0.7 bar) and for different cloud optical thicknesses (ba = 5 or 10). The horizontal axis at the top indicates the minimum orbital inclination angle i* required to be able to observe the exoplanet at the corresponding αcc (an inclination angle of 90° indicates a planetary orbit that is “seen” edge-on).

In the text

	Fig. B.1 Similar to Fig. 4, but for 865 nm and different wind speeds and whitecap properties. The speeds and the respective whitecap properties are 7 m s−1 and afoam = 0.22 (solid lines) (cf. Fig. 4), 13 m s−1 and afoam = 0.22 (dotted lines), 13 m s−1 and afoam = 0.33 (dashed lines), and 13 m s−1 but no whitecaps, i.e., q = 0 (dashed-dotted lines).
In the text

Fig. C.1 Comparison between F (left), Q (middle), and P (right) of the reflected light as computed using the weighted sum method (dashed lines) and as computed for pixelated, horizontally inhomogeneous ocean planets with patchy clouds (solid lines). The solid lines are the same as those shown in the second rows and second columns of Figs. 8–10. The cloud coverage fraction fc, which isthe weighting factor in the weighted sum approach, is 0.50 and the wind speed v is 7 m s−1. The wavelengths are the same as those used in Figs. 8–10, i.e., 350 nm (pink), 443 nm (blue), 550 nm (green), 670 nm (red), 865 nm (brown).

In the text