implementing geometrical optics for rough interfaces within ice column

[HeWhoLikesWaffles]

My colleagues and I are interested in implementing a snow pack over rough ice interface, at least for the top layer of an ice column. Our goal is to model polarimetric brightness temperature at several microwave frequencies.

As a first cut, I've attempted to implement GO to describe the rough ice interface via:

<import smrt libraries>
from smrt.core.interface import make_interface

<boiler plate variables, build snow pack, ice column>

rough_interface = make_interface("geometrical_optics",\
                                 mean_square_slope=mean_square_slope)
								 
ice_column = make_ice_column(ice_type              = ice_type,
                             thickness             = thickness,
                             temperature           = temperature,
                             microstructure_model  = "sticky_hard_spheres",
                             brine_inclusion_shape = "spheres", 
                             salinity              = salinity, 
                             brine_volume_fraction = .02,
                             radius                = radius,
                             stickiness            = stickiness,
                             density               = density,
                             add_water_substrate   = "ocean",
                             interface = rough_interface
                            )

<rest of code>

When executing this code I saw some divergent results. It was only after doing so that I noticed the following admonition within /smrt/interface/__init.py__:

""" This module contains different type of boundary conditions between the layers.
Currently only flat interfaces are implemented. 

.. admonition::  **For developers**

    All the different type of interface must defined the methods: `specular_reflection_matrix` and `coherent_transmission_matrix`.

    It is currently not possible to implement rough interface, a (small) change is needed in DORT. Please contact the authors.

"""

QUESTIONS:
1.) What "small change" is needed for the DORT solver? Would this be something we could help with?
2.) There is a warning message originating from geometrical_optics.py:

to be optimised

Does this refer to a speed optimization or a convergence issue?

[JulienMeloche]

1. I'm pretty sure the __init__.py is outdated and the small change have been done to DORT since we can now run interface with IEM and GO.

2. I also notice the print(" to be optimised") in my work. I'm not sure what it is so I will let @ghislainp answer....

For your divergent results, are you sure you are within the domain of GO?

from geometrical_optics.py
This approximation is suitable for surface with roughness much larger than the roughness scales, typically k*s >> 1 and k*l >> 1, where s the rms heigth and l the correlation length. The precise validity range must be investigated by the user, this code does not raise any warning.

[HeWhoLikesWaffles]

Thanks for your reply. The surfaces we're considering have ks>>1 & kl>>1, so should be well within the GO limit. However, in the example usage within test_geometrical_optics.py and within the tutorials page, geometrical_optics is only called by specifying the mean square slope (MSS) of the surface. Is there a way to call GO by explicitly passing ks and kl that I am missing?

I've tried using MSS values ranging from 0.01 to 2.

[JulienMeloche]

No you need to pass mss with GO. IEM needs s and l separately. In the end you are still using the same parameter from mean_square_slope = 2*s**2/l**2.

If you still have unexpected results... The roughness might be to large for your application and sensor... I know sea ice roughness is usually in the mm (Landy et al 2014). Maybe consider using IEM with smaller roughness. But at this point, I doubt this is an issue with the code in SMRT...

[HeWhoLikesWaffles]

I've found the problem. Whether it's an issue within SMRT is open to debate, but there seem to be cases where the GO routine does not work well if more than one interface within an ice column is roughened, or if any interface other than the top-most interface is roughened. I believe adding a new error or warning message is warranted. Considering the following:

# import libraries
import numpy as np

# import smrt and related functions
from smrt import make_ice_column, make_snowpack, make_model, sensor
from smrt import PSU # g/kg -> kg/kg conversion
from smrt.atmosphere.simple_isotropic_atmosphere import SimpleIsotropicAtmosphere
from smrt.core.interface import make_interface # import lib for roughness

# ice inputs
l           = 2                             # number ice layers
thickness   = np.array([.5,.5])             # ice thickness in m
radius      = np.array([.0001,.0001])       # particle radius
stickiness  = np.array([0.3, 0.3])          # 'tau' 
temperature = np.array([260, 265])          # ice temperature in K
salinity    = np.array([5.1, 5.4])*PSU      # ice salinity profile in kg/kg
density     = np.array([0.924,0.924])       # ice density profile in kg/m^3
mean_square_slope = 0.1                       # MSS (unitless)

# atmosphere
atmos = SimpleIsotropicAtmosphere(10., 3., 0.90)

# create a first-year sea ice column:
ice_type = 'firstyear' # first-year or multi-year sea ice

rough_interface = make_interface("geometrical_optics",\
                                 mean_square_slope=mean_square_slope)

ice_column = make_ice_column(ice_type              = ice_type,
                             thickness             = thickness, # meters
                             temperature           = temperature,
                             microstructure_model  = "sticky_hard_spheres",
                             brine_inclusion_shape = "spheres", 
                             salinity              = salinity, 
                             brine_volume_fraction = .02,
                             radius                = radius,
                             stickiness            = stickiness,
                             density               = density,
                             add_water_substrate   = "ocean" #see comment below
                            )
ice_column.interfaces[0]=rough_interface
#ice_column.interfaces[1]=rough_interface

# snow inputs:
l_s           = 2                            # number of layers
thickness_s   = np.array([.12,.12])          # thickness per layer
density_s     = np.array([216,100])          # density profile in kg/m^3
radius_s      = np.array([.0004, .00025])    # particle radius
stickiness_s  = np.array([0.20, 0.20])       # 'tau' (stickiness)
temperature_s = np.array([253, 252])         # temperature

# create the snowpack
snowpack = make_snowpack(thickness            = thickness_s,
                         microstructure_model = "sticky_hard_spheres",
                         density              = density_s,
                         temperature          = temperature_s,
                         radius               = radius_s,
                         stickiness           = stickiness_s)

snowpack.atmosphere = atmos# from test

#add snowpack on top of ice column:
medium = snowpack + ice_column

# create geometric and EM parameters
theta        = 55           # Earth incidence angle
freq         = 6E9     # range of frequencies  (Hz)
sensor=sensor.passive(freq, theta)

# make model
n_max_stream = 256           # TB calcis more accurate if # of streams increased
                            # (currently: default = 32); 
                            # needs to be increased when using > 1 snow layer 
m = make_model("iba", "dort", rtsolver_options ={"n_max_stream": n_max_stream},
               emmodel_options=dict(dense_snow_correction="auto"))

# run the model for bare sea ice:
res1 = m.run(sensor, ice_column)
# run the model for snow-covered sea ice:
res2 = m.run(sensor, medium)

# print TBs at horizontal and vertical polarization:
print(res1.TbH(), res1.TbV())
print(res2.TbH(), res2.TbV())

As-is, the code above yields:

207.19735077300595 245.69465007105475
200.60817061769922 227.73977475803218

A reasonable result for H and V TBs. But, un-commenting line 42 (ice_column.interfaces[1]=rough_interface) yields:

5427.083559428849 4236.1671434738155
-16228.632555060074 -12573.256466463112

Defining only the bottom interface with GO yields a similarly non-physical TB prediction. Increasing n_max_streams yields higher-magnitude TBs, which is what led me to say the result diverges. Something similar happens if I only describe the bottom-most interface with GO (and leave the remaining interfaces flat).

RECOMMENDATION: Create a warning, error message, or admonition which advises users against using GO for layers other than the top-most layer.

[ghislainp]

note the the ice density values are in g/cm3, they must be in kg/m3