Inquiry about strong fluctuations in limb radiance profiles from Smart-G Monte Carlo simulations

[20713]

Dear Mustapha,
I am currently using the Smart-G radiative transfer model to simulate limb-scattered radiance for ozone retrieval. While the aerosol optical inputs appear physically reasonable, I am observing unexpectedly strong fluctuations (noise-like oscillations) in the simulated limb radiance profiles（See figure below）.

Below is a summary of my current configuration and key diagnostic values:

• Vertical grid: 100 layers, 1 km spacing
• Aerosol extinction cross sections (per particle):
  • [3.67441e-09 3.67464e-09 3.68171e-09 3.91627e-09 4.05730e-09 4.32481e-09]cm²
• Integrated aerosol optical depth (per wavelength):
  • [0.37159 0.37161 0.37233 0.39605 0.41031 0.43736]
• Layer-by-layer dtau (per wavelength):
  • ·min = [0, 0, 0, 0, 0, 0]
  • ·max = [0.14389 0.1439 0.14417 0.15336 0.15888 0.16936]
• Photon number: default (1e6 per sensor)
• Aerosol phase matrix: from full Mie calculation
• Even with these reasonable aerosol parameters, the resulting limb radiance profiles show large oscillations with altitude, which are much stronger than expected. These do not resemble physical atmospheric variability, and appear to be dominated by Monte Carlo statistical noise.
  My questions:

1. Is such strong fluctuation in limb radiance typical when single-layer dtau reaches ~0.15–0.18, especially for limb geometry?

2. Does Smart-G require significantly more photons (e.g., 5×–10× higher than 1e6) to reduce noise in limb-viewing simulations with non-negligible aerosol loading?

3. Are there recommended settings for:

• the photon count (nph),
• angular resolution (mu sampling),
• or vertical smoothing of aerosol profiles,
  to improve numerical stability of limb radiance?

4. Is there a known issue or best practice regarding the use of full Mie aerosol scattering matrices in Smart-G?

If needed, I can provide my codes or sample outputs for inspection.

Thank you very much for your help, and for providing such a powerful radiative transfer tool.
Best regards,
Bin Hou

The main portion of the code is as follows:

# ---- Run Monte Carlo RTM with SMART-G ----

    # Limb geometry: one Sensor per tangent height (scheme-1)
    HTOA = 120 #714# km
    RTER = 6371.0 # km

    def zt2thv(zt, RTER=RTER, HTOA=HTOA):
        return np.degrees(np.arcsin((RTER + zt) / (RTER + HTOA)))

    ZT = 0.5 * (z_grid[:-1] + z_grid[1:])

    sensors = [
        Sensor(
            POSX=0.0, POSY=0.0, POSZ=RTER + HTOA,
            THDEG=180.0 - zt2thv(zt),
            PHDEG=0.0,
            LOC='ATMOS', FOV=0.0, TYPE=0
        )
        for zt in ZT
    ]
   # print(180.0 - zt2thv(ZT))

    surface = LambSurface(ALB=Albedo_cst(alb))

    le = dict(
    th_deg=[sza_deg],
    phi_deg=[saa_deg],
)

    S = Smartg(back=True, pp=False, double=True)

    NPHOT_PER_SENSOR = 1e6
    NB = int(NPHOT_PER_SENSOR * len(sensors))  

    rt = S.run(
        wl=wav_nm,
        surf=surface,
        le=le,
        atm=atm,
        NBPHOTONS=NB,
        refraction=True,
        sensor=sensors,
        reflectance=False,  
        #SEED=-1,
        FFS=True,
        NF=1e4,  #number_of_inverse_phase_function_discretisations
        DEPO = 0.0279 #depolarization_factor 
    )

More samples:
[INFO] Sample 31/50: SZA=50.0°, SAA=80.0°, ALB=1.00, IRANA=5081, REF=0.181μm, SIG=1.629
[CHK] AEROSOL_ext_cross_sections : [1.26161e-09 1.27696e-09 1.32029e-09 1.48749e-09 1.52678e-09 1.53276e-09]
sum dtau_aer: [0.13895 0.14064 0.14542 0.16383 0.16816 0.16882]
Per-wavelength dtau_aer min: [0. 0. 0. 0. 0. 0.]
Per-wavelength dtau_aer max: [0.05464 0.05531 0.05718 0.06442 0.06613 0.06639]

[INFO] Sample 45/50: SZA=15.0°, SAA=20.0°, ALB=0.10, IRANA=18662, REF=0.254μm, SIG=1.535
[CHK] AEROSOL_ext_cross_sections : [2.73014e-09 2.73387e-09 2.80005e-09 3.05651e-09 3.20762e-09 3.38168e-09]
sum dtau_aer: [0.41108 0.41164 0.42161 0.46022 0.48298 0.50918]
Per-wavelength dtau_aer min: [0. 0. 0. 0. 0. 0.]
Per-wavelength dtau_aer max: [0.14637 0.14657 0.15012 0.16387 0.17197 0.1813 ]

[INFO] Sample 50/50: SZA=30.0°, SAA=160.0°, ALB=0.40, IRANA=21052, REF=0.092μm, SIG=1.459
[CHK] AEROSOL_ext_cross_sections : [5.48855e-10 5.57321e-10 5.73128e-10 5.16408e-10 4.60935e-10 4.03870e-10]
sum dtau_aer: [0.15458 0.15696 0.16141 0.14544 0.12981 0.11374]
Per-wavelength dtau_aer min: [0. 0. 0. 0. 0. 0.]
Per-wavelength dtau_aer max: [0.05998 0.06091 0.06263 0.05643 0.05037 0.04414]

[Mustapha-hygeos]

Hi Bin Hou,

"Is such strong fluctuation in limb radiance typical when single-layer dtau reaches ~0.15–0.18, especially for limb geometry?"

I tried to repeat the same simulations but using desert aerosols, which have a very strong forward peak (worst aer scenario), and I don't reproduce abnormal strong fluctuations. For even higher tau than 0.18 we produce stable limb radiances:

And here the optical prop:

"Does Smart-G require significantly more photons (e.g., 5×–10× higher than 1e6) to reduce noise in limb-viewing simulations with non-negligible aerosol loading?"

Not for your case. It seems to be an issue with the input optical properties, may be from the phase matrices?

"Are there recommended settings for:

• the photon count (nph),
• angular resolution (mu sampling),
• or vertical smoothing of aerosol profiles,
  to improve numerical stability of limb radiance?"

You can evaluate the MC noise by setting the parameter stdev to true. If you set stdev=True make sure NBLOOP is sufficiently smaller than NBPHOTONS (at least 30 times smaller should be nice for a good variance estimate). Then you can adjust the number of photons.

However I don't think it will solve your issue since the problem seems more related to the phase matrices.

"Is there a known issue or best practice regarding the use of full Mie aerosol scattering matrices in Smart-G?"

You don't show how you incorporate the optical properties in SMART-G, but if you have used the class AerOPAC where the parameter phase is also used then you should know that a correction (also in the docstring) has been made very recently (in v1.1.2 or v1.1.3).

If you need to incorporate phase matrices we recommend to mimic the aerosol or cloud auxdata files (netcdf), then you specify your file(s) when using AerOPAC or Cloud classes. Do not forget to modify pfgrid parameter in AtmAFGL class if you have more than one aerosols (if you set pfgrid = z_grid you make sure to have at each layer a unique phase matrix).

Also another known issue is when using the local estimate with clouds, more specifically ice clouds with a very high forward peak in the phase function. In this case you must use at least 2 orders of magnitude more photons... variance reduction methods exists but not yet implemented in SMART-G. However, it is planned to implement such a variance reduction method next months.

Regards,
Mustapha

[20713]

Dear Mustapha,

Thank you very much for your previous reply and for your helpful explanations regarding SmartG. They were very useful for clarifying my understanding of the aerosol scattering implementation.

I am currently using SmartG for limb-scattering simulations with stratospheric aerosols. In my implementation, I do not use the built-in AerOPAC class. Instead, aerosol scattering properties are computed externally using Mie theory for spherical particles, and the resulting phase functions are ingested into SmartG via the prof_phases argument of the AtmAFGL class. I would like to briefly explain this workflow and kindly ask whether it is consistent with SmartG’s intended design.

In my setup, aerosol particles are assumed to be spherical liquid droplets with a single-mode lognormal particle size distribution. For each wavelength, an external Mie module provides the scattering cross section and the angle-resolved Mueller matrix:

mu= np.linspace(-1, 1, 201)
saer, Pp = eng.Mie_compute_wrapper(lam, mu, Mie, ref, sig, nargout=2)# Pp shape: (nlam, Nθ, 4, 4)

Following SmartG’s handling of spherical particles in the calculF() routine (called inside SmartG.run()), I extract only the four independent angular components required for spherical scattering and reorganize them as:

theta_deg = np.linspace(0, 180, len(mu))
P11 = Pp[il, :, 0, 0]
P22 = Pp[il, :, 1, 1]
P33 = Pp[il, :, 2, 2]
P43 = Pp[il, :, 3, 2]
P_stk_theta = np.stack([P11, P22, P33, P43], axis=0)  # shape = (4, Nθ)

These four-component phase functions are then stored as lookup tables with axes (stk, theta_atm) and passed to SmartG through:

# --- 1. Build LUT for a single wavelength ---
phaLUT = LUT(
    P_stk_theta,
    axes=[np.arange(4, dtype=np.int32), theta_deg],
    names=['stk', 'theta_atm']
)

# --- 2. Append LUT to the list ---
Aer_phase_list.append(phaLUT)

# --- 3. Pass to AtmAFGL via prof_phases ---
atm = AtmAFGL(
    ...,
    prof_phases=(iphase, Aer_phase_list)
)

As I understand from the SmartG source code, the particle type is automatically identified as spherical in the calculF() routine (invoked inside SmartG.run() to convert phase-function profiles into cumulative phase functions for Monte Carlo sampling) when the phase function contains four polarization components:

phase = profile['phase_atm'][ipha, :, :]
if (len(phase[:,0]) == 4): # spherical particle
    # parameters equally spaced in scattering probability
    # phase_H['p_P11'][idx, :] = interp1d(scum, phase[1,:])(z)  # I par P11
    # phase_H['p_P22'][idx, :] = interp1d(scum, phase[0,:])(z)  # I per P22
    phase_H['p_P11'][idx, :] = interp1d(scum, phase[0,:])(z)  # I par P11
    phase_H['p_P22'][idx, :] = interp1d(scum, phase[1,:])(z)  # I per P22
    phase_H['p_P33'][idx, :] = interp1d(scum, phase[2,:])(z)  # U P33
    phase_H['p_P43'][idx, :] = interp1d(scum, phase[3,:])(z)  # V P43
    phase_H['p_P44'][idx, :] = interp1d(scum, phase[2,:])(z)  # V P44= P33
    phase_H['p_ang'][idx, :] = interp1d(scum, angles)(z) # angle

    # parameters equally spaced in scattering angle [0, 180]
    phase_H['a_P11'][idx, :] = f1(angN)  # I par P11
    phase_H['a_P22'][idx, :] = f2(angN)  # I per P22
    phase_H['a_P33'][idx, :] = f3(angN)  # U P33
    phase_H['a_P43'][idx, :] = f4(angN)  # V P43
    phase_H['a_P44'][idx, :] = f3(angN)  # V P44=P33
else: # non spherical particle
    f5 = interp1d(angles, phase[4,:])
    f6 = interp1d(angles, phase[5,:])
......

Based on this logic, SmartG then applies the spherical-scattering branch in calculF().

Could you please let me know whether this approach—using externally computed spherical Mie phase functions instead of AerOPAC, and providing exactly these four components in this order—is fully consistent with SmartG’s expectations?
If there are any additional normalization or formatting requirements that I should be aware of, I would greatly appreciate your guidance.

Thank you very much for your time and for developing and maintaining SmartG.

Best regards,
Bin Hou

[20713]

Dear Mustapha,

Thank you very much for your previous reply and for your helpful explanations regarding SmartG. They were very useful for clarifying my understanding of the aerosol scattering implementation.

I am currently using SmartG for limb-scattering simulations with stratospheric aerosols. In my implementation, I do not use the built-in AerOPAC class. Instead, aerosol scattering properties are computed externally using Mie theory for spherical particles, and the resulting phase functions are ingested into SmartG via the prof_phases argument of the AtmAFGL class. I would like to briefly explain this workflow and kindly ask whether it is consistent with SmartG’s intended design.

In my setup, aerosol particles are assumed to be spherical liquid droplets with a single-mode lognormal particle size distribution. For each wavelength, an external Mie module provides the scattering cross section and the angle-resolved Mueller matrix:

mu= np.linspace(-1, 1, 201)
saer, Pp = eng.Mie_compute_wrapper(lam, mu, Mie, ref, sig, nargout=2)# Pp shape: (nlam, Nθ, 4, 4)
Following SmartG’s handling of spherical particles in the calculF() routine (called inside SmartG.run()), I extract only the four independent angular components required for spherical scattering and reorganize them as:

theta_deg = np.linspace(0, 180, len(mu))
P11 = Pp[il, :, 0, 0]
P22 = Pp[il, :, 1, 1]
P33 = Pp[il, :, 2, 2]
P43 = Pp[il, :, 3, 2]
P_stk_theta = np.stack([P11, P22, P33, P43], axis=0) # shape = (4, Nθ)
These four-component phase functions are then stored as lookup tables with axes (stk, theta_atm) and passed to SmartG through:

--- 1. Build LUT for a single wavelength ---

phaLUT = LUT(
P_stk_theta,
axes=[np.arange(4, dtype=np.int32), theta_deg],
names=['stk', 'theta_atm']
)

--- 2. Append LUT to the list ---

Aer_phase_list.append(phaLUT)

--- 3. Pass to AtmAFGL via prof_phases ---

atm = AtmAFGL(
...,
prof_phases=(iphase, Aer_phase_list)
)
As I understand from the SmartG source code, the particle type is automatically identified as spherical in the calculF() routine (invoked inside SmartG.run() to convert phase-function profiles into cumulative phase functions for Monte Carlo sampling) when the phase function contains four polarization components:

phase = profile['phase_atm'][ipha, :, :]
if (len(phase[:,0]) == 4): # spherical particle
# parameters equally spaced in scattering probability
# phase_H['p_P11'][idx, :] = interp1d(scum, phase[1,:])(z) # I par P11
# phase_H['p_P22'][idx, :] = interp1d(scum, phase[0,:])(z) # I per P22
phase_H['p_P11'][idx, :] = interp1d(scum, phase[0,:])(z) # I par P11
phase_H['p_P22'][idx, :] = interp1d(scum, phase[1,:])(z) # I per P22
phase_H['p_P33'][idx, :] = interp1d(scum, phase[2,:])(z) # U P33
phase_H['p_P43'][idx, :] = interp1d(scum, phase[3,:])(z) # V P43
phase_H['p_P44'][idx, :] = interp1d(scum, phase[2,:])(z) # V P44= P33
phase_H['p_ang'][idx, :] = interp1d(scum, angles)(z) # angle

# parameters equally spaced in scattering angle [0, 180]
phase_H['a_P11'][idx, :] = f1(angN)  # I par P11
phase_H['a_P22'][idx, :] = f2(angN)  # I per P22
phase_H['a_P33'][idx, :] = f3(angN)  # U P33
phase_H['a_P43'][idx, :] = f4(angN)  # V P43
phase_H['a_P44'][idx, :] = f3(angN)  # V P44=P33

else: # non spherical particle
f5 = interp1d(angles, phase[4,:])
f6 = interp1d(angles, phase[5,:])
......
Based on this logic, SmartG then applies the spherical-scattering branch in calculF().

Could you please let me know whether this approach—using externally computed spherical Mie phase functions instead of AerOPAC, and providing exactly these four components in this order—is fully consistent with SmartG’s expectations? If there are any additional normalization or formatting requirements that I should be aware of, I would greatly appreciate your guidance.

Thank you very much for your time and for developing and maintaining SmartG.

Best regards, Bin Hou

We reran the code following the original approach, and the results look normal so far. Of course, given your expertise, it would be greatly appreciated if you could help review and assess them.

[Mustapha-hygeos]

Dear Bin Hou,

Good to see that you were able to reproduce consistent results.

Concerning you previous approach, it's almost correct. If you provide only 4 phase components (only for spherical particles) then you have to provide (and in this oder) : P11, P21, P33 and P34. And if you provide 6 phase components (for spheroidal, but also work for spherical particles): P11, P21, P33, P34, P22 and P44. You have to ignore what is written in calculF function because this is after the conversion to Iper Ipar convention as in the book of S. Chandrasekhar (radiative transfer 1960). By default, in the calc method you have the parameter "conv_Iparper" set to True. So normally you don't have to convert anything.

But the safest way is as I already explained previously: you create an MLUT (if you are using the lut lib) and mimic the aerosol auxiliary data but with your own calculated optical properties, and save it. Then you can use AerOPAC as usually but by specifying the path of your created netcdf4 file.

Best regards,
Mustapha