Skip to content

API Reference

Two-stream shortwave radiative transfer model for sea ice containing oil droplets.

constants

Constants used in the program

cts_wavelength

irradiances

Classes to store the solution of the continuous wavelength two stream model. Spectral irradiances and the integrated irradiances noralised by the incident.

CtsWavelengthIrradiance dataclass

One dimensional Arrays containing the upwelling and downwelling irradiances at each depth integrated over wavelength.

Irradiances are non-dimensional and need to be multiplied by the incident spectral radiation.

Parameters:

Name Type Description Default
z NDArray

vertical grid specified in dimensional units (m)

 required
upwelling NDArray

1D array of integrated upwelling irradiances

 required
downwelling NDArray

1D array of integrated downwelling irradiances

 required
Source code in oilrad/cts_wavelength/irradiances.py 
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
@dataclass(frozen=True)
class CtsWavelengthIrradiance:
    """One dimensional Arrays containing the upwelling and downwelling irradiances at each
    depth integrated over wavelength.

    Irradiances are non-dimensional and need to be multiplied by the incident spectral radiation.

    Args:
        z (NDArray): vertical grid specified in dimensional units (m)
        upwelling (NDArray): 1D array of integrated upwelling irradiances
        downwelling (NDArray): 1D array of integrated downwelling irradiances
    """

    z: NDArray
    upwelling: NDArray
    downwelling: NDArray

    _ice_base_index: int = 0

    @property
    def net_irradiance(self) -> NDArray:
        """Calculate net irradiance"""
        return self.downwelling - self.upwelling

    @property
    def albedo(self) -> NDArray:
        """Calculate albedo"""
        return self.upwelling[-1]

    @property
    def transmittance(self) -> NDArray:
        """Calculate transmittance at the ice ocean interface or the bottom
        of the domain if the domain is entirely ice."""
        return self.downwelling[self._ice_base_index]
albedo property

Calculate albedo

net_irradiance property

Calculate net irradiance

transmittance property

Calculate transmittance at the ice ocean interface or the bottom of the domain if the domain is entirely ice.

CtsWavelengthSpectralIrradiance dataclass

Two dimensional arrays containing the upwelling and downwelling irradiances at each depth and wavelength.

Irradiances are non-dimensional and need to be multiplied by the incident spectral radiation.

Parameters:

Name Type Description Default
z NDArray

vertical grid specified in dimensional units (m)

 required
wavelengths NDArray

array of wavelengths in nm

 required
upwelling NDArray

2D array of upwelling irradiances

 required
downwelling NDArray

2D array of downwelling irradiances

 required
Source code in oilrad/cts_wavelength/irradiances.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
@dataclass(frozen=True)
class CtsWavelengthSpectralIrradiance:
    """Two dimensional arrays containing the upwelling and downwelling irradiances at
    each depth and wavelength.

    Irradiances are non-dimensional and need to be multiplied by the incident spectral
    radiation.

    Args:
        z (NDArray): vertical grid specified in dimensional units (m)
        wavelengths (NDArray): array of wavelengths in nm
        upwelling (NDArray): 2D array of upwelling irradiances
        downwelling (NDArray): 2D array of downwelling irradiances
    """

    z: NDArray
    wavelengths: NDArray
    upwelling: NDArray
    downwelling: NDArray

    _ice_base_index: int = 0

    @property
    def net_irradiance(self) -> NDArray:
        return self.downwelling - self.upwelling

    @property
    def albedo(self) -> NDArray:
        """Calculate spectral albedo"""
        return self.upwelling[-1, :]

    @property
    def transmittance(self) -> NDArray:
        """Calculate spectral transmittance at the ice ocean interface or the bottom
        of the domain if the domain is entirely ice."""
        return self.downwelling[self._ice_base_index, :]
albedo property

Calculate spectral albedo

transmittance property

Calculate spectral transmittance at the ice ocean interface or the bottom of the domain if the domain is entirely ice.

model

Solve the two-stream shortwave radiation model for a layer of ice with constant optical properties aside from absorption varying due to the oil mass ratio.

Arbitrary wavelength resolution

CtsWavelengthModel dataclass

Class containing all the necessary parameters to solve the two-stream shortwave radiative transfer model in a domain with continuously varying liquid fraction and oil mass ratio.

If no array is provided for liquid fraction, it is assumed to be zero everywhere. This corresponds to a completely frozen domain.

Oil mass ratio is provided in ng oil / g ice and, along with the median droplet radius for the oil droplet distribution, is used to calculate the absorption coefficient by interpolating data for Romashkino oil from Redmond Roche et al. 2022.

Parameters:

Name Type Description Default
z NDArray

vertical grid in meters

 required
wavelengths NDArray

array of wavelengths in nm

 required
oil_mass_ratio NDArray

array of oil mass ratio in ng oil / g ice on the vertical grid

 required
ice_scattering_coefficient float

scattering coefficient for ice in 1/m

 required
median_droplet_radius_in_microns float

median droplet radius in microns

 required
absorption_enhancement_factor float

enhancement factor for oil absorption appropriate for the two-stream model

 1
liquid_fraction NDArray

liquid fraction array on the vertical grid

 None
fast_solve Bool

if True, solve the model only for wavelengths below a wavelength cutoff, assume longer wavelengths are absorbed at the surface

 False
wavelength_cutoff float

cutoff wavelength in nm

 None
Source code in oilrad/cts_wavelength/model.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
@dataclass
class CtsWavelengthModel:
    """Class containing all the necessary parameters to solve the two-stream shortwave
    radiative transfer model in a domain with continuously varying liquid fraction and
    oil mass ratio.

    If no array is provided for liquid fraction, it is assumed to be zero everywhere.
    This corresponds to a completely frozen domain.

    Oil mass ratio is provided in ng oil / g ice and, along with the median droplet radius
    for the oil droplet distribution, is used to calculate the absorption coefficient
    by interpolating data for Romashkino oil from Redmond Roche et al. 2022.

    Args:
        z (NDArray): vertical grid in meters
        wavelengths (NDArray): array of wavelengths in nm
        oil_mass_ratio (NDArray): array of oil mass ratio in ng oil / g ice on the vertical grid
        ice_scattering_coefficient (float): scattering coefficient for ice in 1/m
        median_droplet_radius_in_microns (float): median droplet radius in microns
        absorption_enhancement_factor (float): enhancement factor for oil absorption appropriate for the two-stream model
        liquid_fraction (NDArray): liquid fraction array on the vertical grid
        fast_solve (Bool): if True, solve the model only for wavelengths below a wavelength cutoff, assume longer wavelengths are absorbed at the surface
        wavelength_cutoff (float): cutoff wavelength in nm
    """

    z: NDArray
    wavelengths: NDArray
    oil_mass_ratio: NDArray
    ice_scattering_coefficient: float
    median_droplet_radius_in_microns: float
    absorption_enhancement_factor: float = 1
    liquid_fraction: Optional[NDArray] = None
    fast_solve: bool = False
    wavelength_cutoff: Optional[float] = None

    def __post_init__(self):
        # initialise liquid fraction as zero everywhere if not provided
        if self.liquid_fraction is None:
            self.liquid_fraction = np.full_like(self.z, 0)

        # find the index of the ice ocean interface
        self._ice_base_index = np.argmax(self.liquid_fraction < 1)

solve_at_given_wavelength(model, wavelength)

Use the scipy solve_bvp function to solve the two-stream model as a function of depth for a given wavelenght value.

Parameters:

Name Type Description Default
model CtsWavelengthModel

model parameters

 required
wavelength float

wavelength in nm

 required

Returns: tuple[NDArray, NDArray]: upwelling and downwelling irradiances as functions of depth Raises: RuntimeError: if the solver does not converge

Source code in oilrad/cts_wavelength/model.py 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
def solve_at_given_wavelength(model, wavelength: float) -> tuple[NDArray, NDArray]:
    """Use the scipy solve_bvp function to solve the two-stream model as a function of
    depth for a given wavelenght value.

    Args:
        model (CtsWavelengthModel): model parameters
        wavelength (float): wavelength in nm
    Returns:
        tuple[NDArray, NDArray]: upwelling and downwelling irradiances as functions of depth
    Raises:
        RuntimeError: if the solver does not converge
    """
    fun = _get_ODE_fun(model, wavelength)
    solution = solve_bvp(
        fun,
        _BCs,
        np.linspace(model.z[0], model.z[-1], 5),
        np.zeros((2, 5)),
        max_nodes=12000,
    )
    if not solution.success:
        raise RuntimeError(f"{solution.message}")
    return solution.sol(model.z)[0], solution.sol(model.z)[1]

integrate

Spectrally integrate two-stream model solutions

integrate_over_SW(spectral_irradiance, spectrum=None)

Spectrally integrate the two-stream model solution

When integrating the continuous wavelength model need to supply an instance of the BlackBodySpectrum

Integration of the six band model uses the cloudy incident spectrum of Grenfell et al 2004.

Parameters:

Name Type Description Default
spectral_irradiance SpectralIrradiance | SixBandSpectralIrradiance

spectral two-stream model solution

 required
spectrum Optional[BlackBodySpectrum]

normalised incident shortwave spectrum (only needed when integrating SpectralIrradiance). Currently only the BlackBodySpectrum is implemented.

 None

Returns: Irradiance | SixBandIrradiance: spectrally integrated irradiances

Source code in oilrad/integrate.py 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def integrate_over_SW(
    spectral_irradiance: CtsWavelengthSpectralIrradiance | SixBandSpectralIrradiance,
    spectrum: Optional[BlackBodySpectrum] = None,
) -> CtsWavelengthIrradiance | SixBandIrradiance:
    """Spectrally integrate the two-stream model solution

    When integrating the continuous wavelength model need to supply an instance of the
    BlackBodySpectrum

    Integration of the six band model uses the cloudy incident spectrum of Grenfell
    et al 2004.

    Args:
        spectral_irradiance (SpectralIrradiance | SixBandSpectralIrradiance):
            spectral two-stream model solution
        spectrum (Optional[BlackBodySpectrum]): normalised incident shortwave spectrum
            (only needed when integrating SpectralIrradiance).
            Currently only the BlackBodySpectrum is implemented.
    Returns:
        Irradiance | SixBandIrradiance: spectrally integrated irradiances
    """
    if isinstance(spectral_irradiance, CtsWavelengthSpectralIrradiance):
        if spectrum is None:
            raise ValueError("spectrum must be provided")
        wavelengths = spectral_irradiance.wavelengths
        integrate = lambda irradiance: trapezoid(
            irradiance * spectrum(wavelengths), wavelengths, axis=1
        )
        integrated_upwelling = integrate(spectral_irradiance.upwelling)
        integrated_downwelling = integrate(spectral_irradiance.downwelling)
        return CtsWavelengthIrradiance(
            spectral_irradiance.z,
            integrated_upwelling,
            integrated_downwelling,
            _ice_base_index=spectral_irradiance._ice_base_index,
        )
    if isinstance(spectral_irradiance, SixBandSpectralIrradiance):
        integrate = lambda x: np.sum(x * CLOUDY_SKY_FRACTIONS, axis=1)
        return SixBandIrradiance(
            spectral_irradiance.z,
            integrate(spectral_irradiance.upwelling),
            integrate(spectral_irradiance.downwelling),
            np.sum(CLOUDY_SKY_FRACTIONS * spectral_irradiance.albedo),
            _ice_base_index=spectral_irradiance._ice_base_index,
        )
    raise NotImplementedError()

optics

Module to calculate the optical properties for ice containing oil droplets: - The spectral absorption coefficient which depends on oil mass concentration and droplet size - The wavelength-independent scattering coefficient which takes a constant value in the ice and is zero the liquid

Load data for imaginary refractive index against wavelength from doi:10.1029/2007JD009744. This is used to calculate the absorption coefficient of pure ice. To interpolate the data to other wavelengths should interpolate the log of the data linearly.

Oil absorption calculated following Redmond Roche et al 2022 using given data for mass absorption coefficient of oil in ice which is a function of wavelength and droplet radius

calculate_ice_absorption_coefficient(wavelength_in_nm)

calculate ice absorption coefficient from Warren 2008 data at given wavelengths inputted in nano meters from interpolated imaginary refractive index data

Source code in oilrad/optics.py 
70
71
72
73
74
75
76
77
78
79
def calculate_ice_absorption_coefficient(wavelength_in_nm):
    """calculate ice absorption coefficient from Warren 2008 data at given
    wavelengths inputted in nano meters from interpolated imaginary refractive index
    data"""
    wavelengths_in_m = wavelength_in_nm * 1e-9
    imaginary_refractive_index = _calculate_ice_imaginary_refractive_index(
        wavelength_in_nm * 1e-3
    )
    absorption_coefficient = 4 * np.pi * imaginary_refractive_index / wavelengths_in_m
    return absorption_coefficient

calculate_ice_oil_absorption_coefficient(wavelengths_in_nm, oil_mass_ratio, droplet_radius_in_microns, absorption_enhancement_factor=1.0)

Calculate the absorption coefficient in 1/m of ice polluted with oil droplets following roche et al 2022. The oil droplets radii are distributed log-normally with geometric standard deviation e. We specify the median radius for the distribution.

mass ratio in units of ng oil / g ice

This is for Romashkino oil.

The enahncement factor is an ad hoc correction for the two stream model to try and better match the results of redmondroche2022 which used an 8-stream model

Parameters:

Name Type Description Default
wavelengths_in_nm NDArray

wavelengths in nm

 required
oil_mass_ratio float

mass ratio of oil in ice in ng/g

 required
droplet_radius_in_microns float

median droplet radius in microns

 required
absorption_enhancement_factor float

enhancement factor for absorption coefficient. Defaults to 1.0.

 1.0

Returns:

Name Type Description
NDArray

absorption coefficient in 1/m
Source code in oilrad/optics.py 
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
def calculate_ice_oil_absorption_coefficient(
    wavelengths_in_nm,
    oil_mass_ratio,
    droplet_radius_in_microns,
    absorption_enhancement_factor=1.0,
):
    """Calculate the absorption coefficient in 1/m of ice polluted with oil droplets
    following roche et al 2022. The oil droplets radii are distributed log-normally
    with geometric standard deviation e. We specify the median radius for the distribution.

    mass ratio in units of ng oil / g ice

    This is for Romashkino oil.

    The enahncement factor is an ad hoc correction for the two stream model to try and
    better match the results of redmondroche2022 which used an 8-stream model

    Args:
        wavelengths_in_nm (NDArray): wavelengths in nm
        oil_mass_ratio (float): mass ratio of oil in ice in ng/g
        droplet_radius_in_microns (float): median droplet radius in microns
        absorption_enhancement_factor (float, optional): enhancement factor for absorption coefficient. Defaults to 1.0.

    Returns:
        NDArray: absorption coefficient in 1/m
    """
    mass_ratio_dimensionless = oil_mass_ratio * 1e-9
    return absorption_enhancement_factor * (
        calculate_ice_absorption_coefficient(wavelengths_in_nm)
        + mass_ratio_dimensionless
        * 1e3
        * ICE_DENSITY_ROCHE_2022
        * Romashkino_MAC(wavelengths_in_nm, droplet_radius_in_microns)
    )

calculate_scattering(liquid_fraction, ice_scattering_coefficient)

Calculate scattering coefficient in ice and return zero in liquid. (doesn't depend on wavelength)

Smoothly transitions from the ice_scattering_coefficient to zero in the liquid using a tanh function.

Parameters:

Name Type Description Default
liquid_fraction NDArray

liquid fraction as a function of depth

 required
ice_scattering_coefficient float

scattering coefficient of ice

 required

Returns:

Name Type Description
NDArray

scattering coefficient [1/m] as a function of depth
Source code in oilrad/optics.py 
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
def calculate_scattering(liquid_fraction: NDArray, ice_scattering_coefficient: float):
    """Calculate scattering coefficient in ice and return zero in liquid.
    (doesn't depend on wavelength)

    Smoothly transitions from the ice_scattering_coefficient to zero in the liquid using a tanh function.

    Args:
        liquid_fraction (NDArray): liquid fraction as a function of depth
        ice_scattering_coefficient (float): scattering coefficient of ice

    Returns:
        NDArray: scattering coefficient [1/m] as a function of depth
    """

    return ice_scattering_coefficient * np.tanh((1 - liquid_fraction) * 100)

six_band

irradiances

Classes to store spectral and integrated solution of the six band model. All irradiances are normalised by the incident.

SixBandIrradiance dataclass

One dimensional Arrays containing the upwelling and downwelling irradiances at each depth integrated over the wavelength bands.

Irradiances are non-dimensional and need to be multiplied by the incident spectral radiation.

Parameters:

Name Type Description Default
z NDArray

vertical grid specified in dimensional units (m)

 required
upwelling NDArray

1D array of integrated upwelling irradiances

 required
downwelling NDArray

1D array of integrated downwelling irradiances

 required
albedo float

spectrally integrated albedo (including the snow layer and SSL)

 required
Source code in oilrad/six_band/irradiances.py 
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
@dataclass(frozen=True)
class SixBandIrradiance:
    """One dimensional Arrays containing the upwelling and downwelling irradiances at each
    depth integrated over the wavelength bands.

    Irradiances are non-dimensional and need to be multiplied by the incident spectral radiation.

    Args:
        z (NDArray): vertical grid specified in dimensional units (m)
        upwelling (NDArray): 1D array of integrated upwelling irradiances
        downwelling (NDArray): 1D array of integrated downwelling irradiances
        albedo (float): spectrally integrated albedo (including the snow layer and SSL)
    """

    z: NDArray
    upwelling: NDArray
    downwelling: NDArray
    albedo: float

    _ice_base_index: int = 0

    @property
    def net_irradiance(self) -> NDArray:
        return self.downwelling - self.upwelling

    @property
    def transmittance(self) -> NDArray:
        """Calculate transmittance at the ice ocean interface or the bottom
        of the domain if the domain is entirely ice."""
        return self.downwelling[self._ice_base_index]
transmittance property

Calculate transmittance at the ice ocean interface or the bottom of the domain if the domain is entirely ice.

SixBandSpectralIrradiance dataclass

Two dimensional arrays containing the upwelling and downwelling irradiances at each depth and in each wavelength band.

Irradiances are non-dimensional and need to be multiplied by the incident spectral radiation.

Parameters:

Name Type Description Default
z NDArray

vertical grid specified in dimensional units (m)

 required
upwelling NDArray

2D array of upwelling irradiances

 required
downwelling NDArray

2D array of downwelling irradiances

 required
albedo NDArray

1D array of spectral albedo (including snow and SSL)

 required
Source code in oilrad/six_band/irradiances.py 
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
@dataclass(frozen=True)
class SixBandSpectralIrradiance:
    """Two dimensional arrays containing the upwelling and downwelling irradiances at
    each depth and in each wavelength band.

    Irradiances are non-dimensional and need to be multiplied by the incident spectral
    radiation.

    Args:
        z (NDArray): vertical grid specified in dimensional units (m)
        upwelling (NDArray): 2D array of upwelling irradiances
        downwelling (NDArray): 2D array of downwelling irradiances
        albedo (NDArray): 1D array of spectral albedo (including snow and SSL)
    """

    z: NDArray
    upwelling: NDArray
    downwelling: NDArray
    albedo: NDArray

    _ice_base_index: int = 0

    @property
    def net_irradiance(self) -> NDArray:
        return self.downwelling - self.upwelling

    @property
    def transmittance(self) -> NDArray:
        """Calculate spectral transmittance at the ice ocean interface or the bottom
        of the domain if the domain is entirely ice."""
        return self.downwelling[self._ice_base_index, :]

    @property
    def PAR_transmittance(self) -> NDArray:
        """Calculate plane PAR transmittance as the ratio of the
        net irradiant power in the PAR range (400-700nm) to the incident
        irradiative power at the ice / snow surface.
        """
        return np.sum(
            self.net_irradiance[:, 1:4] * CLOUDY_SKY_FRACTIONS[1:4], axis=1
        ) / np.sum(CLOUDY_SKY_FRACTIONS[1:4])

    @property
    def plane_PAR(self) -> NDArray:
        """Calculate plane PAR normalised by the incident broadband shortwave irradiance.

        To convert to micromol-photns m^-2 s^-1 we need to multiply by the incident shortwave
        irradiance in W m^-2.

        To convert to the scalar value for an isotropic downwelling irradiance multiply
        by a factor of 2.
        """
        PAR_weightings = np.array(
            [9.809215701925024e-08, 1.0750608149135324e-07, 1.0006054876321231e-07]
        )
        return (1e6 / (PLANCK * LIGHTSPEED * MOLE)) * np.sum(
            (self.upwelling[:, 1:4] + self.downwelling[:, 1:4]) * PAR_weightings, axis=1
        )

    @property
    def ice_base_PAR_transmittance(self) -> float:
        return self.PAR_transmittance[self._ice_base_index]

    @property
    def ice_base_plane_PAR(self) -> float:
        return self.plane_PAR[self._ice_base_index]
PAR_transmittance property

Calculate plane PAR transmittance as the ratio of the net irradiant power in the PAR range (400-700nm) to the incident irradiative power at the ice / snow surface.

plane_PAR property

Calculate plane PAR normalised by the incident broadband shortwave irradiance.

To convert to micromol-photns m^-2 s^-1 we need to multiply by the incident shortwave irradiance in W m^-2.

To convert to the scalar value for an isotropic downwelling irradiance multiply by a factor of 2.

transmittance property

Calculate spectral transmittance at the ice ocean interface or the bottom of the domain if the domain is entirely ice.

model

Shortwave radiative transfer model for a layer of sea ice with 6 spectral bands.

Optionally the ice may have a melt pond layer, a snow layer and a surface scattering layer (SSL).

SixBandModel dataclass

Class containing all the necessary parameters to solve the two-stream shortwave radiative transfer model in a domain with continuously varying liquid fraction and oil mass ratio with optical properties averaged in the six spectral bands with wavelengths:

300-400nm 400-500nm 500-600nm 600-700nm 700-1200nm 1200-3000nm

Irradiances are scaled by the incident downwelling in each spectral band.

If no array is provided for liquid fraction, it is assumed to be zero everywhere. This corresponds to a completely frozen domain.

Oil mass ratio is provided in ng oil / g ice and, along with the median droplet radius for the oil droplet distribution, is used to calculate the absorption coefficient by interpolating data for Romashkino oil from Redmond Roche et al. 2022.

Parameters:

Name Type Description Default
z NDArray

vertical grid in meters

 required
oil_mass_ratio NDArray

array of oil mass ratio in ng oil / g ice on the vertical grid

 required
ice_scattering_coefficient float

scattering coefficient for ice in 1/m

 required
median_droplet_radius_in_microns float

median droplet radius in microns

 required
absorption_enhancement_factor float

enhancement factor for oil absorption appropriate for the two-stream model

 required
snow_depth float

snow depth in meters

 required
snow_spectral_albedos NDArray

spectral albedos for the snow layer in each band

 required
snow_extinction_coefficients NDArray

spectral extinction coefficient for the snow layer in each band

 required
SSL_depth float

depth of the surface scattering layer in meters

 required
SSL_spectral_albedos NDArray

spectral albedos for the SSL in each band

 required
SSL_extinction_coefficients NDArray

spectral extinction coefficients for the SSL in each band

 required
liquid_fraction Optional[NDArray]

liquid fraction array on the vertical grid

 None
Source code in oilrad/six_band/model.py 
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
@dataclass
class SixBandModel:
    """Class containing all the necessary parameters to solve the two-stream shortwave
    radiative transfer model in a domain with continuously varying liquid fraction and
    oil mass ratio with optical properties averaged in the six spectral bands with
    wavelengths:

    300-400nm
    400-500nm
    500-600nm
    600-700nm
    700-1200nm
    1200-3000nm

    Irradiances are scaled by the incident downwelling in each spectral band.

    If no array is provided for liquid fraction, it is assumed to be zero everywhere.
    This corresponds to a completely frozen domain.

    Oil mass ratio is provided in ng oil / g ice and, along with the median droplet radius
    for the oil droplet distribution, is used to calculate the absorption coefficient
    by interpolating data for Romashkino oil from Redmond Roche et al. 2022.

    Args:
        z (NDArray): vertical grid in meters
        oil_mass_ratio (NDArray): array of oil mass ratio in ng oil / g ice on the
            vertical grid
        ice_scattering_coefficient (float): scattering coefficient for ice in 1/m
        median_droplet_radius_in_microns (float): median droplet radius in microns
        absorption_enhancement_factor (float): enhancement factor for oil absorption
            appropriate for the two-stream model
        snow_depth (float): snow depth in meters
        snow_spectral_albedos (NDArray): spectral albedos for the snow layer in each
            band
        snow_extinction_coefficients (NDArray): spectral extinction coefficient for the
            snow layer in each band
        SSL_depth (float): depth of the surface scattering layer in meters
        SSL_spectral_albedos (NDArray): spectral albedos for the SSL in each band
        SSL_extinction_coefficients (NDArray): spectral extinction coefficients for the
            SSL in each band
        liquid_fraction (Optional[NDArray]): liquid fraction array on the vertical grid
    """

    z: NDArray
    oil_mass_ratio: NDArray
    ice_scattering_coefficient: float
    median_droplet_radius_in_microns: float
    absorption_enhancement_factor: float

    snow_depth: float
    snow_spectral_albedos: NDArray
    snow_extinction_coefficients: NDArray

    SSL_depth: float
    SSL_spectral_albedos: NDArray
    SSL_extinction_coefficients: NDArray

    liquid_fraction: Optional[NDArray] = None

    bands: ClassVar[List[Tuple[int, int]]] = WAVELENGTH_BANDS

    def __post_init__(self):
        # initialise liquid fraction as zero everywhere if not provided
        if self.liquid_fraction is None:
            self.liquid_fraction = np.full_like(self.z, 0)

        # find the index of the ice ocean interface
        self._ice_base_index = np.argmax(self.liquid_fraction < 1)

        # Generate band average ice absorption
        UV_adjusted_bands = [(350, 400)] + self.bands[1:]  # Just average UV in 350-400
        self.band_average_ice_absorption = np.array(
            [
                np.mean(
                    calculate_ice_absorption_coefficient(
                        np.linspace(band[0], band[1], 1000)
                    )
                )
                for band in UV_adjusted_bands
            ]
        )
        # Generate band average oil MAC
        self.band_average_Romashkino_MAC = np.array(
            [
                np.mean(
                    Romashkino_MAC(
                        np.linspace(band[0], band[1], 1000),
                        self.median_droplet_radius_in_microns,
                    )
                )
                for band in UV_adjusted_bands
            ]
        )

solve_a_wavelength_band(model, wavelength_band_index)

Use the scipy solve_bvp function to solve the two-stream model as a function of depth for each wavelength band.

Parameters:

Name Type Description Default
model SixBandModel

model parameters

 required
wavelength_band_index int

index of the wavelength band to solve

 required

Returns: tuple[NDArray, NDArray]: upwelling and downwelling irradiances as functions of depth Raises: RuntimeError: if the solver does not converge

Source code in oilrad/six_band/model.py 
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
def solve_a_wavelength_band(
    model: SixBandModel, wavelength_band_index: int
) -> tuple[NDArray, NDArray]:
    """Use the scipy solve_bvp function to solve the two-stream model as a function of
    depth for each wavelength band.

    Args:
        model (SixBandModel): model parameters
        wavelength_band_index (int): index of the wavelength band to solve
    Returns:
        tuple[NDArray, NDArray]: upwelling and downwelling irradiances as functions of depth
    Raises:
        RuntimeError: if the solver does not converge
    """
    # Add radiaition absorbed in SSL into the top of the ice
    if model.snow_depth == 0:
        absorbed_in_SSL = (
            1
            - calculate_band_SSL_albedo(model, wavelength_band_index)
            - calculate_band_surface_transmittance(model, wavelength_band_index)
        )
    else:
        absorbed_in_SSL = 0

    # In high wavelength band just assume all radiation is absorbed at ice surface
    if wavelength_band_index == 5:
        upwelling = np.zeros_like(model.z)
        downwelling = np.zeros_like(model.z)
        downwelling[-1] = (
            calculate_band_surface_transmittance(model, 5) + absorbed_in_SSL
        )
        return upwelling, downwelling

    fun = _get_ODE_fun(model, wavelength_band_index)
    BCs = _get_BC_fun(model, wavelength_band_index)
    solution = solve_bvp(
        fun,
        BCs,
        np.linspace(model.z[0], model.z[-1], 5),
        np.zeros((2, 5)),
        max_nodes=12000,
    )
    if not solution.success:
        raise RuntimeError(f"{solution.message}")

    upwelling = solution.sol(model.z)[0]
    downwelling = solution.sol(model.z)[1]
    downwelling[-1] = downwelling[-1] + absorbed_in_SSL
    return upwelling, downwelling

top_surface

Calculate surface transmittance for six band model through a snow layer and SSL

solve

Provide a function to solve the two-stream model in the case of continuously varying optical properties which implements a faster solve approximation for long wavelengths if the fast_solve parameter of the model is set to True for the continuous wavelength case.

solve_two_stream_model(model)

Solve the two-stream model and return an object containing the solution at all specified wavelengths

Args (InfiniteLayerModel | SixBandModel): model: two-stream model parameters

Returns:

Type Description
CtsWavelengthSpectralIrradiance | SixBandSpectralIrradiance

SpectralIrradiance | SixBandSpectralIrradiance: object containing the solution of the two-stream model at each wavelength
Source code in oilrad/solve.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
def solve_two_stream_model(
    model: CtsWavelengthModel | SixBandModel,
) -> CtsWavelengthSpectralIrradiance | SixBandSpectralIrradiance:
    """Solve the two-stream model and return an object containing the solution at all
    specified wavelengths

    Args (InfiniteLayerModel | SixBandModel):
        model: two-stream model parameters

    Returns:
        SpectralIrradiance | SixBandSpectralIrradiance: object containing the solution of the two-stream model at each wavelength
    """

    if isinstance(model, CtsWavelengthModel):
        upwelling = np.empty((model.z.size, model.wavelengths.size))
        downwelling = np.empty((model.z.size, model.wavelengths.size))
        if model.fast_solve:
            cut_off_index = (
                np.argmin(np.abs(model.wavelengths - model.wavelength_cutoff)) + 1
            )
            is_surface = np.s_[cut_off_index:]
            is_interior = np.s_[:cut_off_index]
            for i, wavelength in enumerate(model.wavelengths[is_interior]):
                col_upwelling, col_downwelling = solve_at_given_wavelength(
                    model, wavelength
                )
                upwelling[:, i] = col_upwelling
                downwelling[:, i] = col_downwelling

            upwelling[:, is_surface] = 0
            downwelling[:, is_surface] = 0
            downwelling[-1, is_surface] = 1
        else:
            for i, wavelength in enumerate(model.wavelengths):
                col_upwelling, col_downwelling = solve_at_given_wavelength(
                    model, wavelength
                )
                upwelling[:, i] = col_upwelling
                downwelling[:, i] = col_downwelling
        return CtsWavelengthSpectralIrradiance(
            model.z, model.wavelengths, upwelling, downwelling, model._ice_base_index
        )

    if isinstance(model, SixBandModel):
        upwelling = np.empty((model.z.size, 6))
        downwelling = np.empty((model.z.size, 6))
        for index in WAVELENGTH_BAND_INDICES:
            col_upwelling, col_downwelling = solve_a_wavelength_band(model, index)
            upwelling[:, index] = col_upwelling
            downwelling[:, index] = col_downwelling

        ice_albedo = upwelling[-1, :]
        surface_albedo = np.array(
            [calculate_band_surface_albedo(model, i) for i in WAVELENGTH_BAND_INDICES]
        )
        surface_transmittance = np.array(
            [
                calculate_band_surface_transmittance(model, i)
                for i in WAVELENGTH_BAND_INDICES
            ]
        )
        albedo = surface_albedo + surface_transmittance * ice_albedo
        return SixBandSpectralIrradiance(
            model.z,
            upwelling,
            downwelling,
            albedo,
            model._ice_base_index,
        )

    raise NotImplementedError("Model type not recognized")

spectra

Module to provide the spectrum of incident downwelling shortwave radiation at the top of the domain.

Currently only the black body spectrum is integrated which follows the solar spectrum at the top of the atmosphere, but is normalised to integrate to one in the shortwave range.

Data used to compute the Planck function is from: https://www.oceanopticsbook.info/view/light-and-radiometry/level-2/light-from-the-sun

BlackBodySpectrum dataclass

Spectrum with blackbody shape that integrates to 1 between minimum and maximum wavelength specified in nm

Once initialised the spectrum can be called with an array of wavelengths in nm.

Parameters:

Name Type Description Default
min_wavelength float

minimum wavelength in nm

 required
max_wavelength float

maximum wavelength in nm

 required

Raises:

Type Description
ValueError

if wavelength is not in the shortwave range
Source code in oilrad/spectra.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
@dataclass(frozen=True)
class BlackBodySpectrum:
    """Spectrum with blackbody shape that integrates to 1 between minimum and maximum
    wavelength specified in nm

    Once initialised the spectrum can be called with an array of wavelengths in nm.

    Args:
        min_wavelength: minimum wavelength in nm
        max_wavelength: maximum wavelength in nm

    Raises:
        ValueError: if wavelength is not in the shortwave range
    """

    min_wavelength: float
    max_wavelength: float

    @cached_property
    def _total_irradiance(self) -> float:
        return quad(
            self._top_of_atmosphere_irradiance, self.min_wavelength, self.max_wavelength
        )[0]

    @classmethod
    def _top_of_atmosphere_irradiance(cls, wavelength_in_nm):
        """For wavelength in nm and temperature in K return top of atmosphere solar
        irradiance in W/m2 nm
        https://www.oceanopticsbook.info/view/light-and-radiometry/level-2/blackbody-radiation
        """
        return (
            _solar_planck_function(wavelength_in_nm * 1e-9, T=5782)
            * (SUN_RADIUS**2 / AU**2)
            * np.pi
            * 1e-9
        )

    def __call__(self, wavelength_in_nm: NDArray) -> NDArray:
        if np.any(wavelength_in_nm > self.max_wavelength) or np.any(
            wavelength_in_nm < self.min_wavelength
        ):
            raise ValueError(
                f"wavelength not in shortwave range {self.min_wavelength}nm - {self.max_wavelength}nm"
            )
        return (
            self._top_of_atmosphere_irradiance(wavelength_in_nm)
            / self._total_irradiance
        )