the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 06 May 2026
A new Monte Carlo treatment of multiple scattering of light by black carbon aggregates with varying levels of compactness
Abstract. We use a specially designed Monte Carlo (MC) model that includes individual Mueller matrix selection in each scattering event to conduct a comprehensive theoretical investigation of the intensity and polarization of light multiply scattered by an atmospheric layer embedded with black carbon (BC) particles with varying compactness. We find that the spatial distributions of the degree of linear and circular polarization are sensitive to BC particle shape, and that depolarization is more dominant for spherical BC particles than for fractal aggregates. We also find that the layer transmittance is highest for the more spherical BC particle shapes and lowest for the extended fractal aggregates, with an overall difference of approximately 10 % relative to the incident intensity between the highest and lowest values of transmittance. We find that the reflectance exhibits a similar tendency, and that correspondingly, the layer absorptance is lowest for the more spherical BC shapes and highest for the extended fractal aggregates. In addition, we compare the results obtained from our MC simulations with those obtained with the delta-Eddington approximation for the same optical thickness, single scattering albedo, asymmetry factor, and incident zenith angle. We find that there is a relatively small difference between our results and the delta-Eddington results with respect to the layer transmittance and the layer absorbance, but that the layer reflectance obtained with our MC simulations is up to 50 % lower than that obtained with the delta-Eddington approximation. These results could have important implications for radiative forcing estimations, climate modeling, and remote sensing implementations.
- Preprint
(5032 KB) - Metadata XML
- BibTeX
- EndNote
Status: open (until 17 Sep 2026)
-
RC1: 'Comment on egusphere-2026-1619', Anonymous Referee #1, 12 Aug 2026
reply
-
AC1: 'Response to the comments of Referee #1', Ynon Hefets, 16 Sep 2026
reply
Response to the comments of Referee #1
We thank the anonymous referee for these important comments, and we address each point below.
Comment #1: My main concern is the treatment of particle orientation in the proposed “individual Mueller matrix selection” method. Rather than using a Mueller matrix averaged over particle orientations, the model randomly selects a particular particle realization and orientation for each scattering event and applies the corresponding Mueller matrix. The authors argue that this treatment is more representative of the physical interaction between a photon and an individual particle.
Based on my understanding, if the individual-event method and the conventional pre-averaged method represent exactly the same random particle population, I would expect the two methods to converge toward the same ensemble-mean result when the orientation-dependent extinction, scattering, and phase-matrix properties are weighted consistently. However, the manuscript reports that the two approaches give systematic differences, which can reach to even 10% as shown in Appendix A. The difference might represent genuine physics, but it might instead arise because the authors are not weighting the orientations in exactly the same way as the conventional orientation average.
In particular, the authors should demonstrate that the orientation-sampling procedure gives the correct statistical weighting to particles with different orientation-dependent optical properties. Otherwise, the difference you see compared with the standard orientation-averaged method may just come from giving different particle orientations the wrong statistical weights.
Response: We agree with the referee that generally speaking, if the two methods are intended to represent the same random particle population, then orientation selection and a rigorously implemented pre-averaged treatment should converge to the same ensemble-mean result. What we should have presented more clearly are the following three points. (1) We are interested in a general framework that can in principle handle cases in which the distribution of particle shapes (in the present work, aggregates) is not necessarily mirror symmetric. Individual nonspherical aggregates can possess shape chirality, and the population need not contain mirror-related aggregate forms in exactly equal statistical abundance. Even if the underlying population is mirror symmetric in expectation, a finite sparse realization of an optically thin layer need not contain equal numbers of mirror-related forms. In such cases, pre-averaging can still be used, but it must be performed at the level of the full dimensional vector operators. A dimensional scattering matrix averaged over realization and orientation must be constructed, and the corresponding full extinction matrix must likewise be averaged over realization and orientation rather than replaced by a scalar mean extinction cross section. These averaged matrices can retain elements that vanish for a perfectly mirror-symmetric ensemble. In this general case, the integrated scattering cross section can also depend on the incident polarization state. A rigorous pre-averaged Monte Carlo treatment must therefore allow the full Stokes vector to influence both scattering and propagation between scattering events. Our explicit-selection method, on the other hand, naturally retains the mutually correlated extinction, scattering strength, scattering matrix, handedness, and polarization response for each selected realization and orientation. (2) Our selection method additionally retains information about the variances and distributions of the Stokes-vector components and polarization states across photon histories, as well as correlations between particle handedness/orientation and the resulting polarization. These trajectory-level statistics are removed when the particle ensemble is pre-averaged. (3) In still more general situations (for example, when the orientations of particles encountered at successive interactions are correlated), the conventional independently pre-averaged radiative-transfer description would no longer be sufficient. Our explicit-selection method, on the other hand, can be extended naturally to such cases by drawing each encountered orientation from the appropriate spatially varying and conditional or joint orientation distribution, rather than treating successive orientations as independent draws from a uniform distribution.
In the revised manuscript, we will remove the phrasing “more representative” and instead state these three more specific and more accurately phrased points clearly. In addition, we will remove the results currently shown in Appendix A (in which the averaging requirements stated in point 1 above were not fulfilled in the pre-average and additionally the weighting of the orientation selection was not consistent, though proper and consistent weighting was employed in all of the other simulations with orientation selection in the manuscript) and replace them with a demonstration of point 2, namely how our method retains information on the distribution of polarization states across photon histories. In addition, we address the referee’s specific additional points below.
The manuscript states that one of 102 orientations is selected uniformly, followed by one of 360 azimuthal rotations. How are the 102 orientations distributed in three-dimensional orientation space, and do they provide an unbiased representation of randomly oriented particles?
Response: Please refer to our detailed response in the supplementary material.
How is the probability of selecting a particle orientation determined? Is every orientation equally probable, or should orientations be weighted according to orientation-dependent extinction/scattering cross section?
Response: For the cases examined this study, we consider the physical orientation distribution to be uniform over rotations, i.e., we consider each physical orientation to be equally probable, similar to Mishchenko's et al. (2000) formal definition of randomly oriented nonspherical particles. It is true that the probability of interaction of a photon with a given orientation is also proportional to that orientation’s extinction cross section relative to the ensemble-averaged extinction cross section. However, we found that, unlike the differences among the orientation-dependent Mueller matrices and among the orientation-dependent scattering cross sections, the differences among the orientation-dependent extinction cross sections are negligible (with relative differences <0.03), and thus we take the probability of interaction of a photon with a given orientation to also be uniform over rotations. We will clarify this very important point in the revised manuscript.
Why is the phase function for each selected orientation normalized separately? Please show mathematically that this normalization preserves the correct ensemble scattering probability.
Response: Please refer to our detailed response in the supplementary material.
What should happen theoretically if both methods describe exactly the same randomly oriented particle population? If the two methods should be equivalent, why does Appendix A show a difference? If they should not be equivalent, what physical assumption makes them different?
Response: Please refer to our detailed response to the referee’s general comment at the top of this document. As we stated there, in the revised manuscript, we will remove the phrasing “more representative” and instead state the three more specific and more accurately phrased points written there clearly. In addition, we will remove the results currently shown in Appendix A (in which the averaging requirements stated in point 1 above were not fulfilled in the pre-average and additionally the weighting of the orientation selection was not consistent, though proper and consistent weighting was employed in all of the other simulations with orientation selection in the manuscript) and replace them with a demonstration of how our method retains information on the distribution of polarization states across photon histories.
Comment #2: The main novelty of the manuscript concerns multiple scattering, but the comparison between the new individual-selection method and the conventional orientation-averaged method in Appendix A is limited to single scattering and one realization of the most extended aggregate. This makes it difficult to evaluate how much the proposed treatment actually affects the multiple-scattering results presented in the main text. I suggest repeating a representative subset of the full Monte Carlo calculations using both methods. For example, comparisons for one compact and one extended aggregate at an optical thickness of approximately 1, considering unpolarized, linearly polarized, and circularly polarized incident light, would be useful. Comparing quantities such as transmittance, reflectance, and selected polarization distributions would provide a much clearer demonstration of the practical impact of the proposed method.
Response: Given that we agree with the referee that if the two methods are intended to represent the same random particle population, then orientation selection and a rigorously implemented pre-averaged treatment should converge to the same ensemble-mean result (again, refer to our detailed response to the referee’s general comment at the top of this document), conducting additional comparison (with appropriate weighting) over additional cases would not demonstrate the novelty of our method. Instead, as stated above, we will remove the results currently shown in Appendix A and replace them with a demonstration of how our method retains information on the distribution of polarization states across photon histories.
In addition, we will emphasize better that the results themselves with the explicit anisotropic shapes and unique polarization patterns over the various cases are also part of the novelty of the work and not solely our individual Mueller matrix selection method.
Comment #3: What is the justification of using Delta-Eddington method as a validation benchmark. I am not fully convinced that the delta-Eddington approximation provides an adequate benchmark for validating the Monte Carlo model, since delta-Eddington is itself an approximate radiative transfer treatment. The polarized-Monte-Carlo methodology cited by the authors has historically been compared against adding-doubling solutions, and such comparisons can achieve approximately percent-level agreement. Established polarized adding-doubling approaches provide a natural independent benchmark.
Response: This is an important clarification. Our comparison to delta-Eddington is not intended as a validation benchmark. Rather, it is intended to show how the more explicit Monte Carlo solution differs from a typical approximate solution often employed in large-scale (climate or global circulation) models. We mentioned this on lines 230-234 of the original submitted manuscript, but based on the referee’s comment, we will emphasize it better in the revised manuscript.
Comment #4: Please carefully check the manuscript to ensure that all symbols and variables are clearly defined when they first appear and are used consistently throughout the text.
Response: As suggested, we went back over the text to check that all symbols and variables are defined and used consistently. We indeed found that the parameters q_(n,j), u_(n,j), and v_(n,j) in Eq. 6 were not defined explicitly in the text. In fact, we realized that they are not used anywhere after Eq. 6, and we have therefore decided to remove them. All other symbols and variables should be presented satisfactorily.
List of references that do not appear in the manuscript
Yurkin, M. A., & Hoekstra, A. G. (2011). The discrete-dipole-approximation code ADDA: Capabilities and known limitations. Journal of Quantitative Spectroscopy and Radiative Transfer, 112(13), 2234-2247, https://doi.org/10.1016/j.jqsrt.2011.01.031.
-
AC1: 'Response to the comments of Referee #1', Ynon Hefets, 16 Sep 2026
reply
-
RC2: 'Comment on egusphere-2026-1619', Anonymous Referee #2, 08 Sep 2026
reply
The paper has relevance to the modeling of radiative forcing by black carbon particles in the atmosphere. The writing is clear, the review of previous work is good, and the model and objectives are, for the most part, well defined.
My sole issue with this paper regards the use of a procedure, in the Monte Carlo algorithm for solving the RTE, which relies on the sampling of scattering directions using a distribution function corresponding to fixed orientation of the aggregate particle. The orientation of the particle from which this scattering direction is obtained is, in itself, randomly sampled from a discrete set of orientations. The authors submit that this approach is more representative (or more accurate) than that based upon the sampling from an orientation--averaged distribution, and they demonstrate in the appendix that the two approaches yield small yet significant differences. I am at a loss to see how there can be any difference from a physical point of view. If their RTE model is based upon random orientations + uncorrelated orientations among different particles + macroscopic conditions (the number of particles in a layer -> infinity) -- which I believe it is -- the only relevant features of the particle scattering characteristics are those corresponding to the orientation--averaged properties. Given this, the results shown in the appendix are, to me, more indicative of an error, inconsistency, or lack-of-convergence in their numerical method as opposed to some sort of improved representation of physics. A specific point, regarding these results (paragraph at line 540), is: what is the distinction between the 10^7 photons for each simulation and the 10 simulations? Would this be the same as 10*10^7 photons for a single simulation? The authors need to rectify this, either by demonstrating analytically that their method offers distinct numerical advantages and/or a better representation of the physics (which, personally, I don't see how it can), or by reverting to a numerical method in which scattering angles are sampled from orientation--averaged distribution functions.
Citation: https://doi.org/10.5194/egusphere-2026-1619-RC2 -
AC2: 'Response to the comments of Referee #2', Ynon Hefets, 16 Sep 2026
reply
Response to the comments of Referee #2
The paper has relevance to the modeling of radiative forcing by black carbon particles in the atmosphere. The writing is clear, the review of previous work is good, and the model and objectives are, for the most part, well defined.
Response: We thank the referee for this positive statement.
My sole issue with this paper regards the use of a procedure, in the Monte Carlo algorithm for solving the RTE, which relies on the sampling of scattering directions using a distribution function corresponding to fixed orientation of the aggregate particle. The orientation of the particle from which this scattering direction is obtained is, in itself, randomly sampled from a discrete set of orientations. The authors submit that this approach is more representative (or more accurate) than that based upon the sampling from an orientation--averaged distribution, and they demonstrate in the appendix that the two approaches yield small yet significant differences. I am at a loss to see how there can be any difference from a physical point of view. If their RTE model is based upon random orientations + uncorrelated orientations among different particles + macroscopic conditions (the number of particles in a layer -> infinity) -- which I believe it is -- the only relevant features of the particle scattering characteristics are those corresponding to the orientation--averaged properties. Given this, the results shown in the appendix are, to me, more indicative of an error, inconsistency, or lack-of-convergence in their numerical method as opposed to some sort of improved representation of physics… The authors need to rectify this, either by demonstrating analytically that their method offers distinct numerical advantages and/or a better representation of the physics (which, personally, I don't see how it can), or by reverting to a numerical method in which scattering angles are sampled from orientation--averaged distribution functions.
Response: The referee brings up an important point that is in a similar vein to the main comment of Referee #1, and we thus include the same response here. What we should have presented more clearly are the following three points. (1) We are interested in a general framework that can in principle handle cases in which the distribution of particle shapes (in the present work, aggregates) is not necessarily mirror symmetric. Individual nonspherical aggregates can possess shape chirality, and the population need not contain mirror-related aggregate forms in exactly equal statistical abundance. Even if the underlying population is mirror symmetric in expectation, a finite sparse realization of an optically thin layer need not contain equal numbers of mirror-related forms. In such cases, pre-averaging can still be used, but it must be performed at the level of the full dimensional vector operators. A dimensional scattering matrix averaged over realization and orientation must be constructed, and the corresponding full extinction matrix must likewise be averaged over realization and orientation rather than replaced by a scalar mean extinction cross section. These averaged matrices can retain elements that vanish for a perfectly mirror-symmetric ensemble. In this general case, the integrated scattering cross section can also depend on the incident polarization state. A rigorous pre-averaged Monte Carlo treatment must therefore allow the full Stokes vector to influence both scattering and propagation between scattering events. Our explicit-selection method, on the other hand, naturally retains the mutually correlated extinction, scattering strength, scattering matrix, handedness, and polarization response for each selected realization and orientation. (2) Our selection method additionally retains information about the variances and distributions of the Stokes-vector components and polarization states across photon histories, as well as correlations between particle handedness/orientation and the resulting polarization. These trajectory-level statistics are removed when the particle ensemble is pre-averaged. (3) In still more general situations (for example, when the orientations of particles encountered at successive interactions are correlated), the conventional independently pre-averaged radiative-transfer description would no longer be sufficient. Our explicit-selection method, on the other hand, can be extended naturally to such cases by drawing each encountered orientation from the appropriate spatially varying and conditional or joint orientation distribution, rather than treating successive orientations as independent draws from a uniform distribution.
In the revised manuscript, we will remove the phrasing “more representative” and instead state these three more specific and more accurately phrased points clearly. In addition, we will remove the results currently shown in Appendix A (in which the averaging requirements stated in point 1 above were not fulfilled in the pre-average and additionally the weighting of the orientation selection was not consistent, though proper and consistent weighting was employed in all of the other simulations with orientation selection in the manuscript) and replace them with a demonstration of point 2, namely how our method retains information on the distribution of polarization states across photon histories.A specific point, regarding these results (paragraph at line 540), is: what is the distinction between the 10^7 photons for each simulation and the 10 simulations? Would this be the same as 10*10^7 photons for a single simulation?
Response: The 10^7 photons are initialized identically in each of the 10 simulations. The idea was to repeat the simulation 10 times to verify that 10^7 photons are sufficient, such that the differences among the 10 simulations are negligible. Indeed, we found that the relative differences in the number of photons scattered into each possible scattering angle over the 10 simulations are <0.01. At the same time, as the reviewer suggests, it would be equivalent to regard the 10 simulations of 10^7 photons each as a single simulation with 10^8 photons.
Citation: https://doi.org/10.5194/egusphere-2026-1619-AC2
-
AC2: 'Response to the comments of Referee #2', Ynon Hefets, 16 Sep 2026
reply
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 252 | 40 | 19 | 311 | 26 | 22 |
- HTML: 252
- PDF: 40
- XML: 19
- Total: 311
- BibTeX: 26
- EndNote: 22
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
This manuscript investigates polarized multiple scattering by black carbon aggregates with different morphologies using a Monte Carlo radiative-transfer model in which a specific Mueller matrix is selected for each scattering event. The topic is relevant to aerosol optics, polarized radiative transfer, and potentially remote sensing. The manuscript contains a substantial amount of numerical work and several interesting polarization patterns. However, I still have several concerns regarding the physical basis and validation of the proposed methodology. In my view, these issues should be addressed before the main conclusions can be fully supported. I therefore recommend major revision.
Comment #1: My main concern is the treatment of particle orientation in the proposed “individual Mueller matrix selection” method. Rather than using a Mueller matrix averaged over particle orientations, the model randomly selects a particular particle realization and orientation for each scattering event and applies the corresponding Mueller matrix. The authors argue that this treatment is more representative of the physical interaction between a photon and an individual particle.
Based on my understanding, if the individual-event method and the conventional pre-averaged method represent exactly the same random particle population, I would expect the two methods to converge toward the same ensemble-mean result when the orientation-dependent extinction, scattering, and phase-matrix properties are weighted consistently. However, the manuscript reports that the two approaches give systematic differences, which can reach to even 10% as shown in Appendix A. The difference might represent genuine physics, but it might instead arise because the authors are not weighting the orientations in exactly the same way as the conventional orientation average.
In particular, the authors should demonstrate that the orientation-sampling procedure gives the correct statistical weighting to particles with different orientation-dependent optical properties. Otherwise, the difference you see compared with the standard orientation-averaged method may just come from giving different particle orientations the wrong statistical weights.
Based on this concern, there are several additional questions listed below:
Comment #2: The main novelty of the manuscript concerns multiple scattering, but the comparison between the new individual-selection method and the conventional orientation-averaged method in Appendix A is limited to single scattering and one realization of the most extended aggregate. This makes it difficult to evaluate how much the proposed treatment actually affects the multiple-scattering results presented in the main text. I suggest repeating a representative subset of the full Monte Carlo calculations using both methods. For example, comparisons for one compact and one extended aggregate at an optical thickness of approximately 1, considering unpolarized, linearly polarized, and circularly polarized incident light, would be useful. Comparing quantities such as transmittance, reflectance, and selected polarization distributions would provide a much clearer demonstration of the practical impact of the proposed method.
Comment #3: What is the justification of using Delta-Eddington method as a validation benchmark. I am not fully convinced that the delta-Eddington approximation provides an adequate benchmark for validating the Monte Carlo model, since delta-Eddington is itself an approximate radiative transfer treatment. The polarized-Monte-Carlo methodology cited by the authors has historically been compared against adding-doubling solutions, and such comparisons can achieve approximately percent-level agreement. Established polarized adding-doubling approaches provide a natural independent benchmark.
Comment #4: Please carefully check the manuscript to ensure that all symbols and variables are clearly defined when they first appear and are used consistently throughout the text.