# Estimating Underwater Light Regime under Spatially Heterogeneous Sea Ice in the Arctic
## 2 Material and Methods
### 3D Monte Carlo Numerical Simulations of Radiative Transfer  
#### Theory and Geometry  
3D numerical Monte Carlo simulation is a convenient approach for modelling the light field under spatially heterogeneous sea surfaces [5; 12; 13; 14]. They are simple to understand and versatile, and incident light, IOPs and geometry can be easily changed. In this study, we used SimulO, a 3D Monte Carlo software program that simulates the propagation of light in optical instruments or in ocean waters [15]. Our objective was to simulate the propagation of sunlight underneath heterogeneous ice-covered ocean waters. Simulations were performed in an idealized ocean described by a cylinder of 120 m radius and 150 m depth (Figure 1). The water IOPs were selected to reflect pre-bloom conditions in the green-blue spectral region (\(a=b=0.05\) m\({}^{-1}\)). These typical averaged values were measured during the GreenEdge 2016 campaign using an in situ spectrophotometer (ac-s from Sea-Bird Scientific) and represent the contribution of both pure water and the water constituents. The scattering phase function was described by a Fournier-Forand analytic form with a 3% backscatter fraction [16; 17]. The inclusion of a 3D sea ice layer at the upper boundary of the ocean would require extensive computing power because of the high scattering properties of sea ice. Instead, sea ice was incorporated at the upper boundary of the ocean using a 2D light-emitting surface with a radius of 100 m. The angular distribution and magnitude of the light field emitted by the surface was chosen to mimic observed field data [18]. SimulO does not allow the use of arbitrary angular distribution for photon-emitting surfaces. To overcome this problem, two sources of photons were summed up in order to reproduce an observed under-ice light field (Figure 2). The first source was a regular Lambertian emitting surface while the second was a Lambertian emitting surface but restricted to an emission within 60 degrees of the zenith angle. A 5-m radius melt pond was set up at the center of the emitting surface (Figure 1). The melt pond had the same emitting angular distribution as the surrounding ice.  
Its intensity was four times higher than the surrounding ice, which corresponds to typical conditions found in the Arctic during summer [19].  
Figure 1: Spatial configuration used for the 3D Monte Carlo numerical simulations. (**A**) Surface view showing the percentage of the total area covered by the melt pond over the areas described by the black lines. For each of these areas, light profiles were averaged (see Figure 7). For visualization purpose, lines of the horizontal sampling distances from the centre of the melt pond have been plotted only at 5 m intervals. (**B**) 2D side view showing the 3D volume for which simulated data were extracted and how photon detectors were placed in the water column. Orange arrows indicate incident light sources.  
Given our interest in surface light profiles, 2D horizontal software detectors were placed vertically every 0.5 m, from 0.5 m up to a depth of 25 m. Detectors include 1 m\({}^{2}\) pixels measuring downward irradiance and upward radiance (5-degree half angle of acceptance). In order to avoid the effect of the boundary (i.e., absorption by the side of the cylinder used to simulate the water column), data outside a radius of 50 m were not used (see the green box in Figure 1). A total number of \(7.14\times 10^{10}\) photons were simulated to obtain a sufficient number of upwelling photons. The simulation took approximately 6000 h distributed over 2000 CPU cores. Because the geometry was symmetrical azimuthally, irradiance and radiance were averaged over the azimuth in order to increase the signal-to-noise ratio. Because of the low scattering coefficients used to reproduce in situ conditions observed during the sampling campaign, radiance profiles were noisy because a small number of upward photons could be captured. To address this issue, radiance profiles were smoothed using a Gaussian fit (Figure A2 included in the Appendix B).