Mie Scattering Databases#

import sasktran2 as sk
import numpy as np
import matplotlib.pyplot as plt

SASKTRAN2 includes the ability to generate scattering properties for Mie scatterers. Generally a Mie scattering atmospheric constituent is described by a particle size distribution, and a refractive index. Since Mie scattering calculations can be slow, SASKTRAN2 does not contain perform Mie calculations alongside the radiative transfer calculation, instead databases are generated using internal tools that can then be re-used.

Mie scatterers are parameterized by two quantities, the refractive index, and a particle size distribution. SASKTRAN2 contains several built in refractive index functions (see Refractive Index Containers), letโ€™s create one for sulfates, sasktran2.mie.refractive.H2SO4

refrac = sk.mie.refractive.H2SO4()

wavelengths = np.arange(280, 2000, 0.1)
plt.plot(wavelengths, np.real(refrac.refractive_index(wavelengths)))
plt.xlabel("Wavelength [nm]")
plt.ylabel("Real Part of Refractive Index")
Text(0, 0.5, 'Real Part of Refractive Index')
../_images/be844dea75dfc40a8900fe8e4182ba5d7cf919cfca7056d7c9e08b70af8b1bb9.png

it is also possible to define your own refractive index entirely using the sasktran2.mie.refractive.RefractiveIndex object.

The second quantity that we need is a particle size distribution (see Particle Size Distributions). Here we will create a log-normal particle size distribution

distribution = sk.mie.distribution.LogNormalDistribution()

Every particle size distribution is parameterized by several parameters, we can look at them explicity,

distribution.args()
['median_radius', 'mode_width']

And we can also inspect the distribution for specific values of the parameters (median_radius, mode_width)

radii = np.arange(1, 500, 0.1)
plt.plot(radii, distribution.distribution(median_radius=160, mode_width=1.6).pdf(radii))
plt.xlabel("Radius [nm]")
plt.ylabel("PDF")
Text(0, 0.5, 'PDF')
../_images/026c2654030787bfc8c202649f124ddec84589fd34051bda4d3eb49346f215bd.png

Now that we have a particle size distribution and a refractive index, we can create our Mie scattering database.

mie_db = sk.database.MieDatabase(
    distribution,
    refrac,
    wavelengths_nm=np.arange(270, 1000, 50.0),
    median_radius=[100, 150, 200],
    mode_width=[1.5, 1.7],
)

The database is a function of wavelength, and any arguments of the particle size distribution, median_radius and mode_width in our case. When it is created for the first time the local database will be generated. Any subsequent instantiations of the object will re-use the cached database.

We can look at the created database

mie_db.load_ds()
<xarray.Dataset> Size: 648kB
Dimensions:        (wavelength_nm: 15, median_radius: 3, mode_width: 2,
                    cos_angle: 128, legendre: 64)
Coordinates:
  * wavelength_nm  (wavelength_nm) float64 120B 270.0 320.0 ... 920.0 970.0
  * median_radius  (median_radius) int64 24B 100 150 200
  * mode_width     (mode_width) float64 16B 1.5 1.7
  * cos_angle      (cos_angle) float64 1kB -0.9993 -0.9964 ... 0.9999 1.0
Dimensions without coordinates: legendre
Data variables: (12/13)
    xs_scattering  (wavelength_nm, median_radius, mode_width) float64 720B ...
    xs_total       (wavelength_nm, median_radius, mode_width) float64 720B ...
    xs_absorption  (wavelength_nm, median_radius, mode_width) float64 720B ...
    p11            (wavelength_nm, cos_angle, median_radius, mode_width) float64 92kB ...
    p12            (wavelength_nm, cos_angle, median_radius, mode_width) float64 92kB ...
    p33            (wavelength_nm, cos_angle, median_radius, mode_width) float64 92kB ...
    ...             ...
    lm_a1          (wavelength_nm, legendre, median_radius, mode_width) float64 46kB ...
    lm_a2          (wavelength_nm, legendre, median_radius, mode_width) float64 46kB ...
    lm_a3          (wavelength_nm, legendre, median_radius, mode_width) float64 46kB ...
    lm_a4          (wavelength_nm, legendre, median_radius, mode_width) float64 46kB ...
    lm_b1          (wavelength_nm, legendre, median_radius, mode_width) float64 46kB ...
    lm_b2          (wavelength_nm, legendre, median_radius, mode_width) float64 46kB ...

Which contains the scattering and absorption cross sections, as well as the Legendre expansion moments of the phase function that are necessary for the radiative transfer calculation.

Including Mie Scatterers in the Radiative Transfer Calculation#

The Mie database can be included and used as an optical property with the standard sasktran2.constituent.NumberDensityScatterer and sasktran2.constituent.ExtinctionScatterer constituents.

First we set up the example calculation without aerosol

config = sk.Config()

altitudes_m = np.arange(0, 65001, 1000)

model_geometry = sk.Geometry1D(cos_sza=0.6,
                               solar_azimuth=0,
                               earth_radius_m=6372000,
                               altitude_grid_m=altitudes_m,
                               interpolation_method=sk.InterpolationMethod.LinearInterpolation,
                               geometry_type=sk.GeometryType.Spherical)

viewing_geo = sk.ViewingGeometry()

for alt in [10000, 20000, 30000, 40000]:
    ray = sk.TangentAltitudeSolar(tangent_altitude_m=alt,
                                    relative_azimuth=45*np.pi/180,
                                    observer_altitude_m=200000,
                                    cos_sza=0.6)
    viewing_geo.add_ray(ray)

wavel = np.arange(280.0, 800.0, 1)

atmosphere = sk.Atmosphere(model_geometry, config, wavelengths_nm=wavel)

sk.climatology.us76.add_us76_standard_atmosphere(atmosphere)

atmosphere['rayleigh'] = sk.constituent.Rayleigh()
atmosphere['ozone'] = sk.climatology.mipas.constituent("O3", sk.optical.O3DBM())

engine = sk.Engine(config, model_geometry, viewing_geo)

radiance_no_aerosol = engine.calculate_radiance(atmosphere)

Now we can add in aerosol, letโ€™s set the extinction to be a constant 1e-7 per metre until 30 km, and then zero above that. When we specify our Mie scatterer, we have to make sure to also specify any arguments of the particle size distribution (median_radius and mode_width) as a function of altitudes_m.

aero_ext = np.zeros(len(altitudes_m))
aero_ext[0:30] = 1e-7

atmosphere['aerosol'] = sk.constituent.ExtinctionScatterer(
    mie_db,
    altitudes_m=altitudes_m,
    extinction_per_m=aero_ext,
    extinction_wavelength_nm=745,
    median_radius=np.ones_like(aero_ext)*120,
    mode_width=np.ones_like(aero_ext)*1.6
)

radiance_with_aerosol = engine.calculate_radiance(atmosphere)

And we can plot the results

plt.plot(
    wavel, radiance_no_aerosol['radiance'].isel(los=0, stokes=0),
    wavel, radiance_with_aerosol['radiance'].isel(los=0, stokes=0)
)
[<matplotlib.lines.Line2D at 0x785b540e87a0>,
 <matplotlib.lines.Line2D at 0x785b540e8a10>]
../_images/a36d58bd477d9fe9a2b8210025ce5cc27d52ad80e1863682a7cac135bb96f36c.png

If we look at the dataset with aerosol included,

radiance_with_aerosol
<xarray.Dataset> Size: 7MB
Dimensions:                   (wavelength: 520, los: 4, stokes: 1,
                               aerosol_altitude: 66, altitude: 66,
                               ozone_altitude: 50)
Coordinates:
  * wavelength                (wavelength) float64 4kB 280.0 281.0 ... 799.0
  * stokes                    (stokes) <U1 4B 'I'
Dimensions without coordinates: los, aerosol_altitude, altitude, ozone_altitude
Data variables:
    radiance                  (wavelength, los, stokes) float64 17kB 0.000584...
    wf_aerosol_extinction     (aerosol_altitude, wavelength, los, stokes) float64 1MB ...
    wf_temperature_k          (altitude, wavelength, los, stokes) float64 1MB ...
    wf_aerosol_mode_width     (aerosol_altitude, wavelength, los, stokes) float64 1MB ...
    wf_pressure_pa            (altitude, wavelength, los, stokes) float64 1MB ...
    wf_aerosol_median_radius  (aerosol_altitude, wavelength, los, stokes) float64 1MB ...
    wf_specific_humidity      (altitude, wavelength, los, stokes) float64 1MB ...
    wf_ozone_vmr              (ozone_altitude, wavelength, los, stokes) float64 832kB ...

We can see that the model has also calculated the derivatives with respect to the particle size parameters, median_radius and mode_width.

Freezing Distribution Parameters#

The above method of creating Mie scattering databases is useful if you want to run calculations over a wide range of the particle size distribution parameters (median_radius and mode_width) in our example. But there are many applications where some or all of these values are fixed and known. In this case we can โ€œfreezeโ€ the parameters of the distribution,

distribution = sk.mie.distribution.LogNormalDistribution().freeze(median_radius=160, mode_width=1.6)

Now when we look at the args property of the distribution,

distribution.args()
[]

We see that it is empty, the distribution is not a function of any parameters because they are all frozen.

We can construct the Mie scattering database and inspect it,

mie_db = sk.database.MieDatabase(
    distribution,
    refrac,
    wavelengths_nm=np.arange(270, 1000, 50.0),
)

mie_db.load_ds()
<xarray.Dataset> Size: 109kB
Dimensions:        (wavelength_nm: 15, cos_angle: 128, legendre: 64)
Coordinates:
  * wavelength_nm  (wavelength_nm) float64 120B 270.0 320.0 ... 920.0 970.0
  * cos_angle      (cos_angle) float64 1kB -0.9993 -0.9964 ... 0.9999 1.0
    distribution   int64 8B ...
Dimensions without coordinates: legendre
Data variables: (12/13)
    xs_scattering  (wavelength_nm) float64 120B ...
    xs_total       (wavelength_nm) float64 120B ...
    xs_absorption  (wavelength_nm) float64 120B ...
    p11            (wavelength_nm, cos_angle) float64 15kB ...
    p12            (wavelength_nm, cos_angle) float64 15kB ...
    p33            (wavelength_nm, cos_angle) float64 15kB ...
    ...             ...
    lm_a1          (wavelength_nm, legendre) float64 8kB ...
    lm_a2          (wavelength_nm, legendre) float64 8kB ...
    lm_a3          (wavelength_nm, legendre) float64 8kB ...
    lm_a4          (wavelength_nm, legendre) float64 8kB ...
    lm_b1          (wavelength_nm, legendre) float64 8kB ...
    lm_b2          (wavelength_nm, legendre) float64 8kB ...

We can see that the Mie properties are no longer a function of the particle size distribution parameters since they are assumed to be known. In addition when constructing the database we did not have to specify them.

Similarly, to add aerosol in to our SASKTRAN2 atmosphere we no longer have to specify the frozen particle size parameters,

atmosphere['aerosol'] = sk.constituent.ExtinctionScatterer(
    mie_db,
    altitudes_m=altitudes_m,
    extinction_per_m=aero_ext,
    extinction_wavelength_nm=745,
)

In this example we froze all of the distribution parameters, but we could have frozen only one of them

sk.mie.distribution.LogNormalDistribution().freeze(mode_width=1.6).args()
['median_radius']

Using this distribution we would have to specify median_radius when constructing the database and when adding scatterers to the atmsophere.