A Semianalytic Monte Carlo Simulator for Spaceborne Oceanic Lidar: Framework and Preliminary Results

Spaceborne lidar (light detection and ranging) is a very promising tool for the optical properties of global atmosphere and ocean detection. Although some studies have shown spaceborne lidar’s potential in ocean application, there is no spaceborne lidar specifically designed for ocean studies at present. In order to investigate the detection mechanism of the spaceborne lidar and analyze its detection performance, a spaceborne oceanic lidar simulator is established based on the semianalytic Monte Carlo (MC) method. The basic principle, the main framework, and the preliminary results of the simulator are presented. The whole process of the laser emitting, transmitting, and receiving is executed by the simulator with specific atmosphere–ocean optical properties and lidar system parameters. It is the first spaceborne oceanic lidar simulator for both atmosphere and ocean. The abilities of this simulator to characterize the effect of multiple scattering on the lidar signals of different aerosols, clouds, and seawaters with different scattering phase functions are presented. Some of the results of this simulator are verified by the lidar equation. It is confirmed that the simulator is beneficial to study the principle of spaceborne oceanic lidar and it can help develop a high-precision retrieval algorithm for the inherent optical properties (IOPs) of seawater.


Introduction
The marine ecosystems are extremely complex and play an essential role in the biosphere. The phytoplankton is at the basis of the marine food web [1], and their annual net photosynthetic carbon fixation is nearly equivalent to that of all terrestrial plants [2,3]. The utilization and protection of the ocean is crucial to the sustainable development of human society. The ocean has been globally monitored by the spaceborne sensors for the past several decades.
Passive remote sensing of Ocean Color Radiometry (OCR) has provided a global view of the concentration of phytoplankton total suspended matter, colored dissolved organic matter (CDOM), and so on [4][5][6]. The performance of OCR spaceborne sensors has gradually improved since the In this paper, a spaceborne oceanic lidar simulator is built using a semianalytical Monte Carlo method, which can obtain the lidar returns from different atmosphere and ocean conditions. This is the first spaceborne oceanic lidar simulator that can be used to calculate the optical properties and the lidar signals of both atmosphere and ocean simultaneously. The simulator is helpful to explore the shapes of spaceborne lidar returns as a function of the seawater optical properties, to analyze the effects of multiple scattering, and to pave the way for the estimation of seawater inherent optical properties (IOPs), (i.e., the particulate back-scattering coefficient [4]). Furthermore, more application related to the ecology such as the study of polar phytoplankton biomass [34] and measurements of global ocean carbon sticks [35] can be done with accurate lidar data. The principle and method of the simulator is in Section 2. In Section 3, we use this simulator to obtain lidar returns from various ocean conditions, such as different types of scattering phase functions and the inhomogeneity of seawaters. Discussions on the preliminary results of the simulator are given in Section 4. A conclusion of this paper along with an outlook of possible improvements are outlined in the last section.

Principle of the Spaceborne Oceanic Lidar
The detection principle of the spaceborne oceanic lidar is given in Figure 1. The lidar system is mainly composed of the transmitting and the receiving systems. The transmitting system emits a pulsed laser into the atmosphere and ocean. Then, the laser pulse goes through absorption and scattering effects by aerosols and clouds in the atmosphere and by water molecules, CDOM, total suspended matter (TSM), and phytoplankton in ocean, as well as the refraction and reflection of the air-ocean interface. Finally, the detector in the receiving system gets the lidar return signals containing optical information of the atmosphere and ocean. The lidar return signal can be given by the lidar equation [19].
AT 2 a T 2 s (nH + z) 2 β π (z) exp(−2 z 0 k lidar (z )dz ) (1) where P(z) is the power received by the detector from depth z, P 0 is the transmitted power, K is the constant of the lidar system instrument, v is the speed of light in vacuum, τ is the pulse width, n is the refractive index of seawater, T s is the Fresnel transmittance of the air-ocean interface, A is the area of the detector, H is the altitude of the lidar above the sea surface, β π is the volume-scattering coefficient at a scattering angle of π rad, k lidar is the effective attenuation coefficient, and T a is the transmission through the atmosphere, which can be expressed as where α a is the extinction coefficient of the atmosphere. The lidar equation reveals that the power is related to the hardware parameters of lidar, the height of the lidar, as well as the optical properties of seawater. The simulator is designed to solve the lidar equation exactly. β π and k lidar are essential parameters for retrieving the optical properties of seawater. As is explained in reference [36], if the field of view (FOV) of the lidar receiver is large enough, the backscattering signal from the seawater is attenuated at a rate determined by the absorption coefficient. The radius of the receiving footprint on the sea surface from the spaceborne lidar system is usually tens of meters; that means that the FOV is large enough, so k lidar is close to the absorption coefficient a. Moreover, the scattering phase function described in Section 2.3 is also an important parameter. The particulate backscattering coefficient b bp can be obtained from the volume scattering function of particle β p (θ) by using Equation (3), which is a key biogeochemical parameter [35].

Simulation Method
The traditional MC method is very time-consuming. To improve simulation efficiency, several methods have been applied including a semianalytic MC algorithm, photon weight method, and multithreading operation method [19,26].
Based on the photon weight method, the semianalytic MC algorithm combines stochastic simulation and an analytical algorithm to calculate the expected proportion of the photon packet [26], which can return to the detector without further interaction. As shown in Figure 2 where p( ) is the scattering phase function, which represents the part of the photon packet scattered into an element of solid angle ∆Ω about direction [26]. The flow chart of the simulator is shown in Figure 3. To run the semianalytic MC simulation process with the simulator, lidar parameters (laser pulse energy, pulse width, height of the lidar, divergence angle of laser, etc.), environment parameters (absorption and scattering coefficients and

Simulation Method
The traditional MC method is very time-consuming. To improve simulation efficiency, several methods have been applied including a semianalytic MC algorithm, photon weight method, and multithreading operation method [19,26].
Based on the photon weight method, the semianalytic MC algorithm combines stochastic simulation and an analytical algorithm to calculate the expected proportion of the photon packet [26], which can return to the detector without further interaction. As shown in Figure 2, the photon packet moves from point O to point B in the direction → ε . The expected value of the collected photon packet can be described as where p(θ) is the scattering phase function, which represents the part of the photon packet scattered into an element of solid angle ∆Ω about direction θ [26].

Simulation Method
The traditional MC method is very time-consuming. To improve simulation efficiency, several methods have been applied including a semianalytic MC algorithm, photon weight method, and multithreading operation method [19,26].
Based on the photon weight method, the semianalytic MC algorithm combines stochastic simulation and an analytical algorithm to calculate the expected proportion of the photon packet [26], which can return to the detector without further interaction. As shown in Figure 2 where p( ) is the scattering phase function, which represents the part of the photon packet scattered into an element of solid angle ∆Ω about direction [26]. The flow chart of the simulator is shown in Figure 3. To run the semianalytic MC simulation process with the simulator, lidar parameters (laser pulse energy, pulse width, height of the lidar, divergence angle of laser, etc.), environment parameters (absorption and scattering coefficients and The flow chart of the simulator is shown in Figure 3. To run the semianalytic MC simulation process with the simulator, lidar parameters (laser pulse energy, pulse width, height of the lidar, divergence angle of laser, etc.), environment parameters (absorption and scattering coefficients and  In the semianalytical MC simulation, the Cartesian coordinates are used to describe the position and direction of the photon. The position of the photon is represented as (x, y, z), and the initial position is (0, 0, 0). The direction of the photon is represented by the cosine of the angle between the direction of the photon and the coordinate axes (ux, uy, uz). The schematic diagram of the photon scattering direction is shown in Figure 4. The origin direction is OA ⃗ , and the new direction after scattering is OB ⃗ . There are two angles that determine the scattering direction of the photon, including the scattering angle and the azimuth angle . The scattering angle represents the angle between the new direction after scattering and the previous original direction, ∈ (0, π). The azimuth angle represents the angle at which the photon is projected on the plane with respect to a reference direction, ∈ [0,2π). The step size s is defined as the integration of the attenuation coefficient c over the photon pathway, which is calculated by = −ln (1 − )/ , where is a random number uniformly distributed in [0,1]. The position is transformed from (x, y, z) to a new one (x', y', z') when the photon moved. The conversion equation is given by ' ; ' ; '  When the photon interacts with particles, the photon weight is updated by ′ = × , where W is the weight of the photon packet before collision and ω0 is the single scattering albedo (ω0 = b/c). The reflectance of the sea surface can be calculated according to the Fresnel formula if the sea surface In the semianalytical MC simulation, the Cartesian coordinates are used to describe the position and direction of the photon. The position of the photon is represented as (x, y, z), and the initial position is (0, 0, 0). The direction of the photon is represented by the cosine of the angle between the direction of the photon and the coordinate axes (u x , u y , u z ). The schematic diagram of the photon scattering direction is shown in Figure 4. The origin direction is → OA, and the new direction after scattering is → OB. There are two angles that determine the scattering direction of the photon, including the scattering angle ϕ and the azimuth angle θ. The scattering angle represents the angle between the new direction after scattering and the previous original direction, θ ∈ (0, π). The azimuth angle represents the angle at which the photon is projected on the plane with respect to a reference direction, ϕ ∈ [0, 2π). The step size s is defined as the integration of the attenuation coefficient c over the photon pathway, which is calculated by s = − ln(1 − ξ)/c, where ξ is a random number uniformly distributed in [0,1]. The position is transformed from (x, y, z) to a new one (x , y , z ) when the photon moved. The conversion equation is given by Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 16 phase functions of atmosphere and seawater, etc.) and output file parameters (vertical resolution of atmosphere and seawater, path of output file, etc.) need to be set. In the semianalytical MC simulation, the Cartesian coordinates are used to describe the position and direction of the photon. The position of the photon is represented as (x, y, z), and the initial position is (0, 0, 0). The direction of the photon is represented by the cosine of the angle between the direction of the photon and the coordinate axes (ux, uy, uz). The schematic diagram of the photon scattering direction is shown in Figure 4. The origin direction is OA ⃗ , and the new direction after scattering is OB ⃗ . There are two angles that determine the scattering direction of the photon, including the scattering angle and the azimuth angle . The scattering angle represents the angle between the new direction after scattering and the previous original direction, ∈ (0, π). The azimuth angle represents the angle at which the photon is projected on the plane with respect to a reference direction, ∈ [0,2π). The step size s is defined as the integration of the attenuation coefficient c over the photon pathway, which is calculated by = −ln (1 − )/ , where is a random number uniformly distributed in [0,1]. The position is transformed from (x, y, z) to a new one (x', y', z') when the photon moved. The conversion equation is given by When the photon interacts with particles, the photon weight is updated by ′ = × , where W is the weight of the photon packet before collision and ω0 is the single scattering albedo (ω0 = b/c). The reflectance of the sea surface can be calculated according to the Fresnel formula if the sea surface When the photon interacts with particles, the photon weight is updated by W = W × ω 0 , where W is the weight of the photon packet before collision and ω 0 is the single scattering albedo (ω 0 = b/c). The reflectance of the sea surface can be calculated according to the Fresnel formula if the sea surface Remote Sens. 2020, 12, 2820 6 of 16 is considered as flat. If the random number is less than the reflectance, the photon will be reflected; otherwise, it will be refracted into the seawater according to the refraction law.

Scattering Phase Function Models
The backscattered lidar signals from the atmosphere and ocean are closely related to their scattering characteristics. The scattering coefficient and the scattering phase function are important to represent the scattering characteristics of the media. The scattering phase function is the ratio of the volume scattering function (VSF) to the scattering coefficient, which specifies the angular dependence of the scattering without regard for its magnitude [37]. The phase functions of atmosphere and seawater included in this simulator are introduced below.

Scattering Phase Function of the Atmosphere
The phase function for atmospheric molecules is approximated by Rayleigh scattering, which can be given by [38] β with the polarization parameter p = 1. The phase function of the clouds and aerosols in the atmosphere can be calculated with Mie theory [39], as the radius distribution and optical properties are given. It is expressed as dr dr, and σ(r) is the scattering cross-section of the particle. i 1 and i 2 are the parallel and vertical components of scattered light intensity, respectively. dN(r)/dr is the distribution function of the particle radius.
The software package optical properties of aerosols and clouds (OPAC) [40] provides the microphysical and optical properties of 6 water clouds, 3 ice clouds, and 10 aerosol components, which are included in this simulator as typical and available cases. According to OPAC [40], Gamma and lognormal distributions are two main size distribution functions for particles in atmosphere. The equation of the Gamma distribution is [40] dN(r)/dr = Nar α exp(−Br γ ) (8) where B = α/γr γ mod , N is the total number density in particles per cubic centimeter, r mod is the mode radius in micrometers, α and γ are two constants that describe the slope of the size distribution, and a is the normalization constant. The equation of lognormal distribution is where b m and b p are the scattering coefficients of atmospheric molecules and particles, respectively, and b atm = b m + b p .

Scattering Phase Function of Seawater
The scattering phase function of seawater influences the shape of the lidar return signal. Limited by the measurement technology, the scattering phase function of seawater at small and large scattering angles has not been accurately measured. The typical measured phase function is the averaged-particle phase function derived from three measurements of VSF of the Bahama Islands, San Pedro Channel, and San Diego Harbor by Petzold (Petzold phase function for short) [41]. The Petzold phase function is given in form of discrete data and can be used by the lookup table method in MC simulation. At present, this simulator includes the averaged Petzold phase function, One-Term Henyey-Greenstein (OTHG) phase function [42], two-term Henyey-Greenstein (TTHG) phase function [43,44] and Fournier-Forand (FF) phase function [45]. Meanwhile, the user-defined phase function is also accepted, allowing the user to make any possible modification to the current ones or calculating with their own one. The OTHG phase function is a one-parameter function that is widely used in ocean optics because of its mathematical simplicity: [46] β OTHG (θ) = 1 4π where g is asymmetry factor, ranging from −1 to 1 [46]. The OTHG phase function gives a poor description of the particulate phase functions at large and small scattering angles [43]. Therefore, a weighted sum of the OTHG phase function, which is called the TTHG phase function, is proposed [43,44]: The parameter g 1 is given a value near 1, which makes the TTHG phase function increase more strongly at small θ than does the OTHG phase function. The g 2 parameter is given a negative value, which makes the TTHG increase with the scattering angle θ approaching π rad. The expression for parameters g 2 and α can be given as functions of g 1 [44]: The FF phase function is derived from a collection of Mie scattering with a Junge particle size distribution [45], which is based on the anomalous diffraction approximation [47] β FF (θ) = 1 where v = (3 − µ)/2, δ = 4/ 3(n − 1) 2 ] sin 2 (θ/2) , n is the true refractive index of particles, µ is the slope of hyperbolic distribution, and δ 180 • is the value of δ at θ = 180 • . The particle backscattering fraction B p represents the probability that a photon will be scattered through an angle ≥ 90 • , which is given by B p = b bp /b p . The subscript p represents particles, and b b is the backscatter coefficient, which is the integral of the VSF over the hemisphere of the backscattering directions. These mentioned phase functions are plotted in Figure 5 with B p = 0.0183, which is consistent with that of the Petzold phase function [46]. These seawater phase functions are all strongly peaked in the forward direction. In the backward scattering direction, both FF and OTHG are similar to Petzold. However, the TTHG phase function has a characteristic of high backscattering and low forward scattering compared to the other three phase functions according to Figure 5.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 16 function [43,44] and Fournier-Forand (FF) phase function [45]. Meanwhile, the user-defined phase function is also accepted, allowing the user to make any possible modification to the current ones or calculating with their own one. The OTHG phase function is a one-parameter function that is widely used in ocean optics because of its mathematical simplicity: [46] 2 OTHG 2 3 2 where g is asymmetry factor, ranging from −1 to 1 [46]. The OTHG phase function gives a poor description of the particulate phase functions at large and small scattering angles [43]. Therefore, a weighted sum of the OTHG phase function, which is called the TTHG phase function, is proposed [43,44]: The parameter 1 g is given a value near 1, which makes the TTHG phase function increase more strongly at small θ than does the OTHG phase function. The 2 g parameter is given a negative value, which makes the TTHG increase with the scattering angle θ approaching π rad. The expression for parameters 2 g and α can be given as functions of 1 g [44]: The FF phase function is derived from a collection of Mie scattering with a Junge particle size distribution [45], which is based on the anomalous diffraction approximation [  , n is the true refractive index of particles, is the slope of hyperbolic distribution, and ° is the value of at θ = 180°. The particle backscattering fraction Bp represents the probability that a photon will be scattered through an angle ≥ 90°, which is given by p bp p The subscript p represents particles, and bb is the backscatter coefficient, which is the integral of the VSF over the hemisphere of the backscattering directions. These mentioned phase functions are plotted in Figure 5 with Bp = 0.0183, which is consistent with that of the Petzold phase function [46]. These seawater phase functions are all strongly peaked in the forward direction. In the backward scattering direction, both FF and OTHG are similar to Petzold. However, the TTHG phase function has a characteristic of high backscattering and low forward scattering compared to the other three phase functions according to Figure 5.

Results
Some preliminary simulation results are presented here to demonstrate the performance of the simulator. The parameters of the elastic lidar system used in the simulations in this paper are given in Table 1. The air-ocean interface in this simulator is assumed to be flat, and the strong reflection from the surface is ignored, because that the lidar usually operates at an oblique incidence for ocean observation [48,49]. In order to save the simulation time and ensure a low level of standard deviation of the statistical results [19], all the simulations in this paper are executed with 10 8 photon packets.

Results of Atmosphere-Ocean Simulation
In order to give an overall view of the performance of this simulator, firstly, we simulate the spaceborne lidar return signal from the atmosphere and ocean. The atmosphere within 30 km (atmosphere above 30 km is clean air) is equally divided into 300 horizontal layers. The Rayleigh-scattering theory is adopted to calculate the optical properties of air molecules [50]. Two typical atmosphere cases are shown in Figure 6, where the blue dashed line is for aerosol and the red solid line is for cloud.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 16 increases with the decreasing altitude because of the increase in the attenuation coefficient of the atmosphere. The lidar signals of the aerosol and cloud layers are much stronger than that of the air molecule. The high attenuation of the cloud layer (with an optical depth of 1.2) causes the lidar signals of the molecule layer under the cloud layer and seawater to be weak. The signal of seawater is much stronger than that of atmosphere because of the high backscattering, and it attenuates quickly because of the high absorption and scattering effects.  Table 3. Optical properties of typical seawaters [51].

Verification of the Simulator's Accuracy
The normalized lidar signal of seawater obtained by the simulator and theoretical lidar equations are compared to verify the accuracy of the simulator. The atmosphere is set to vacuum, which means T a = 1 in Equation (1) to simplify the calculation. Three types of seawater (with IOPs shown in Table 3) with the Petzold phase function are used for simulation. The simulation results in Figure 7 (the solid lines) show that the signal attenuation slope is related to the optical properties of the seawater. The signal attenuation slope of the turbid water with a high attenuation coefficient is larger than that of the other two types of seawater. For the aerosol case, the atmosphere within 4 km above the ocean surface is assumed to contain maritime clean aerosols with a relative humidity of 50% called Maritime-50 [40]. The distribution of aerosol particles with height is described as [40] Remote Sens. 2020, 12, 2820

of 16
where h is the altitude above the ocean surface. Maritime-50 is composed of three types of aerosol components, including water-soluble, accumulate sea salt (acc.), and coarse sea salt (coa.). The properties of each aerosol component at wavelength of 532 nm are shown in Table 2. The absorption coefficient a and scattering coefficient b are normalized to a number density of 1 particle per cm 3 . The phase functions of these three aerosols are calculated by the Mie theory which is explained in Section 2.3 with the logarithmic distributions of the particle radius. The values of r mod and µ used in Equation (9) are shown in Table 2 according to the data given by OPAC [40]. The total values of the a, b, and phase functions of the mixed aerosol Maritime-50 can be calculated by where a i and b i are the absorption and scattering coefficients of aerosol component i, respectively. For the cloud case, the marine cloud with a number density of 50 cm −3 , droplet sizes of Gamma distribution, and the scattering coefficient of 0.0058 m −1 is assumed to be uniformly distributed in an altitude range of 600-800 m. The open ocean water with the IOPs given in Table 3 are used as an example [51]. The Petzold phase function is used in the simulations. The vertical distribution of the optical properties of atmosphere and seawater is shown in Figure 6a. The absorption of aerosol and cloud is weak, so the total attenuation mainly consists of scattering. The attenuation coefficient of the cloud layer is peaked. The lidar return signals from atmosphere and ocean obtained by the simulator are shown in Figure 6b. The vertical resolution of the signal in atmosphere is 30 m with reference to CALIOP [17], and the vertical resolution in seawater is set to 0.3 m. The spaceborne lidar return signal increases with the decreasing altitude because of the increase in the attenuation coefficient of the atmosphere. The lidar signals of the aerosol and cloud layers are much stronger than that of the air molecule. The high attenuation of the cloud layer (with an optical depth of 1.2) causes the lidar signals of the molecule layer under the cloud layer and seawater to be weak. The signal of seawater is much stronger than that of atmosphere because of the high backscattering, and it attenuates quickly because of the high absorption and scattering effects.

Verification of the Simulator's Accuracy
The normalized lidar signal of seawater obtained by the simulator and theoretical lidar equations are compared to verify the accuracy of the simulator. The atmosphere is set to vacuum, which means T a = 1 in Equation (1) to simplify the calculation. Three types of seawater (with IOPs shown in Table 3) with the Petzold phase function are used for simulation. The simulation results in Figure 7 (the solid lines) show that the signal attenuation slope is related to the optical properties of the seawater. The signal attenuation slope of the turbid water with a high attenuation coefficient is larger than that of the other two types of seawater.
To verify the accuracy of the simulator, we compare the normalized lidar return signal (dotted lines in Figure 7) with the simulated ones (solid lines in Figure 7) for three typical seawaters. It shows great consistence between the signals calculated by the lidar equation and the simulator, which verifies the accuracy of the simulator. In addition, the multiple scattering effect can be explained further. As the FOV is large enough, the interaction volume of the photon packet is within the FOV of the receiver. The attenuation of the lidar signal is almost caused by absorption without no scattering loss. The receiver can collect all the multiple scattered signals. Thus, the lidar effective attenuation coefficient is very close to the absorption coefficient. For further explanation, we have published a specialized paper about the relationship between the klidar value of the spaceborne oceanic lidar signal and the IOPs of case 1 waters under the fully multiple scattering condition [19].

Discussion
The purposes of this simulator are to obtain the lidar returns and analyze the law of signals under different conditions. Eventually, we will find the optimal lidar system parameters according to detection conditions and apply them to the actual system design. In this section, the lidar signal with different scattering phase functions and inhomogeneous seawaters are simulated and discussed.

Influence of Different Scattering Phase Functions
This simulator includes four types of oceanic scattering phase functions, as mentioned in Section 2.3. To compare the influence of different phase functions, the lidar return signals of homogeneous coastal ocean waters with four scattering phase functions ( Figure 5) are analyzed, and the results are presented in Figure 8. The lidar signal for water with the TTHG phase function is the strongest because of the higher backscattering probability of the TTHG phase function. According to the phase function in Figure 5, the intensity of the TTHG phase function for the scattering angle near 180° is much stronger than the other three, resulting in a higher lidar return signal. The lidar signal for water with the FF phase function is similar to that for the Petzold phase function. The lidar signal for water with the OTHG phase function near the air-ocean interface is lower, because the value of the HG phase function is the smallest at the scattering angle of 180°. Therefore, the backscattering characteristic of the phase function is the key factor in the lidar signal formation. The FOV in these simulations is 0.15 mrad, and the height of lidar is 700 km, so the radius of the footprint on the sea surface of the lidar receiving system is about 52.5 m. The effective attenuation coefficient k lidar is close to the absorption coefficient a of seawater, as explained in Section 2.1. Substituting the relationship k lidar ≈ a into Equation (1) and rewriting it, we have the normalized lidar signal P n (z) expressed by [19] P n (z) = P(z)(nH + z) 2 To verify the accuracy of the simulator, we compare the normalized lidar return signal (dotted lines in Figure 7) with the simulated ones (solid lines in Figure 7) for three typical seawaters. It shows great consistence between the signals calculated by the lidar equation and the simulator, which verifies the accuracy of the simulator. In addition, the multiple scattering effect can be explained further. As the FOV is large enough, the interaction volume of the photon packet is within the FOV of the receiver. The attenuation of the lidar signal is almost caused by absorption without no scattering loss. The receiver can collect all the multiple scattered signals. Thus, the lidar effective attenuation coefficient is very close to the absorption coefficient. For further explanation, we have published a specialized paper about the relationship between the k lidar value of the spaceborne oceanic lidar signal and the IOPs of case 1 waters under the fully multiple scattering condition [19].

Discussion
The purposes of this simulator are to obtain the lidar returns and analyze the law of signals under different conditions. Eventually, we will find the optimal lidar system parameters according to detection conditions and apply them to the actual system design. In this section, the lidar signal with different scattering phase functions and inhomogeneous seawaters are simulated and discussed.

Influence of Different Scattering Phase Functions
This simulator includes four types of oceanic scattering phase functions, as mentioned in Section 2.3. To compare the influence of different phase functions, the lidar return signals of homogeneous coastal ocean waters with four scattering phase functions ( Figure 5) are analyzed, and the results are presented Remote Sens. 2020, 12, 2820 11 of 16 in Figure 8. The lidar signal for water with the TTHG phase function is the strongest because of the higher backscattering probability of the TTHG phase function. According to the phase function in Figure 5, the intensity of the TTHG phase function for the scattering angle near 180 • is much stronger than the other three, resulting in a higher lidar return signal. The lidar signal for water with the FF phase function is similar to that for the Petzold phase function. The lidar signal for water with the OTHG phase function near the air-ocean interface is lower, because the value of the HG phase function is the smallest at the scattering angle of 180 • . Therefore, the backscattering characteristic of the phase function is the key factor in the lidar signal formation.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 16 Lidar signals for OTHG, FF, and Petzold phase functions tend to be equal, and the difference between the signal for the TTHG phase function and the other three tends to be small with the increasing depth. The reason for this tendency lies in the multiple scattering effect. The emitted laser beam maintains good collimation in the shallow water, and the backscattering angle is nearly 180°. As the depth increases, the laser beam gradually diverges due to the multiple scattering, resulting in a backscattering angle less than 180° [52]. With the increases of depth, multiple scattering shows stronger impacts on the lidar signal. The backward TTHG phase function decreases with the decrease of the angle, while the other three backward phase functions increase with the decrease of the angle, so that the gap of the signal decreases. More detailed information about the influence of the phase function on the lidar signal is shown in [52].

Lidar Signals from Inhomogeneous Seawaters
In reality, seawater is normally inhomogeneous, and the IOPs of case 1 water are closely related to the chlorophyll a concentration [Chla] . The vertical distribution of chlorophyll a concentration conforms to the Gaussian distribution given by [53]  Uitz et al. divided the world's stratified case 1 waters into 9 types (named S1-S9) based on the near-surface chlorophyll a concentration [53]. Three types (S1, S4, and S7 mentioned in reference [53]) of case 1 water are utilized here. The absorption coefficient can be calculated based on the bio-optical model of case 1 water [54] where ( ) w a λ is the absorption coefficient of pure sea water and * ( ) c a λ is a nondimensional statistically derived chlorophyll-specific absorption coefficient. The scattering coefficient is given by [54] ( ) where w ( ) b λ is the scattering coefficient of pure seawater, and its value is 0.022 m −1 at a wavelength of 532 nm. The vertical distributions of the inherent optical properties (IOPs) and chlorophyll a concentration for these three waters are shown in Figure 9. Lidar signals for OTHG, FF, and Petzold phase functions tend to be equal, and the difference between the signal for the TTHG phase function and the other three tends to be small with the increasing depth. The reason for this tendency lies in the multiple scattering effect. The emitted laser beam maintains good collimation in the shallow water, and the backscattering angle is nearly 180 • . As the depth increases, the laser beam gradually diverges due to the multiple scattering, resulting in a backscattering angle less than 180 • [52]. With the increases of depth, multiple scattering shows stronger impacts on the lidar signal. The backward TTHG phase function decreases with the decrease of the angle, while the other three backward phase functions increase with the decrease of the angle, so that the gap of the signal decreases. More detailed information about the influence of the phase function on the lidar signal is shown in [52].

Lidar Signals from Inhomogeneous Seawaters
In reality, seawater is normally inhomogeneous, and the IOPs of case 1 water are closely related to the chlorophyll a concentration [Chla]. The vertical distribution of chlorophyll a concentration conforms to the Gaussian distribution given by [53] where ζ is the dimensionless depth obtained by dividing the geometric depth z by the depth of the euphotic layer, z eu , ζ = z/z eu . C(ζ) = [Chla(ζ)]/Chla Zeu is the dimensionless concentration of chlorophyll a, and Chla Zeu is the average concentration of chlorophyll within the euphotic layer. C b is the background chlorophyll a concentration, which decreases linearly from the surface value [Chla] s with the slope s. C max is the maximum concentration, ζ max is the depth of the concentration, and maximum ∆ζ represents the width of the peak. Uitz et al. divided the world's stratified case 1 waters into 9 types (named S1-S9) based on the near-surface chlorophyll a concentration [53]. Three types (S1, S4, and S7 mentioned in reference [53]) of case 1 water are utilized here. The absorption coefficient can be calculated based on the bio-optical model of case 1 water [54] where a w (λ) is the absorption coefficient of pure sea water and a * c (λ) is a nondimensional statistically derived chlorophyll-specific absorption coefficient. The scattering coefficient is given by [54] b(λ) = b w (λ) + (550/λ) × 0.3 × [Chla] 0.62 (19) where b w (λ) is the scattering coefficient of pure seawater, and its value is 0.022 m −1 at a wavelength of 532 nm. The vertical distributions of the inherent optical properties (IOPs) and chlorophyll a concentration for these three waters are shown in Figure 9.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 16 The simulation results of these waters are shown in Figure 10. Different concentrations lead to different optical properties of seawater, which makes different slopes of the lidar signal. The attenuation of the lidar signal is stable, because the optical properties change slowly with depth. These results show that this simulator can be used for inhomogeneous seawater and even for the study of the subsurface scattering layers.

Conclusions and Outlook
This paper introduces the principle, process, and method of the spaceborne oceanic lidar signal simulator. The simulation model of the lidar signal from the atmosphere and ocean is established for different scattering phase functions and water types, and some preliminary results are presented. We compare the simulated lidar signal with the theoretical lidar equation to verify the accuracy of the simulator, and the results show great consistency between them. The simulator can get the spaceborne oceanic lidar signals for both homogeneous and inhomogeneous seawaters. The comparisons indicate that the scattering phase function has an obvious impact on the lidar signal, which is consistent with the conclusion about the phase function in reference [52]. The larger the The simulation results of these waters are shown in Figure 10. Different concentrations lead to different optical properties of seawater, which makes different slopes of the lidar signal. The attenuation of the lidar signal is stable, because the optical properties change slowly with depth. These results show that this simulator can be used for inhomogeneous seawater and even for the study of the subsurface scattering layers. The simulation results of these waters are shown in Figure 10. Different concentrations lead to different optical properties of seawater, which makes different slopes of the lidar signal. The attenuation of the lidar signal is stable, because the optical properties change slowly with depth. These results show that this simulator can be used for inhomogeneous seawater and even for the study of the subsurface scattering layers.

Conclusions and Outlook
This paper introduces the principle, process, and method of the spaceborne oceanic lidar signal simulator. The simulation model of the lidar signal from the atmosphere and ocean is established for different scattering phase functions and water types, and some preliminary results are presented. We compare the simulated lidar signal with the theoretical lidar equation to verify the accuracy of the simulator, and the results show great consistency between them. The simulator can get the spaceborne oceanic lidar signals for both homogeneous and inhomogeneous seawaters. The comparisons indicate that the scattering phase function has an obvious impact on the lidar signal, which is consistent with the conclusion about the phase function in reference [52]. The larger the

Conclusions and Outlook
This paper introduces the principle, process, and method of the spaceborne oceanic lidar signal simulator. The simulation model of the lidar signal from the atmosphere and ocean is established for different scattering phase functions and water types, and some preliminary results are presented. We compare the simulated lidar signal with the theoretical lidar equation to verify the accuracy of the simulator, and the results show great consistency between them. The simulator can get the spaceborne oceanic lidar signals for both homogeneous and inhomogeneous seawaters. The comparisons indicate that the scattering phase function has an obvious impact on the lidar signal, which is consistent with the conclusion about the phase function in reference [52]. The larger the attenuation coefficient, the higher the slope. The effective attenuation coefficient is commonly weaker in the atmosphere than that in the seawater. The relationship between the effective attenuation coefficient and the lidar system including FOV and lidar height is discussed carefully in reference [19]. The comprehensive simulation of the spaceborne oceanic lidar is a very complex project. Although we have achieved some preliminary results so far, this simulator can be further improved in terms of mechanism, performance, efficiency, and so on.
The semianalytical MC method, as the principle of the simulator, also can be used to simulate the polarized lidar and high-spectral resolution lidar (HSRL) signals (not shown here). Liu et al. [32] used the MC method to obtain the polarization lidar signal and compared it with the experimental data, as the polarization characteristics of the lidar signal include the IOPs information of seawater. Zhou et al. [55] developed a semianalytic Monte Carlo model to simulate the HSRL signals with multiple scattering. We are working on the integration of these simulation models for lidars with different mechanisms into the simulator.
In addition, the air-ocean interface is very complex because of the wind, bubble, and foam [13]. The complex sea surface influences the propagation properties of the laser beam and the form of the lidar signal. The air interface in the current simulator is assumed to be flat, and it can be improved to be a wind-roughed one.
Furthermore, when the mechanism and radiation transmission process are complete and complex, the efficiency of the MC method will be a serious limitation. At present, the simulator runs on the Central Processing Unit (CPU), and the simulation speed is slow. A Graphics Processing Unit (GPU) is especially suitable for processing the parallel tracking of multiple independent photons in MC simulation. Erik et al. [56] first proposed a parallel implementation based on the compute unified device architecture (CUDA) running on a GPU, and the speed is 1000 times faster than that of a CPU. Since then, GPU multithreaded parallel computing has gradually been widely used in biomedical optics and passive remote sensing [57,58]. We will apply the GPU acceleration technology into this simulator to improve the simulation efficiency.
The simulator will also be equipped with data calibration and retrieval abilities in the future. The limitation of the hardware parameters such as dynamic range should be considered in practical applications. We will also evaluate the actual parameters of the spaceborne oceanic lidar system, which is conducive to the design of the future spaceborne oceanic lidar system. In addition, through the simulation of various lidar signal models, a complete inversion algorithm is attempted to be established to obtain the IOPs of the seawater, which can be compared with the data of CALIOP to continuously improve the performance of the simulator [35,48], which is of great significance for global biomass and carbon stock estimation.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.