Overall, it is a great paper with valuable reproducible benchmark results (to the best of my own knowledge, reported first time)  written by leading experts in RT numerical simulations. Below, I dare to provide some criticism, questions, and comments, the number of which may seem high. But it only indicates that the paper is interesting and reports new, practically important, and reproducible results. I'll be happy if the authors will address some of my comments, but will not argue if they skip them all.

• Section 2 “Model Descriptions”

For each model, please indicate if it is publicly available and provide weblinks. E.g., MYSTIC libradtran.org is publicly available without limitations, while distribution of SCIATRAN may require registration, etc.

For each model, please indicate parameters that influence speed and accuracy. E.g., for Monte Carlo it is the number of photons, for successive orders it is the number of scattering orders and Gauss ordinates to integrate over angle (Fourier for  azimuth?), for Gauss-Seidel that would be the number of iterations and Gauss nodes (and Fourier moments?), etc. Line 131: SASKTRAN-HR “has been configured to use the number of diffuse profiles required to obtain approximately 0.2% accuracy…” – what parameters need to be adjusted? Are these parameters come as input (hence can be changed by a user), or they are hard-coded and can be changed only with developers help?

For each model, it would be very helpful to provide a ready to use input (with short instructions), so that an interested user could download and quickly run a model and be able to independently reproduce the reported results;

• Numerical results (benchmarks) are the most useful part of the paper. However, it is not clear from the beginning how to get the numbers. A simple python driver to read the provided netcdf file would also be helpful due to the amount of the data. I would recommend providing structure (tree) for the netcdf file (e.g. in Appendix).
• Different instruments are discussed in the Introduction; however, their measurment accuracy is not indicated. The authors are talking about 0.1, 0.5, 1% error in simulations (e.g. Abstract & Conclusion). However, it is not clear if this level of error is too good or just right or insufficient. Since ratios of intensities are often used, one should keep in minds that sufficient level in simulated intensity may or may not be sufficient for the ratio.
• Fig.2: definition of angles is unclear. Specifically, is the zero relative azimuth correspond to forward or backward scattering? Same for the zenith angle ‘theta’ (bidirectional arrows are confusing). It would be helpful to show observation of another tangent point (I believe the observer will change location, so that local normal at the tangent point would not move);
• Across the text: does the multiple scattering include the first one? Or multiple literally means second and higher?
• Line 230: 11 wavelengths. Why all 11 are needed? It would make since to test the codes for the shortest band, longest band, and something in between (like, Lambertian surface is defined for black, gray and white cases – only three). If an RT code works for 3 mentioned bands, why it would not work for the rest? Does any of the mentioned 11 bands have some interesting features (e.g., peculiarities in O3-profile), not presenting in other bands? Leaving only 3 necessary bands would make it easier to represent the results while pursuing “the purpose of brevity” (line 286).

• Section 2.2.5. is important but very confusing… First, in the parenting Section 2.2, two approaches are discussed: forward and backward. What about local estimation - is it similar to backward (excatly the same... completely different...)? Further, MYSTIC “used the backward MC method” (line 171). However, in lines 210-214, only Siro and SASKTRAN-MC are mentioned as based on  backwards technique, while MYSTIC is mentioned one paragraph above (different technique?)   Condition “not if scattering actually occurs” (line 211) is unclear. Neither is “a desired noise floor” (line 213) term – does that mean level of MC noise? Finally, description of the SMART-G feature (line 215) seems too brief to understand – please either elaborate or reformulate or maybe refer to the SMART-G section. Overall, differences between MC techniques deserves a separate paper, and it is great the authors decided to summaries them in one paragraph. But such a paragraph should be written with extreme care.

• A whole paragraph around line 255, where polarization discussed: what is the frame of reference that defines the Q & U components?
• For aerosol scenario, it would be helpful to plot the phase matrix vs scattering angle and tabulate expansion moments for deterministic RT codes.
• Section 4.3: how percentile deviation for DoLP is defined? At neutral points, DoLP -> 0 and delta_DoLP would be Inf. “LPO” – orientation wrt what reference plane?
• Table 4: Timing – does the scaled runtime account for different number of CPU cores in addition to difference in CPU themselves? E.g., if one runs SASKTRAN on 6 cores (instead of 12 as indicated), will the result be 2 x 1909 ~= 4000?

*********************************************************************

Minor Comments:

Title: I’ve never met the word “scalarial” and would suggest to use “vector” instead of “vectorial”. The word “vector” is actually used in the text “vector test cases” (line 7), “vector modes” (line 11), and other places in the paper including list of references; See also line 102 for “… vectorial radiance…”.

In the Abstract, “fully spherical” sounds confusing and should be either explained or reduced to “spherical”, while the word “fully” should be introduced and explained later.

Line 131: “0.2% accuracy” – I believe, you are targeting 0.2% *error*…

Fig.1: “SSA” could be confused with single scattering albedo. Maybe SCA (scattering angle) sounds better?

Lines 50 & 92: “approximative” – please confirm the word is correct. In line 103, I see “approximate” twice which sounds better.

Line 78 and across the text: “polariSed” is correct, but “polariZed” is used more oftet, including this paper: 6 times for polariSed, 15 for polariZed (as Ctrl+F in Acrobat shows me)

Line 98: “…single plane parallel solution” – I was not able to understand what it means.

Lines 104-106: “The number of solar ... (2015)” – I would suggest to reformulate as “The number of solar zenith angles FOR the multiple scatter source function depends on …”

Line 143: “The comparison shown… SCIATRAN” – would it be better to reformulate it like this: “SCIATRAN uses fully spherical mode in the shown comparison”. In general, I have a feeling that the authors use passive voice too often.

Line 160: what is “three FULL dimensions”?

Line 191: “…wavy interface () or any surface spectral BRDF boundary condition” – I would suggest simplification: “…ocean and land.” It is clear that ocean has waves while land reflection is in general non-Lambertian and spectrally dependent.

Line 204: “While all four models listed above use Monte Carlo” – the sentence is excessive: it is clear from the Section title one line above.

Lines 254 & 310: please confirm “differing” is the right word

Line 261: what is “atmospheric state”? Please list all parameters.

Line 271: what it means to compare models “more directly”? The authors are likely talking about making input as close as possible.

Line 372: “It is possible….”- this should be either further explained or, better, dropped.

Line 435: last sentence in the section is very confusing. It sounds like SMART-G simulates only single scattering in spherical geometry. Please reformulate (or drop, because all necessary reference to SMART-G are in References).

Line 444-445: it is strange to “agree to a level better than the precision of calculation” (because difference in results cannot be treated accurately at that level). Would it be better to say something like “the codes agree within the expected accuracy?”

Line 455: the meaning (not the cause) of “discontinuity in radiance” in SMART-G remains unclear to me, but maybe I am just a bit tired by the end of the paper…

Line 466: “…upcoming ALTIUS mission which is linearly polarized” – sounds like the mission is polarized. I would cut the sentence at “…. ALTIUS mission.”

END-OF-DOCUMENT===================================================================================================

In the netCDF file, 'radiance' (and other arrays of the same shape) has the following diemsnion: (7, 3, 9, 3, 3, 11, 80, 4). Second, fourth and fifth dimension have the same size, 3. How can one make sure that data in those is not messed up? I.e., 2nd dimension is the test case, 4th is composition and 5th is albedo (as Panoply shows me).