Development and Validation of an End ‐ to ‐ End Simulator and Gas Concentration Retrieval Processor Applied to the MERLIN Lidar Mission †

In the context of MERLIN (MEthane Remote LIdar missioN), a French–German spatial lidar mission dedicated  to monitoring  the atmospheric methane content,  two software programs have been developed: LIDSIM (LIDar SIMulator) and PROLID (PROcessor LIDar). The objectives are to assess whether the instrument design meets the performance requirements and to study the sensitivity  of  this  performance  to  geophysical  parameters.  LIDSIM  is  an  end‐to‐end  mission simulator and PROLID is a retrieval processor that provides mole fractions of methane in dry air, averaged over an atmospheric column. These two tools are described in this paper. Results of the validation tests and the first full orbit simulations are reported. Merlin target performance does not seem to be reachable but breakthrough performance is reached.


MERLIN is a joint CNES-DLR (Centre National d'Etudes Spatiales and Deutsches
Zentrum für Luft-und Raumfahrt) satellite mission [1] designed to measure the atmospheric methane content at the global scale. Methane is a potent greenhouse gas largely emitted by human activities and responsible for more than 20% of the radiative forcing induced by well-mixed greenhouse gases [2]. To date, methane sources are estimated from surface observations supplemented by data from many passive instruments, using shortwave infrared (SWIR) absorption spectrometry, on board of satellites in sun-synchronous orbit. Table 1 provides a list of such instruments and the continuity of these observations is planned [3,4]. These passive satellite data have considerably improved our view of the space-time distribution of methane in the atmosphere [5]. However, these missions only provide data on sunlight reflecting areas, and their accuracy and precision are still far from the recommendation of the World Meteorological Organisation (WMO) for in situ surface observations of atmospheric methane, which is to be better than 2 ppb (<0.1%) [6]. The accuracy and precision of satellite data on atmospheric methane concentration must be increased to improve the identification of anthropogenic methane sources. The use of an active observation method, such as lidar, makes it possible to carry out measurements in all seasons, at all latitudes, by day and by night, and its implementation (small laser spot on the ground, differential approach) allows to limit systematic errors [7]. The MERLIN instrument is therefore an Integrated Path Differential Absorption (IPDA) lidar [14,15]. Two beams are pulsed at separate wavelengths (called On and Off) with significantly different molar methane absorptions ( and ), but close enough in frequency and time to minimise differences other than methane absorption in their various interactions with the atmosphere, the ground or the instrument's optics and detector. The principle of the IPDA measurement applied to the MERLIN mission is summarised in Figure 1 where orders of magnitude and values of various instrumental characteristics are given. The MERLIN instrument measures EOn and Eoff, the energies emitted per pulse [16], and, POn and Poff, those backscattered by the ground or by the cloud top, in order to provide the logarithm of the normalised energy ratio between the two frequencies. This quantity is the slant Differential Atmospheric Optical Depth (DAODslant) which represents the difference in the methane cross sections between the two frequencies, weighted by the methane concentration and integrated along the optical path. This relationship between the atmospheric methane concentration and the quantities measured by an IPDA lidar is expressed in Equation (1), where Nair is the number of moles of air per unit of volume and dl is the unit length along the path (1) Spectroscopic data and a priori knowledge of surface pressure, temperature and humidity profiles are used to estimate the dry air methane mole fraction (XCH4) from the DAOD. In addition, an average over many shots is necessary to reduce the individual pulse noise [17].  Table 2 summarises the user requirements for an integrated vertical methane column to improve methane emissions and sinks estimates [18]. The precision on XCH4 is demanding (~2% for random errors at threshold) and the accuracy challenging (~0.2% for systematic errors at threshold). It is also advantageous to give more weight to the lower troposphere, which the MERLIN measurement does by construction, and the specifications for industrial development take this into account. To meet these requirements, all sources of variability and uncertainty in the MERLIN measurements must be tracked. The user requirements expressed in Table 2 can be turned into constraints on the standard deviation of methane concentration measurements averaged over a given time period [18]. For a 7 s averaging time (corresponding to 50 km on the ground, due to the satellite velocity set by the orbit altitude) the limit for the squared deviation between the actual and estimated values of XCH4 is the squared sum of the random and systematic errors shown in Table 2. For longer averaging times, this squared error takes into account the averaging effect on the random noise. Conversely, the performance required for an individual shot (at 20 Hz) is reduced. This is because the random errors at this scale will be averaged over 140 shots to obtain the 50 km data. Figure 2 shows this transcription of Table 2. During the preparatory phase of the MERLIN mission, scheduled for launch in the second half of the decade, it is essential to verify by simulation whether the instrument design meets these challenging performances. This is the objective of the two software packages developed by the LMD (Laboratoire de Météorologie Dynamique): LIDSIM (an end-to-end mission simulator) and PROLID (a gas concentration retrieval processor). Instead of considering typical geophysical cases, the choice is made to be able to simulate data, in a systematic way, on one or several orbits, while having capacity to study specific cases. The software is developed in a modular way, in Fortran 90 with parallelization on the various target points treated, with the aim of being usable for instruments other than MERLIN, including the CHARM-F airborne demonstrator [19].
LIDSIM simulates the digitised lidar signal transmitted from the satellite to the processing centres. Its block diagram is shown in Figure 3a. The physical, chemical and radiative processes of the atmosphere and soil are first defined. (Section 2.1). The backscatter and extinction coefficients are evaluated for each atmospheric and surface level. The transmissions are then integrated at the relevant wavelengths and along the optical path both for the laser and solar fluxes (Section 2.2). The temporal distribution of light power is determined at the detector for atmospheric laser fluxes, solar flux and energy monitoring laser fluxes (Section 2.3). The various stochastic processes in the electronic detection chain are taken into account (Section 2.4). Random draws are performed to add simulated noise. Finally, the noisy power time series are convolved with the detection chain responses defined from the electronic characteristics (see Appendix A) to obtain voltage time series that are sampled and converted to digital counts (Section 2.5).
PROLID retrieves the atmospheric methane columns from MERLIN data simulated by LIDSIM, pending actual space-based measurements. Figure 3b shows its block diagram. PROLID is composed of two parts corresponding to the MERLIN data processing steps for level 1 and level 2 products [1]. An overview of the measurement principle and processing equations is given in Ehret et al. [1]. First of all, in PROLIDL1, the digital counts are pre-processed: subtracting an offset, filtering by a function to give a Gaussian shape to the signals, and determining a suitable filter window (Section 3.1). Then, for each peak in the pulse records, estimates are made for both its position in time to determine the mean scattering surface elevation (SSE) and its total energy to determine the mean slant DAOD (Section 3.2). Then, in PROLIDL2, the DAOD is corrected for the angle of incidence and the contributions of H2O and CO2 absorptions are subtracted. Finally, using a function calculated from auxiliary spectroscopic and meteorological data, the DAOD is converted for each shot into a vertically averaged methane content (Section 3.3). An average over many shots is used to reduce random noise at this stage. To deal with the various biases, the pulse energy and other quantities are averaged before estimating the average of XCH4 (Section 3.4). In yellow the inputs, in blue the processes, in green some variables, in red those related to noise, in grey the impulse responses of the electronic chain, and in purple the outputs.
The physical processes taken into account in LIDSIM and algorithms implemented in PROLID are described in MERLIN's Algorithm Theoretical Basis Document (ATBD) an internal project document [20]. Using LIDSIM and PROLID, the order of magnitude of the major system parameters can be determined (Section 4.1). Each step of these computations is carefully validated and the sensitivity of the results to various parameters is investigated (Section 4.2). The impact of noise is quantified using the standard atmosphere (Section 4.3). Finally, full MERLIN orbits are simulated to determine the precision and accuracy expected from the current instrument design (Section 4.4).

Ground and Atmosphere Description
The purpose of this section is to define the surface and atmosphere characteristics needed to compute the interactions between light and an Earth target.

Physical and Chemical Properties
Air is considered to be a mixture of ideal gases whose composition is constant for the major components, except for water vapour, although the minor components may vary (2) where P is the pressure, Nair is the number of moles of air per unit of volume, R is the ideal gas constant [21], T is the temperature, ρair is the air density, Mair/dry/H2O are the molar mass of the moist air, dry air and water vapour, respectively, qH2O is the specific humidity and Tv is the virtual temperature, which depends on both T and qH2O.
The molar masses (related to isotopic abundances) are in g/mol [22]: Mair = 28.965, MH2O = 18.015, MCO2 = 44.009 and MCH4 = 16.0425. For any G-gas of molar mass MG with NG moles per unit of volume, XG the mole fraction referring to dry air (expressed in %, ppm or ppb) is associated to the dry mass ratio rG referring to dry air and to the total mass ratio qG referring to moist air and 1 For sensitivity studies, the temperature profile is that of the standard atmosphere [23] defined by the pressure and temperature at sea level, and by the gradients, relative to the geopotential, of the virtual temperature, which are constant in the geopotential bands. The humidity profile decreases exponentially with the geopotential to a minimum value. Analytical functions are used, so no vertical interpolation is required to calculate the different quantities at the various altitudes defined from the arbitrarily fixed ground elevation. An alternative is to use temperature and humidity values at a set of pressures, taken from the Thermodynamic Initial Guess Retrieval (TIGR) radiosonde database [24] and to perform vertical interpolations at the desired altitudes. As usual, the interpolations used are linear in pressure for humidity and for the virtual temperature deviation from the standard atmosphere, and linear in pressure logarithm for geopotential deviation from the standard atmosphere. Figure 4 shows the temperature and humidity profiles for the standard atmosphere and for five TIGR profiles.  Regarding aerosols and clouds, one of the five vertical aerosol distributions and one of the seven cloud types from the ESA-defined reference model [25] can be selected. For the orbit simulations, the location of the laser footprints can be determined from the orbital characteristics using the IXION software [26] or taken from any database. Data on meteorological parameters [27], aerosols [28,29] and gas concentrations [30] are taken from operational analyses provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) (Examples of these data sets are shown in Figure 5a-c). Following the habits of the meteorologists, they are interpolated linearly in time and bi-linearly on the horizontal, and then vertically at defined altitudes from the ground elevation. The latter is derived from a Digital Elevation Model (DEM): the EarthEnv DEM [31] at a 90 m resolution ( Figure 5d shows this DEM over southern Italy) supplemented by data from the GTOPO30 database [32] at a 1 km resolution for areas where EarthEnv data are missing. The vertical interpolation consists of first determining the pressure of the target elevation by linear interpolation in geopotential of the logarithm of the pressure, and then interpolating T, q and xG to this pressure by linear interpolation in pressure. Extrapolations below ground are carried out by keeping a constant gradient from the penultimate layer above ground: zero gradient for specific humidity and any molar mass fraction, and standard gradient for temperature (−6.5 K/km). To preserve boundary layer profiles close to the ground, it is possible to use a more elaborate algorithm, developed at Météo-France (but unpublished), which consists of performing a weighted average between the profile found previously and a profile constructed from the potential temperature of the ground, thus preserving for each layer the gradient of potential temperature on the one hand and relative humidity on the other. A known disadvantage of this approach is the change in the geopotential even at higher altitudes, but this is of little importance here, because the geopotential is only involved in terms of the difference between two levels.

Radiative Properties
The laser emissions of the two MERLIN pulses are nominally centred at λOn = 1645.5516 nm and λOff = 1645.8460 nm. The detection optics contain a filter to limit the incoming solar flux to a Δλfilter interval of 2 nm. The spectral range is defined in wavelength as [1644.7 nm; 1646.7 nm] or in wavenumber as [6072.75 cm −1 ; 6080.14 cm −1 ]. The spectroscopic parameters of the gases of interest (CH4 but also CO2 and H2O) are taken from the GEISA database [33], with recent specific improvements for methane [34,35] and for water vapour [36]. For each G-gas, the vertical distributions of the molar cross sections , , are calculated from the thermodynamic description of the atmosphere using the 4A radiative transfer software [37,38] which has been modified to provide absorption cross sections instead of optical thicknesses. The Hartmann-Tran profile [39], recommended by the International Union of Pure and Applied Chemistry, is used for methane lines air-broadening and -shifting and combined with first-order line-mixing approximation. The 4A computations are performed at the maximum spectral resolution of 5 10 -4 cm -1 (about 15 MHz). Figure 6a shows the absorption cross sections on the spectral region of interest for CH4, CO2 and H2O, and Figure 6b shows the absorption coefficients for typical abundances of these gases. The atmospheric vertical distributions of extinction and backscatter coefficients for air molecules (αM and βM), aerosol particles (αA and βA) and cloud droplets (αC and βC) are calculated in accordance with the ESA recommendations [25], so that the following relations are imposed For molecules, the extinction resulting from scattering is added to the molecular absorption computed above. Figure 7 shows the vertical distribution of backscatter and extinction coefficients estimated at a mean frequency for molecules, aerosols and clouds. It is assumed that these coefficients are not wavelength dependent over a few nanometres. For clouds, with Hbot and Htop representing the altitudes of the base and top of the cloud layer, respectively, multiple scattering is taken into account for the downwards solar flux. Scattering in the cloud layer is replaced by scattering at the top of the cloud assimilated to a Lambert surface with a hemispherical reflectance Rc. Rc is fixed using τ the optical thickness of the cloud (computed from the extinction coefficient) and ga = 0.85 the asymmetry factor which characterises the phase function of cloud scattering [41,42] 1 The ground is assumed to be a set of scattering surfaces distributed vertically according to with 1. The reflectance of the scattering surfaces varies with incidence and is prescribed as constant for the wavelengths of interest: Rnad for the laser beam (near nadir incidence) and Rhia for sunlight (high incidence angle). For sensibility studies, a typical value of 0.1 is used and is considered as a set of Dirac or Gaussian distributions. For the simulated orbits, the reflectances (see Figure 8) are derived from MODIS [43][44][45] with auxiliary data from SeaWifs and QuikScat for sea reflectances [46] and is estimated, based on the DEM variability between adjacent data, as a weighted sum of nine Dirac distributions.
(a) (b) Figure 8. Global maps of zenith reflectance around 1.6 μm on a monthly basis (here June). (a) for a nadir incidence (hot spot configuration), the map provided by the company Estellus under a CNES contract, over the continents with a lat/lon resolution of 1/120°. These data are used for the laser fluxes even though, with an incidence angle of −1°, there is no specular reflectivity in the MERLIN context. (b) for a solar incidence angle of 65°, the map provided by the company Noveltis under an ESA contract [46], with resolution of 1/10°. These data are used for sunlight even if the incidence angle (given the choice of orbit) is larger as the database has too many artefacts at higher angles of incidence.

From Radiative Properties to Transmission and Backscatter Coefficients
The light transmission factor between the satellite and an atmospheric layer, as well as the scattering coefficient of light towards the satellite from an atmospheric or ground level, are determined according to the observation geometry and the associated Doppler effects. The calculations are performed in spherical geometry with the assumption of horizontal homogeneity, and considering the variation of the speed of light with air density. Multiple scattering is only taken into account for the downward solar flux.

Observation Geometry and Vertical Coordinates
The lidar pulses are emitted in the ⃗ direction towards the Earth, with an angle θsat to the vertical and an angle to the North direction at the satellite location. r indicates the distance between the satellite and the backscatter target along the beam propagation direction. Figure 9a illustrates the observation geometry as modelled with the location of the different variables.
For the geometric calculations, the geodetic system used is WGS84 [47,48]. Its reference ellipsoid (defined by REq = 6378.137 km its terrestrial radius at the equator and e = 0.081819191 its flatness) is locally approximated by a tangent sphere of radius Rcurv as a function of latitude lat and azimuth direction Due to the density inhomogeneities of the Earth, the geoid differs from the reference ellipsoid. H the altitude (relative to the geoid) therefore differs from dell the distance to the ellipsoid by N the height of the geoid above the ellipsoid. N is fixed arbitrarily for the sensitivity studies and taken from the Earth Gravitational Model 2008 (EGM08) dataset [49] (shown in Figure 9b) for the orbit simulations A spherical atmosphere is used to determine dr, the contribution to the optical path of a vertical layer of thickness dH. With μ(H) the angle at altitude H between the local vertical and the direction of beam propagation, Rcurv the local curvature radius of the ellipsoid, N the height of the geoid above the ellipsoid, Hsat = dsat − N the altitude of the satellite above the geoid and θsat = μ(Hsat) with (8) and (9) The gravity above the geoid is assumed to be the normal gravity above the reference ellipsoid. It varies with latitude and along the vertical. This last dependence is expressed as a Taylor series with altitude. In LIDSIM, the gravity varies as the inverse of H 2 but the parameter Rgrav allows the Taylor series to be adjusted. Then, the gravitational field is defined as and (11) with H the altitude, gEq = 9.780327 m/s 2 the gravity on the geoid at the equator, g1 = 5.3024 10 −3 and g2 = 5.8 10 −6 two constants to approximate Somigliana's formula, rg1 = 1.0068 and rg2 = 6.7056 10 −3 two other constants to define Rgrav, and REq = 6378.137 km the Earth's radius at the equator. Figure 10a shows the latitude dependence of the different radii used. As the satellite is intended to be manoeuvred so that ⃗ forms a constant angle tau = −1° with the direction of the centre of the Earth, θsat varies with latitude (see Figure 10b). To use the meteorological data provided on vertical levels generally defined by their pressure, Z, the geopotential in J/kg, and H, the altitude in m, must be linked The time trt required for light to travel the distance r and return is determined using the actual light speed calculated as a function of the air index nair(H) where c is the light speed in vacuum [21] and nair(H) is the refractive index along the path. nair is calculated according to Ciddor [50] from the number of moles per unit volume of dry air and water vapour (Ndry and NH2O) with two constants ciddry = 6.49 10 −6 and cidH2O = 5.57 10 −6 calculated for a wavelength set to λref = 1645.699 nm and an abundance of CO2 fixed at XCO2 ref 1 (14) Figure 11 shows the variation of the speed of light with pressure and the corresponding error in the distance to the satellite if the speed of light is assumed constant.

Spectral Distribution of the Lidar Emission and Doppler Effects
In the satellite frame of reference, the lidar emission is assumed to follow a normal frequency distribution with a central frequency νo (alternatively νOn and νOff for a pair of pulses) and a standard deviation corresponding to a full width at half maximum (FWHM) of 60 MHz Spectral impurities can be introduced by adding similar Gaussian energy distributions on either side of the emission frequency, at significantly lower intensities and at frequencies shifted by one or more free spectral intervals [51].
The movement of the ground or of any air parcels relative to the satellite generates νDop, a Doppler shift between νsat the wavelength in the reference frame of the satellite and νtar the wavelength in the reference frame of the target. This Doppler shift is proportional to ⃗ the relative speed where ⃗ is the unit vector pointing to the laser beam emission direction and c is the speed of light.
As the wind speed's contribution to the Doppler shift is negligible for near nadir observations (estimated at less than 1 MHz), the generated Doppler shift νDop is constant on the vertical. It is determined by the scalar product ⃗ . ⃗ where ⃗ is the relative speed of the satellite with respect to the ground ( ⃗ = 7.6 km/s). Moreover, as the tilt angle tau is chosen to minimise this Doppler effect (tau = −1°), νDop is negligible (less than 5 MHz or 3.3 10 −5 cm −1 ). In addition, the velocity distribution of the molecules and particles creates a spectral broadening of the light beam scattered by them (illustrated in Figure 12). Assuming a Maxwell-Boltzmann distribution for the velocities of the scatterers (of molar mass Msca), the probability for velocity Vz on the incident direction is with (17) The probability distributions of forward and backward frequencies are linked by the conditional probability GRay of observing νb on the backward path when νf is observed on the forward path. Assuming elastic interactions, GRay depends only on the shift between νb and νf which is directly related to Vz and νref, a medium value used for all the frequencies as their dispersion is really small 2 and As particles are heavier than molecules, their Brownian motion corresponds to lower individual velocities. Therefore, for the time being, the broadening associated with particles scattering is neglected with respect to the spectral width of the pulse. Due to the large number of photons in the considered flows, their interactions with molecules or particles are classically treated in a deterministic way.

Vertical Distribution of Transmission and Backscatter Coefficients
TL(ν,H) refers to the atmospheric transmission at frequency ν, integrated along the direction of the beam propagation between the satellite and the observed scattering target at altitude H. It is calculated from the vertical distribution of the coefficients αM/A/C(H) representing the extinction due to air molecules, aerosols and clouds, respectively, and from the cross sections σG(ν, H) representing the absorption of the minor gases involved where NG(H) is the number of moles of G-gas per unit of volume at altitude H and Htoa ≅ 40 km the upper limit of the atmosphere above which light-matter interactions are negligible.

TL(ν,H)
is calculated for the two bands centred on the wavelengths emitted by the lidar and whose width is defined by the Doppler broadening, with a full spectral accuracy of 5 10 −4 cm −1 (see Figure 12).
A similar computation, but with a spectral accuracy reduced to 5 10 −2 cm −1 to save computing time, is done for TS(ν,H) the atmospheric transmissions of the solar flux on its own pathway defined by θsun the zenith solar angle given at the ellipsoid level where the angle μsun(H) is calculated from the angle θsun and the radius (22) is determined using Equation (6), from , the azimuthal angle of the sun, instead of that of the laser beam.
Furthermore, , , (in m −1 sr −1 ) the backscatter coefficient for the laser fluxes is given, with referring to the Dirac function for the wavenumbers, by , , , , (in m −1 sr −1 ) represents the part of the solar flux scattered towards the satellite from the ground or the cloud top and is defined by where refers to the Dirac function for the altitude.

Radiative Fluxes on the Detector
In this section, the quantities TL(ν,H), TS(ν,H), , , and , , determined in the previous section are used together with the instrument properties to specify not only , the backscattered flux powers, but also , the solar flux power, at both the On and Off frequencies. Simultaneously, , the monitoring fluxes powers, are estimated at both the On and Off frequencies.
The main optical properties of the instrument are EL = 9.5 mJ (alternatively EOn and EOff for a pair of pulses) the total energy of the emitted laser pulses, Eo(ν) the optical emission efficiency factor, A the effective area of the collecting mirror and Do(ν) the optical reception efficiency factor including the bandpass filter used to limit the solar flux reaching the detector. This filter is assumed to be an idealised bandpass filter, with a rectangular curve centred on λref =1645.699 nm with Δλfilter = 2 nm full width. Such a filter does not bias the estimation of DAOD. The effects due to differences in transmissions through the filter for the On and Off pulses are outside the scope of this work but are being studied by Airbus. If necessary, these effects will be integrated into LIDSIM.

From Emitted Lidar Pulses to Backscattered Fluxes on the Detector
The radiative transfer laws provide /m (unit m −1 ) the fraction of energy that returns from altitude H per unit of length after a round trip through the atmosphere in the integrals, νf denotes the frequencies for the forward path and νb for the backward path, as seen by the Earth targets. Do(ν) and Eo(ν) are assumed to be constant over the range of the transmitted frequencies, and are set to 0.77 and 0.952, respectively. /m is then computed using Gat2 (26) where the auxiliary variables Gat1 and Gat2 average the spectral dependencies of the transmissions at the various altitudes Gat1 , Gat2 with 3√2 the parameter fixing the truncation of the Gaussian distributions.
/m /s is estimated on a set of vertical layers where the radiative properties have been previously computed, but the corresponding set of times trt(H) is an irregular basis of time. Therefore, to evaluate the convolution product, /s must be resampled. A resampling procedure is implemented to conserve the energy (the number of photons) per time interval. (See Appendix B).

Solar Flux
Psun(υ) the solar spectrum at the top of the atmosphere (shown in Figure 13) is taken according to Toon [52] and converted in W/m 2 /cm −1 with Ωsun = 68.36 μsr the mean solar solid angle. There is no modelled temporal variation in solar flux on the scale of a pair of pulses, although a displacement of a few metres may induce a significant fluctuation in the presence of clouds. Seasonal and diurnal modulations are only taken into account by the sun zenith angle. The mean Doppler shift due to the motion of the Earth on its elliptical orbit and its rotation on itself is less than 0.5 cm −1 and is not taken into account because the seasonal fluctuations of the solar background are out of the scope of this study. As already mentioned, Doppler broadening is also negligible as the solar spectrum does not show any variability at the scale of a few hundredths of cm −1 .
the solar power on the detector is then constant over time (on a shot scale) and can be obtained as an integral over the scattering altitude 2 GatS (31) where δθ is the angle determining the field of view (the part of the Earth from which the backscattered photons reach the detector through the receiving optics of the instrument), μsun(H) is calculated, using Equation (22) from the angle θsun. The auxiliary variable GatS averaging the spectral dependencies of the transmissions is , ,

Energy Monitoring Fluxes
To handle fluctuations in the energy emitted by the laser, the instrument has a calibration unit that records the energy emitted through the same detector and electronic chain as the signal backscattered from the Earth targets [15]. The equation for , the energy monitoring flux, is simply /s (33) where /s = 0.31 10 −12 is the fraction of emitted energy injected on the detector by the calibration unit and td0 is the delay (calibrated before launch) between the laser emission time and the corresponding detector illumination time. To avoid problems of nonlinearity with the detector sensitivity, the energy of the monitoring fluxes must be of the same magnitude as the atmospheric fluxes, so integrating spheres are needed to reduce the energy taken out during the beam emission [18].

Electronic Chain and Noises
The electronics chain (shown in Figure 14) consists mainly of a detector (APD), a trans-impedance amplifier (TIA) and an analogue-to-digital converter (ADC). The detector is an avalanche photodiode that converts photons into electrons and amplifies their number. The TIA transforms the small current associated with the movement of the electrons into a voltage that can be digitised by the ADC. An anti-aliasing filter (AAF) and other electronics to control the offset (OC) and signal gain (PGA) are located between the TIA and the ADC. Five fluxes have to be considered on the APD: the On and Off laser beams after their return from the atmosphere, the On and Off laser beams for calibration and the solar flux. Due to the large number of photons, the interactions with the atmosphere and the ground can be considered as deterministic, but in the electronic chain, many stochastic processes have to be taken into account to simulate the mean signals and their variability.

Speckle Effects
A first random process is related to interference on the detector resulting from the differences in optical path between different points in the area illuminated by the laser.
Mtav the speckle number is determined from Sc the characteristic area, Tc the characteristic time, A = 3850.51 cm 2 the effective surface of the collecting mirror, δT the time interval considered and P the polarisation index of the flux on the detector The laser fluxes are polarised unlike the solar flux. can be considered infinite for a pulsed laser, because for each pulse there is only one temporal speckle due to the full correlation of the transmitted signal over the pulse duration. Moreover, even if the interference pattern on the detector evolves due to the motion of the satellite relative to the ground, this evolution is negligible over the duration of a laser pulse [54,55]. The speckle noise during a laser pulse is therefore fully correlated. Table 3 gives the values of the parameters of Equation (35), calculated for the MERLIN instrument characteristics, and and the resulting speckle number for the laser and solar fluxes, respectively. 1 ⁄ can be considered as zero due to the large area (field of view of the receiver) from which the solar flux originates and the relatively wide spectral filter. The monitoring fluxes are depolarised by the integrating spheres, allowing the energy level to be reduced. But the complete analysis of the interference is complex and very dependent on the instrument design (integrating spheres and optical fibres). A value of = 1850 is prescribed in LIDSIM following a dedicated study by Airbus. In addition, a specific mechanism is installed to decouple the speckle patterns between each shot [15]. However, decorrelation can only be achieved at the time scale of a shot and not at the time scale of a sample [56]. Thus, the speckle pattern for the monitoring streams is fixed for a full shot but changes from one shot to another. So, the speckle noise on a monitoring pulse is fully correlated.
The effect of speckle noise is to alter the energy estimate of the pulses, each single-shot DAOD can be considered biased. But if there is no correlation on the changes between successive shots, then the final effect on the average energies and an averaged DAOD will be similar to random noise.

Shot Noise
While the speckle noise is due to fluctuations in the optical energy at the detector, arising from interferences reflecting the wave nature of light, the shot noise is due to fluctuations in the number of photons arriving in a given time interval, reflecting the corpuscular nature of light.
The arrival of photons at the detector follows a Poisson law [53]. The mean and variance of the photon number Nph within any δt-time interval are derived from Mandel's formula [57], where h is Planck's constant and υOn/Off is the frequency

Photoelectric Effect and Avalanche Process in the APD
In the APD, the photoelectric effect is modelled as a Bernoulli dilution characterised by a quantum efficiency η. The avalanche process is characterised by an avalanche factor M = 10 and an excess noise factor F [58]. η and F can be computed from Rdet = 0.88 A/W the detector responsivity and κ = 40 the effective ratio of electron and hole ionization coefficients [59] ℎ and 1 2 1 where h is Planck's constant [21], e is the charge of an electron [21] and ν is the frequency. Each photon reaching the detector has the probability pi of emitting a primary electron with <pi> = η and Var(pi) = η (1 − η). Each primary electron is transformed into mi electrons during the avalanche process with <mi> = M and Var(mi) = (F − 1) M 2 . Using Burgess' variance theorem [60,61] for k independent random variables with the same distribution The mean and variance of the number of primary electrons produced by the APD are and those of the total electron number are

Electronics Noises
Based on the previous paragraph, the APD receives a photon flux P(t) = δE(t)/δt of mean power <P>(t) and generates an electron flux corresponding to the intensity ip(t) = δq(t)/δt, with mean <ip>(t) and variance Var(ip)(t) where e is the charge of the electron and δ(t) is the Dirac distribution in time.
The TIA and AAF transform this current into a voltage v1 but add various electronic noises to the signal. Appendix A defines in(t) and ucn(t) as white noises characterised by în and ^ as their Amplitude Spectral Densities (ASD) and amplified by Ri (like the signal) and by Ru, respectively. ⊗ ⊗ ⊗ , the voltage after the AAF thus has the following mean and variance

From Radiative Fluxes to Digitised Signal
According to the analysis described above, all noises are simulated as Gaussian noises with the appropriate variances and no covariance between samples, except for the speckle noise for which the interference patterns are held constant for a full pulse time by performing a single draw for all discretisation times of that pulse.
Random draws are performed and the electric current time series for the five processed streams are given on a discrete regular time basis {tj} with tj+1 − tj = δt by With Equations (30), (31), (33), (43) and (44) and their notations, the digital counts expression is The first terms provide the mean value, with the offset and solar terms causing a bias towards the expected mean value. The other terms correspond to random Gaussian noise. To complete the simulation of the MERLIN data, these digital counts are averaged for the times corresponding to high altitude measurements in order to emulate the reduction in the volume of data to be transmitted from the satellite to the ground that will be implemented on board. Figure 15 shows an example of a simulated signal with and without noise. For further analysis, the following expressions give the mean and covariance of a sum of digital counts over a W-window

PROLID: Processing Lidar Data
The processing of the lidar signals follows the MERLIN ATBD [20] but may differ to some extent. PROLID is a research prototype to help define the algorithms that will be used in the official MERLIN ground segment. It is divided into two parts. Firstly, each peak is dated to obtain an estimate of the scattering surface elevation (SSE). The average energy of each peak and its variability are also estimated to derive the differential absorption optical depth (DAOD) along the optical path. Secondly, using ancillary meteorological data, a weighting function, then XCH4, are computed for each shot. However, due to the noise level, an average over several shots is required to obtain an estimate whose accuracy meets the mission objectives.

Preliminary Treatments
As the sampling of the signal varies in time with its altitude range, the signal is first reconstructed with a uniform time base by simply assigning the mean value for each corresponding data item.
For each shot, Rout, a region with NR values without signal influence, is defined and the average signal value on this region is computed ⁄ can be considered as an offset and expressed in terms of the LIDSIM parameters When the instrument is in flight, the electronic characteristics, which may drift with time, will be unknown. However, it will be possible, using Equation (51), to estimate UD the normalised electronic impulse response convolved with the time emission. To avoid broadening due to the vertical distribution of ground reflectance, only data from the energy calibration path will be used. The digital counts of the individual shots used to determine UD will be averaged after removing the amplitude offset and aligning them to an appropriate time base (the signals will be shifted to fit each peak to its maximum) To improve the estimation algorithms that follow, the digital counts are filtered by Ke a function defined such that its convolution with UD gives G a Gaussian centred function with standard deviation σG (σG = 6.38 following a recommendation from the ATBD but this value will be adjusted in future work) To obtain a smooth Ke function, the low values of G(i) and Ke(i) are forced to zero. After various tests, the limits defining the "low values" are empirically set to a fraction 10 −9 and 10 −6 of the maxima of G and K, respectively. For the processing, ⁄ and ⁄ are computed as As long as the photon signals on the detector can be considered as quasi-Dirac pulses, the signals received at the ground convolved with Ke have a quasi-Gaussian shape. Despite the different noises, it will then be easier to estimate the time location of the signal, its energy and also its signal-to-noise ratio (SNR).

Estimations
From the NCFC(i) digital counts, a peak (or two if there is a transparent cloud) is identified corresponding to the monitoring pulse or return signal record, for each of the two wavelengths On and Off. Each peak corresponds to a photon pulse reaching the detector convolved by the response function. Each peak is associated with the index of its maximum ∈ (55) To account for the temporal dispersion introduced by the vertical distribution of the target reflectance, the window size used for the estimates is adjusted on the data. The search for the half-value width of the signal is performed on the νOff -frequency signal (as it is more energetic). A search interval is set and a first window is constructed by searching, from the limits of the max interval and moving towards the centre, for the first indices left and right for which the signal is greater than the half of its maximum . This window is then expanded on either side by half the length of the noise correlation. This procedure is illustrated in Figure 16. The same window size

Elevation of the Scattering Target
⁄ the recording time of a peak for backscattered or monitoring signals is estimated as centroids on the laser pulse at the νOff-frequency The difference between the peak recording times of the Earth and calibration signals is converted into dist, the distance between satellite and target. This is computed by taking into account tdo the delay between a laser emission and its recording (calibrated before launch and arbitrarily fixed at 1760 ns in PROLID for the time being) and with c the average speed of light on the path estimated using Href the height scale of the atmosphere where P00 = 1013.25 hPa, T00 = 288.15 K, g00 = 9.80665 m s −2 are the standard values of pressure, temperature and gravity on Earth, R is the ideal gas constant, Mair is the molar mass of dry air and ciddry = 6.49 10 −6 is from Ciddor [50] 2 The scattering surface elevation (SSE) is then determined by a geometric computation using a sphere of constant radius Rterre = 6371.226 km instead of the curvature of the geoid (but not a flat Earth, as would induce a bias of 6 m) The target is identified as the ground or a cloud by comparing the SSE and the DEM (not yet studied in detail, especially in the case of clouds close to the ground). The digital counts can be divided into a part of photonic origin xp that is unsteady but results from a convolution by a Gaussian function, and a statistically stationary dark component xd.

Energy Peaks and Their Variability
By using the following property of Gaussian distributions can be estimated using orbit sections without solar illumination, and K can be estimated [20] using sections with different constant illuminations to establish a regression between the variance of the signal over the section and the average of the Gvar-filtered signal over the same section. At present, these quantities are computed in PROLID directly from LIDSIM parameters

DAOD along the Optical Path
DAODslant on the optical path is estimated from ⁄ ⁄ the total energy peaks which are assumed to have a Gaussian distribution but for the estimation of its mean, due to the non-linear transformation, a statistical bias coming from the dispersion of the peak energy distribution has to be taken into account [17] (70) < > 1 2

XCH4 Retrieval
The actual vertical average of XCH4, i.e., the XCH4 column, is But, the IPDA method actually averages the observed quantity over the vertical in a specific way. The weighting function WF(p) is used to define a reference observable value corresponding to what the instrument actually observes : (72) The weighting function is estimated using spectroscopic and meteorological auxiliary data , , where g(lat, H(p)) is the vertical acceleration of gravity for latitude lat and altitude H corresponding to pressure p, Mdry is the molar mass of dry air, MH2O is the molar mass of water vapour, XH2O is the volume mixing ratio of water vapour to dry air, q is the specific humidity and , , is the effective absorption cross section of one mole of methane at ν-frequency, p-pressure and T-temperature, which is an average value on the frequencies of the absorption cross sections computed with 4A software , , , with , calculated from Taking into account the residual differences in absorption due to the other gases, PROLID retrieves XCH4 from the estimation of the mean slant DAOD, with a correction on the optical path length depending on μ, the angle of incidence relative to nadir with IWF the integral of WF, the weighting function, to the SSE pressure, and , , To determine the sources of anthropogenic methane, it is preferable to have an estimate of its concentration averaged over the vertical with greater weight given to the boundary layer. This is why , the actual methane column defined by Equation (71), is not the objective of a mission such as MERLIN, which aims only to find the concentration defined by Equation (72). Indeed, source identification favours knowledge of a tropospheric column or even a boundary layer column, rather than a complete column, especially since the methane content of the stratosphere is largely decoupled from the sources and governed by its own physics.

Averaging
Several shots must be averaged to reduce random noise and achieve the mission objectives. Using Equation (70), the expected DAOD over Nshot is obtained by averaging the power ratio of the pulses according to Tellier et al. [17], a weighting factor is defined that gives more weight to shots with more signal ∑ Then, , the methane column averaged over Nshot, is obtained as follows A correction must be made for the variability of geophysical conditions for the different pairs of pulses during the averaging interval. This is an iterative correction formula, but only one iteration is required [17] 2 0.5 ∑ ∑ (81) Table 4 gives the order of magnitude obtained from the above equations for several important instrumental parameters of the MERLIN mission using the LIDSIM and PROLID software. Table 4. Estimated value for the detector quantum efficiency, the detector noise factor, the speckle number both for the Earth signal and the calibration path, the response function norm, the offset value and the standard deviation of the electronic noise (in digital counts or in mV), and the K parameter defined by Equations (62) and (68). In summary, a photon incident on the detector gives a mean current over a 13.3 ns sampling interval of 0.086 nA which becomes, in steady state, a potential difference at the output of the detection chain of 85.8968 μV and finally a digital count of 10.425.

Order of Magnitude of Some Characteristic Values
For a numerical simulation with a standard atmosphere, no aerosol and no clouds, and a standard deviation for the dispersion of the scattering ground levels of 15 m, the transmissions are 0.984 for the Off laser pulse, 0.297 for the On laser pulse and 0.853 for the sunlight in the spectral filter. The DAOD is then 0.61319. The window size for estimating the pulse energy is 35 samples for the monitoring pulse and is expanded to 43 samples for the ground signal due to the vertical distribution of the scattering surface. Table 5 gives the total number of photons per pulse reaching the detector, the number of sunlight photons per 13.3 ns time interval, the equivalent number of photons for the noise in the same interval (187.3 associated with ^ , the overall electrical intensity noise and 0.9 associated with ^ the voltage noise, both defined in Appendix A) and the SNR for the different recorded pulses. In the absence of any other noise, a shot-level SNR of 4.2 is required to achieve 2% accuracy on methane retrieval for a set of 140 shots. Table 5. For two ground reflectance values: the number of photons in the laser pulses (calibration, Off and On), the number per samples of photons in the solar flux and their contribution in mV to both the signal offset and the additional variability (for a sun zenith angle of 70° and 80°), the electronic noise (both for the amplified part as signal or current noise and for the differently amplified part or voltage noise) expressed in photons per sample and by its contribution in mV to the bias and dispersion, and finally, the 3 SNRs of the laser peak energy: calibration, Off and On (with and without solar flux).

Signal Validation
Numerous experiments were carried out without noise, to verify the developed software, by varying technical parameters such as: the vertical resolution and upper limit in altitude for radiative computations, vertical and spectral resolutions for transmission calculations, the temporal resolution for electronic simulations and ground reflectance representation. In all noise-free cases, the SSE (defined by Equation (59)) is retrieved to better than 0.1 m and XCH4 is retrieved to better than 0.5 ppb. Sensitivity experiments were also conducted for methane content, ground elevation, weather patterns, aerosol amounts and cloud types. Figure 17 shows the impact of methane content, atmosphere type and aerosol content on the digitised signal. These experiments use as a reference point a standard atmosphere without aerosol with a methane content of 1780 ppb. For the methane sensitivity: the methane content varies from 1760 to 1800 ppb; for the atmosphere type sensitivity: the typical atmospheres from the TIGR database are used, with the temperature and humidity profiles given in Figure 4; for the aerosol sensitivity: the five distributions recommended by ESA [25] are used (see Figure 7). To quantify the sensitivity of the digital counts to methane concentration, it can be noted that there was an increase of 40 ppb of methane results in a decrease of 23 counts out of 2788 for the peak of the On signal (no change for the Off signal), i.e., only about 2 photons less per sample. When changing the type of atmosphere, even for times with no laser signal, the signal level varies by about the same amount of 20 counts. This is mainly due to the average humidity of the atmosphere which, as it increases, decreases the solar radiation and the laser radiation for the Off beam, which is placed on a water line (see Figure 6). The variation for the On beam is differenthumidity variations are dominated by variations in the temperature profile. Indeed, at the methane multiplet where the On wavelength is located, the absorption increases when the temperature decreases (see Figure 6b). The signal is more perturbed by the type of atmosphere than by a change of 40 ppb of methane, highlighting the importance of good temperature and humidity profiles for recovering methane columns. The aerosols, for their part, can significantly attenuate the signal to a level that will not exceed that of noise when taken into account. In addition, it should be noted that at ground level, the aerosols contribute to the bias in the determination of the SSE. This is because photons backscattered from aerosols close to the ground mix with those returned from the ground itself and distort the return pulse. This effect will be investigated in a future study on the accuracy of SSE determination. For all cloud types, except polar stratospheric clouds and sub-visible cirrus clouds, the laser signal is fully scattered, so there is no echo from the ground (see Figure 18). But in the absence of noise, it is still possible to retrieve XCH4 content above the cloud layer with similar precision and accuracy as over cloudless ground. Note the importance of the change in the background signal with sunlight.  Figure 19 shows the weighting functions computed from Equation (73) for representative atmospheres from the tropics to the poles. They are mainly determined by the choice of the On and Off wavelengths. As desired, these functions are not too sensitive to the air mass and are practically uniform over a large part of the atmosphere. They decay in the stratosphere, which means that MERLIN only observes a tropospheric methane column. The maximum sensitivity is around 700 hPa. A maximum in the boundary layer would have been preferred, but it is not physically possible. Finally, for the error budget, it should be noted that the use of the vertical resolution of the ECMWF model in PROLID results in a quadrature error in the estimation of these weight functions, leading to a bias of about 0.2 ppb in the methane retrieval.

Noise Impact
Simulations were performed one hundred times for the standard atmosphere [23] with 140 random realisations of the noise each time. Figure 20 shows the shot-by-shot distribution of the DAOD retrieved in PROLID minus the actual DAOD computed in LIDSIM, as well as the distribution of XCH4 retrieved minus the actual value of 1780 ppb. As expected, with a bias of 21 ppb and a standard deviation of 303 ppb, the shot-by-shot XCH4 is not usable. Therefore, averaging several shots is necessary. Moreover, it should be noted that the distribution is not Gaussian, which implies that the averaging has to be performed with care, as explained above in Section 3.4, in line with what is studied in Tellier et al. [17].  Figure 21 shows the distribution of XCH4 averaged over 140 shots before and after subtracting the statistical bias on the DAOD. Before subtracting, the average XCH4 distribution has a bias of 30 ppb and a standard deviation of 27 ppb. After subtracting the statistical bias, estimated from the SNRs estimates (see. Equation (78)), the experimental bias is reduced to 1 ppb and the standard deviation becomes 26 ppb.

Full Orbit Simulations
CNES provided the orbital data (sun-synchronous orbits at 6.00 UTC) and ECMWF the atmospheric data for 18 June 2019. The orbit starts north of Canada, descends over the Pacific, passes over Antarctica, ascends over Africa, Tunisia, Sardinia, Corsica and the Alps, before ending over Greenland. Figure 22 shows the ground track and the DEM interpolated on it. Complete orbits of 810 cells of 140 shots were simulated, without aerosols or clouds, first without any noise and then with the instrumental noise. Some experiments were also performed with vertically uniform XCH4, but the results presented here correspond to experiments with realistic methane profiles.  Figure 23 shows the actual ground reflectivity used along the ground track. At high latitudes, the reflectivity data show some anomalies. At the same latitudes, due to the nature of the ground (snow or ice), the reflectivity is close to zero according to the reflectivity data set used. There is therefore no signal. The corresponding points are eliminated in the processing.  Figure 24 shows along the ground track: in green, the full column average of XCH4 computed (Equation (71)) as the ratio of the number of moles of methane to the number of moles of air; in red, the XCH4 reference computed (Equation (72)) as the vertical average concentration obtained by applying the MERLIN weighting function; in black, the XCH4 retrieved shot-by-shot provided by PROLID (Equation (81)). As mentioned earlier, the objective of the MERLIN mission is to retrieve the XCH4 reference, not the full column average of XCH4. Without noise, the retrievals are almost perfect. For this reason, the red curve is hidden by the black dots in Figure 24a. This is neither the object of the Merlin mission nor of this article, but one can observe the differences and their geographical variability between the full column of methane (green curve) and the reference methane that MERLIN will be able to observe (red curve hidden by the black points). For the simulation that was carried out, the differences are positive in the Northern Hemisphere and negative in the Southern Hemisphere, with the major differences being over the relief. Figure 24b shows the highly noisy nature of the shot-by-shot inversion (note that the green and red curves are the same on Figure 24a Table 6 summarises the results on SSE and methane content, shot-by-shot and also averaged over 140 shots (7 s). There are less data for the experiment with instrumental noise than for the one without, because, for points with low reflectivity, the signal is lost in the noise, and it is then impossible to process these points. The determination of the SSE requires only the Off pulse, whereas the XCH4 also requires the On pulse. Therefore, the number of shots considered for the SSE and for the XCH4 is not the same. For the 7 s averages, only those for which all 140 shots are available are considered. This will be improved later on.
The standard deviation of 16 ppb for these 7 s average data is better than the 26 ppb presented in the previous section for the tests with a standard atmosphere, just because here the reflectivities used (see Figure 23) are on average higher than the value of 0.1 taken in the tests of the previous section. It should be noted that this standard deviation clearly depends on the number of shots available for averaging, and will therefore increase when aerosols and clouds are considered. Table 6. For the experiments with and without instrumental noise, per shots and per 50 km cell, the differences between retrieved SSE and XCH4 and their references in terms of mean bias and mean standard deviation calculated over the number of points for which XCH4 can be retrieved (points with sufficient signal due to sufficient reflectivity).

SSE (in m)
XCH4 (in ppb) For the experiments with instrumental noise, Figure 25 shows, for each 50 km cell, the differences between the retrieved SSE and XCH4 and their references. In the experiments without noise, a slight modulation of errors was observed as a function of the ground altitude and the temperature profile. However, the experiments with instrumental noise give homogeneous results in quality over the whole orbit, except of course for a strong dependence on the target reflectivity. This is consistent with the MERLIN random noise budgets, established by Airbus and CNES, which show that the instrumental noise is predominant [2]. . XCH4 retrieved by PROLIDL2, the actual XCH4 and the reference XCH4 (green and red curve as in Figure 24) also averaged over 140 shots.

Experiment
Various sources of error are not taken into account in this study, such as non-linearities of the laser emission or detection, or the lack of knowledge of certain parameters: in particular, the emitted wavelengths, the satellite attitude or the state of the atmosphere. The main source of error not taken into account is that induced by uncertainties in operational meteorological analyses. Nevertheless, preliminary studies (not detailed here) have indicated that it is less than 0.9 pp. Thus, although there is no geophysical noise in these simulations (the weighting function is computed with the actual values used to simulate the MERLIN data in LIDSIM), the present study can conclude on the performance achieved by the MERLIN mission. Figure 26 shows the performance of the modelled instrument compared to that expected by the users. Figure 26. MERLIN performances estimated with LIDSIM and PROLID (the red/blue/green curves correspond to the threshold/breakthrough/target performances, respectively). The cyan dots represent the tests described in Section 4.3 (at the shot and cell level before and after removing the statistical bias), and the orange diamonds represent the results obtained in this section (for one shot and for one cell).
The target performance does not seem to be achievable as the speckle noise alone gives a standard deviation of 60 ppb for one shot and 5 ppb for a 50 km average [54]. However, breakthrough performance is reached and optimisation of the processing parameters can further improve it.

Conclusions
The LIDSIM and PROLID software packages have been developed and tested, respectively, to generate MERLIN lidar signals and to produce XCH4 data from them. They are proving to be powerful tools for preparing satellite missions such as MERLIN. They have been developed with enough generality to be applicable to other IPDA lidar missions and many parts of these packages can be useful for preparing other types of satellite missions or even for simulating ground-based remote sensing instruments. In this line of work, PROLID has already been used on CHARM-F data [19] provided by the DLR.
In the framework of MERLIN, numerous sensitivity tests have been performed on both technical parameters and geophysical conditions. Orbital simulations investigating the impact of the clouds and aerosols remain to be carried out, as well as the updating of the detection chain in order to follow the choices of the instrument manufacturer (Airbus). In the future, the level of performance estimated by our tools will be systematically compared with that estimated by the performance model used by Airbus and CNES to validate the design of the MERLIN instrument. LIDSIM and PROLID will also be used to optimise the pre-processing of digital counts and to study the possibility of retrieving cloud and canopy heights from MERLIN data (secondary products).
The simulations carried out so far and presented in this paper confirm that the current MERLIN design, as simulated here, achieves the ambitious breakthrough performance set for the mission. In practice, any instrument measures the one-sided power spectral density (OPSD), so the value provided must be divided by two when the power spectral density (PSD) is required. In particular, white noise is characterized by its constant amplitude spectral density (ASD), with the following relation Electrically, the APD is characterized by a resistance Rd = 1 MΩ and a capacitance Cd = 2.5 pF, while the TIA is characterized by an open loop for the operational amplifier with a continuous gain A0 = 1778 and a cut-off frequency = 230 MHz, and by a resistance Rf = 1 MΩ and a capacitance Cf = 0.2 pF. Current and voltage noise are considered as white noise characterised by their ASD, the APD dark current îd = 1.3 fA Hz −0.5 , the TIA intensity noise î0 = 1.3 fA Hz −0.5 , the thermal or Johnson-Nyquist noise in the shunt resistance ^ 4 (Tf = 280 K) and the TIA voltage noise ûn = 7 nV Hz −0.5 . The output voltage of the amplifier ^ is modelled as follows where ^ is the voltage at the point P (see Figure A1) and is the voltage noise. According to Kirchhoff laws, the sum of the currents at any node is zero. Applying this to node P (see Figure A1) and taking into account the noise sources, the following relationship is obtained and V/A/s 2 , respectively. In practice, the inverse Laplace transforms are computed separately for the TIA and the AAF using the partial fraction decomposition, since then, for each fraction, the inverse Laplace transform is known analytically. Ri ( Figure A2a ) and Ru ( Figure A2b) are then computed by convolution.