Impact of Meteorological Uncertainties in the Methane Retrieval Ground Segment of the MERLIN Lidar Mission

: MERLIN (MEthane Remote sensing LIdar missioN) is a Franco-German space mission designed to provide weighted columns of atmospheric methane through an inversion of the lidar signal using a priori information on the atmospheric state. Uncertainties about the meteorological parameters of the observed scene used in the ground segment contribute to the error budget on the retrieved methane column. With the LIDSIM (LIDar SIMulator) data simulator and the PROLID (PROcessor LIDar) inversion processor developed for MERLIN, we perform an impact experiment using ECMWF (European Centre for Medium Weather Range Forecast) ensemble forecast data. In addition, we estimate the standard deviation of the error in the methane column due to the meteorological uncertainties to be about 0.6 ppb. In addition, we innovate by discussing the impact of interpolations both in time and space, focusing on vertical extrapolations under the topography by using state-of-the-art methods to determine from the scatter between these methods the range in which the actual proﬁle should be. We conclude that, in areas where the topography variations exceed 10 m over 10 km, an additional random error of 0.1 ppb is due to our lack of knowledge of the adjustment of atmospheric proﬁles to terrain. Finally, we point out that further work needs to be performed on temporal interpolation. Indeed, the 3 h time interpolation of atmospheric tides can create regional biases of up to 2 ppm (which is a major problem for models trying to identify methane sinks and sources).


Introduction
Even if not very abundant in the atmosphere (less than 1.9 ppm), but with 60% of its emissions being anthropogenic, methane (CH 4 ) is responsible for 20% of the global warming produced by all greenhouse gases so far [1]. With a global warming potential 28 times larger than carbon dioxide (CO 2 ) over 100 years [2], methane is the second most important greenhouse gas contributing to human-induced climate change. Its monitoring on a global scale is a challenge to which Germany and France are contributing through the development by their respective space agencies, DLR (D eutsches Zentrum für Luft-und Raumfahrt) and CNES (Centre National d'Etudes Spatiales), with the MERLIN (MEthane Remote sensing LIdar missioN) mission planned for the end of this decade [3]. The MERLIN mission uses the CNES Myriade Evolution platform developed by Airbus France; meanwhile, Airbus Germany is developing the MERLIN instrument: an IPDA (Integrated Path Differential Absorption) lidar [4]. The data processing is described in the MERLIN

Meteorological Data in IPDA Lidar Methane Retrieval
An IPDA lidar such as MERLIN emits two beams with close ν On and ν O f f wavenumbers but very different methane absorptions. It records the signals corresponding to the energies emitted in these beams and those backscattered by the ground or a cloud [9]. The ground processing of these records allows us to estimate the DAOD f ull (Differential Atmospheric Optical Depth over the entire optical path) from the normalised ratio between the energies of the two beams and the SSE (Scattering Surface Elevation) from the assessment of the time elapsed between the emission and reception of the signals (for a description of how such estimates are made, see MERLIN ATBD [5] or Cassé et al. [6]). Then, XCH r 4 , the mole fraction of methane to dry air averaged over the atmospheric column, is retrieved from the following relations derived from the Beer-Lambert law: with, for any G-gas and any p-pressure: and where X G (p) denotes the mole fraction profile with respect to dry air, g(lat, p) the gravity at the target location (determined by its latitude and pressure), M d the molar mass of dry air, T(p) the temperature profile, q(p) the specific humidity profile and σ e f f G (ν 0 , p, T) an average (over the spectral distribution of the laser emission around wavenumber ν 0 ) of the absorption cross sections [6]. The absorption cross sections are computed from the GEISA2015 spectroscopic database [10] with specific improvements made by work carried out in the framework of the MERLIN mission [11,12].

Meteorological Data Availability
The meteorological data required by the computation of WF G (p) (the weighting functions) as well as of P SSE (the pressure at the SSE) are provided by operational analysis such as those of the European Centre for Medium-Range Weather Forecasts (ECMWF), accessible in the Meteorological Archival and Retrieval System (MARS) [13]. The data we need are available there every three hours on the native model grid from the 4D-VAR assimilations performed twice a day on the 00 UT and 12 UT networks [14].

Time and Horizontal Interpolations
For each vertical level of the meteorological model, a linear interpolation in time and a bilinear horizontal interpolation provide the temperature and humidity at the MERLIN target location. Similarly, Φ NWP (the interpolated model topography) and Ps NWP (the interpolated surface pressure) are obtained at the target point. Figure 1a shows a diagram of these time and space interpolations. Figure 1b is an illustration of the density of the points where meteorological data are available and of the target points over the Corsica region. Bilinear interpolations in longitude and latitude are an extremely simple, cheap and efficient implementation for horizontal interpolations. It avoids the discontinuities encountered with the use of the nearest point and requires only four original points. On each direction, f a (x) is the interpolated value at the location x defined with k = x ∆ x by: where ∆ x is the distance between the points and f k is the value of the function at the location k ∆ x . Bilinear interpolations accuracy for a particular field can be estimated by comparing their expression f a with the exact value f e defined by its expression as a linear combination of spherical harmonics Ψ i which form an orthogonal basis for the decomposition of the different horizontal fields in the NWP models [16].
, Equation (5) shows the accuracy of the linear interpolation as the factor term of the a i is at least as ∆ x 2 : This means that the interpolation error is maximal in the middle of the interpolation interval and depends on the second derivatives of the basis functions at the sample points (i.e., their curvature at these points), and the error decreases quadratically with the size of the sampling interval.

Different Methods
Vertical interpolations are performed to adapt the meteorological data to Φ DEM : the actual topography provided by a DEM (Digital Elevation Model). In order to determine Ps DEM , the pressure corresponding to Φ DEM (the ground elevation at the target point), on the one hand, and the temperature and humidity values at the predefined pressure levels as a function of Ps DEM , on the other hand, interpolations are used between the vertical levels. These interpolations are chosen linear in pressure for temperature and specific humidity, linear in log pressure for geopotential and reciprocally linear in geopotential for log pressure. They are performed in terms of departures from the standard atmosphere (Appendix A) for increased accuracy [17]. Extrapolations are required to extend the profiles below the model topography when the target is below.
Four extrapolation methods are considered in this paper. Three of them differ in the gradient used to extrapolate the temperature and humidity profiles. For the first one, the gradient is chosen to be zero: it is assumed that the values of the layer just above the model floor are representative of the layer below the ground up to the target altitude. The last layer just at the level of the model ground is extended up to the actual ground. For the second one, the gradient is computed from the difference between the lowest layers of the atmospheric model. The extension is carried out taking into account the gradient between the second last and the last layer.
These first two extrapolations are physically relevant only if the difference between Φ NWP and Φ DEM is small. Unfortunately, this is not the case. Figure 2 shows the differences in topography for a full orbit between the ECMWF model at a resolution of about 10 km and the EarthEnv data at a resolution of 90 m. In a region such as the Alps, where the relief varies greatly, the differences can be over 1 km, while for flatter regions, the differences remain in the order of 20 m. It should be noted that over the sea, the differences are generated by Gibbs waves due to the spectral fit of the model topography, but they do not exceed one metre. For the third extrapolation method, detailed in Appendix B, the temperature gradient below the ground is that of the standard atmosphere and the specific humidity gradient is zero. This third method is close to the way operational meteorological centres proceed [18] to post-process their data. When the difference in topography is large, it is better to use a standard value for the gradients instead of zero or an estimate computed very near the ground.
The fourth method, described in Appendix C, uses an unpublished approach developed at Météo-France for interpolations between different resolutions of the operational model ARPEGE (Action de Recherche Petite Echelle et Grande Echelle), hence its name APACHE (ARPEGE Pour ARPEGE CHamp Echangé) [19]. Its objective is to preserve a boundary layer close to the ground, both when the target topography is lower than the initial topography and in the opposite case. The idea is to retain certain quantities in the near-ground layers that are characteristic of the physical interactions between the atmosphere and the ground. The APACHE method preserves the potential temperature gradient related to the stability of the atmosphere and the relative humidity which is assumed to be more independent of temperature than the specific humidity and more characteristic of the balance between physical processes. Table 1 summarises the main features of the four vertical interpolation methods used. To document the differences between the different vertical extrapolation, let us examine the temperature and humidity profiles produced by the different extrapolation methods in two cases, one above Africa and the other above the Alps, when the target elevation is set arbitrary 200 m below (Figure 3a,c) or above (Figure 3b,d) the model topography.
When the actual ground is below the model ground, the use of a non-standard gradient leads to unrealistic temperatures extrapolations near the actual ground. As expected, only the APACHE approach maintains an inversion near the actual ground whether it is above or below the model ground. Similarly for humidity, Figure 4 shows the different methods of extrapolations. There is a risk of creating oversaturated humidities when using the last gradient (negative humidity can also occur in that case) and more generally when the specific humidity is maintained but the temperature is decreased. This is why APACHE maintains the relative humidity rather than the specific humidity. Finally, the APACHE approach again shows its advantage in preserving the structure of the boundary layer near the ground.

Intrinsic Uncertainties in Weather Forecasts
Furthermore, regardless of these interpolation issues, our knowledge of the meteorological fields is subject to errors that are generally estimated from the statistics of the differences between the model and the observations. A CNES study [20], based on ECMWF data from 2006, estimates for surface pressure, a bias of 0.2 hPa and a random error of 2 hPa; for temperature, a bias between −0.2 K and 0.2 K (varying with vertical) and a random error of 2 K; and for specific humidity, a bias of 4% and a standard error of 15%. With these estimates, in the case of MERLIN, Bousquet et al. [1] attribute to the NWP data 7 ppb (8%) out of the 33 ppb estimated for the total random error and 0.7 ppb (6%) out of the 6.4 ppb estimated for the total systematic error. However, this statistical approach neglects error correlations both between vertical levels and between variables. For this reason, we have chosen another approach based on ensemble forecasting tools [21]. ECMWF generates every 12 H a set of fifty analyses compatible with the uncertainties of the forecast used as a first guess and those of the observations available [22]. These fifty perturbed forecasts were used in this study to provide an estimate of the uncertainty of the actual state of the atmosphere with respect to the unperturbed initial-state forecast used as reference. With our LIDSIM software, we have simulated the MERLIN signals from the atmospheric data of the reference forecast. Then, we have used each of the fifty perturbed forecasts in turn to obtain a estimate of the mole fraction of methane with respect to dry air averaged over the atmospheric column with the PROLID software, by computing WF G (p), the weighting functions for each gas, following Equation (3).

Time Interpolation and Tide Waves
We have decided to use the weather data available every 3 h in the MARS database and to interpolate them linearly. However, it has been highlighted by work carried out by the CNES [20,23] that using model outputs at the hourly frequency leads to systematic differences between the pressure fields predicted at time H and interpolated between H − 2 and H + 1, whether using linear interpolation or a cubic spline. The differences interpreted as being related to the representation of atmospheric tides may induce regional biases up to 3 ppm on the mole fraction of methane with respect to dry air averaged over the atmospheric column (see Figure 5). Indeed, time interpolation creates systematic errors for insufficiently sampled periodic components such as atmospheric tides, as it is well documented in other contexts [24].

Impact on the Methane Retrieval of Vertical Interpolation Methods
We used several methods to adapt the atmospheric profiles corresponding to a particular ground altitude to another ground altitude. With LIDSIM, we compute the temperature and humidity profiles from the ECMWF data for a vertical pressure grid defined from the pressure value estimated at the DEM value, using one of the interpolation methods described above and detailed in Appendices B and C. In addition, we use these profiles to simulate the MERLIN data. Then, with PROLID, we invert the simulated signals using profiles constructed with another method. In all cases, except the "last gradient" case, a "zero gradient" method is used to extrapolate data from the DEM value to the SSE. There is no reason to use a more sophisticated method because the difference between the DEM value and the SSE is small except in the cloud cases, but then the boundary layer should not be moved as the SSE does not represent the ground.
The calculation of the IWF (Integrated Weighting Function, see Equation (3)) for a particular temperature and humidity profile, even if the surface pressure is fixed, depends on the distribution of the levels over the vertical. As explained in Appendices B and C, the "standard gradient" and APACHE methods change the discretisation of the profiles above the pressure determined at the DEM altitude (they resample atmospheric levels). However, the implementation of the "last gradient" and "zero gradient" methods is performed without level resampling. The deviation of the IWF due to the resampling could be measured. Expressed as a maximum error in the mole fraction of methane to dry air averaged over the atmospheric column, it corresponds to 0.08 ppb. This value is in agreement with an estimate we made (not shown) on five representative atmospheres from the TIGR database [25] and which shows that the difference between the exact calculation of the IWF and its approximation with a discretisation on 137 levels (vertical resolution of the ECMWF data used) is of the order of 0.04 ppb. This order of magnitude is also the contribution of the region above 40 km which is not taken into account in LIDSIM although its contribution will be present in the satellite data.
For the two points used for Figures 3 and 4, the methane Weighting Function (see Equation (3)) is computed for the various interpolations both for a target above and below the model topography (see Figure 6). By using a different extrapolation/interpolation method in LIDSIM and in PRO-LID, it is possible to estimate the order of magnitude of the error due to these extrapolations/interpolations. For the points considered above, whether the target is above or below the model topography, variations in the IWF (Integrated Weighting Function) induce uncertainties on the mole fraction of methane to dry air averaged over the atmospheric column of no more than 0.32 ppb. Table 2 shows results for a complete orbit, the one whose topography on its ground track has already been shown in Figure 2a. Table 2. The standard deviation and the absolute maximum of the difference on the mole fraction of methane to dry air averaged over the atmospheric column in ppb according to the extrapolation methods used by LIDSIM and by PROLID. The mean is in all cases the same and its value of 0.08 ppb can be attributed to factors other than the meteorological data. The simulation is for a complete orbit. It is based on meteorological data from 6 May 2019 between 3:30 UT and 5:08 UT and on realistic reflectivity data (for more details, see Cassé et al. [6]) but assumes a constant value of 1780 ppb for methane content. Although there is no simulation of instrumental noise, due to the limited numerical accuracy when encoding the file representing the data transmitted to the ground segment, filtering of shots with a signal-to-noise ratio greater than 1 is necessary to eliminate the shots for which there is no signal. The results are optimal when we use the same method in LIDSIM and in PROLID. This simply means that the more our method matches what is happening in our environment, the better the results. Using the "last gradient" method in LIDSIM leads to thermodynamic profiles with extreme values that cannot be obtained with the other approaches.

Extrapolation in LIDSIM with
The discrepancies are located in areas where the topography is highly variable. Otherwise, there is little difference in topography between the DEM and the meteorological model, which generates no significant difference in the methane estimates. Focusing on situations where the topography gap is large, the "standard gradient" method is preferred to the "zero gradient" method for PROLID because the latter gives strongly different values, as expected. Figure 7a shows, as a function of the difference in topography between the meteorological model and the DEM, the differences between the mole fraction of methane to dry air averaged over the atmospheric column values retrieved by crossing between LIDSIM and PROLID the "standard gradient" method and the APACHE method. The error due to interpolation can be estimated, even if the real profiles are not known, from the dispersion of the results given by the different interpolation methods. This dispersion leads to differences of 1 ppb on the mole fraction of methane to dry air averaged over the atmospheric column when the topography correction reaches 200 m. However, examination of their distribution (Figure 7b) shows that the vertical interpolation error can be estimated at 0.0 ± 0.1 ppb.
(a) (b) Figure 7. Distribution of the differences between the mole fraction of methane to dry air averaged over the atmospheric column retrieved using the "standard gradient" in LIDSIM and APACHE in PROLID and the same quantity retrieved using APACHE in LIDSIM and the "standard gradient" in PROLID. To summarise, the impact on methane retrieval of the uncertainties in the way surface pressure as well as temperature and humidity profiles are adapted to another topography than the model topography is less than 3 ppb everywhere (even when the change in topography is more than 1 km). This value is obtained from Figure 7a before it is truncated at topographic deviations of ±200 m. The additional standard deviation is less than 0.1 ppb but with a systematic part at a particular location.

Sensitivity to the Analysis Errors
As explained above, we simulated the MERLIN data in LIDSIM using the ECMWF description of the atmosphere from the unperturbed initial state and inverted these data in PROLID using the other fifty descriptions of the atmosphere corresponding to forecasts from the perturbed states. We used data from 6 May 2019 between 03 UT and 06 UT. We present the results for the same orbit as above. In LIDSIM, we did not add instrumental noise to see only the impact of the NWP data. In LIDSIM and PROLID, we used the standard method to adjust the data to the DEM value. Figure 8 shows the mean and the standard deviation of the fifty differences between XCH r 4 (the mole fraction of methane to dry air averaged over the atmospheric column) obtained with Equation (1) for a set of meteorological parameters and its reference value [6]. This reference is defined as the pressure-averaged actual methane content provided by ECMWF with the weighting function WF CH 4 (p) calculated for the set of meteorological parameters from the unperturbed initial state. The mean bias is of −0.03 ± 0.39 ppb and the mean standard deviation is of 0.58 ± 0.27 ppb. Along the orbit, the dispersion of the differences in surface pressure ( respectively in total water vapour) with respect to the values of the forecast from the initial unperturbed state is shown in Figure 9 (respectively in Figure 10).
The fifty initial perturbations at 00 TU are symmetrical in pairs, but the integration of the model up to the observation time introduces many non-linearities that break this symmetry. This is why the average difference is not zero. Over the whole simulated orbit, for the surface pressure, the mean bias is 0.95 ± 12.65 Pa and the mean standard deviation is 17.96 ± 9.48 Pa, and for the total amount of water vapour expressed in Pa, the mean bias is −0.38 ± 7.22 Pa and the mean standard deviation is 6.82 ± 5.62 Pa.  We calculated the correlation coefficient (r) between the methane error standard deviation series (X(i) i=1,N ) and the pressure error standard deviation series (Y(i) i=1,N ).
We obtained a correlation factor of 0.86. Then, we compute the correlation coefficient between the methane error standard deviation series minus the pressure error series multiplied by the previous regression coefficient and the total water vapour error standard deviation series. This new correlation is also 0.86. These correlations can be interpreted in terms of explained variance. The variance from the meteorological uncertainties in XCH 4 is explained to 50% by the uncertainty in the surface pressure and to 25% by that in the water vapour. The rest is due to variations in temperature and vertical humidity distribution.

Forecasting Errors
In order to compare our results regarding the impact of NWP data accuracy on the methane retrieval accuracy with those of Bousquet et al. [1], one must take into account the large difference in the estimates of systematic and random errors in the surface pressure used as input: 1 Pa and 18 Pa (see Figure 9) in this study instead of 20 Pa and 200 Pa in Bousquet et al. (i.e., almost 10 times less). Revisiting the 2006 data on the differences between observations and meteorological analyses, which were used to determine the Bousquet et al. values, we find that: globally at sea, compared to buoy measurements in the northern hemisphere, the pressure bias was 0.2 Pa and the standard deviation was 44 Pa, while on land, compared to aeronautical station data, the pressure bias was −11 Pa and standard deviation 60 Pa. Therefore, the values established in the CNES study [20] and used by Bousquet et al. [1] seem pessimistic. In 2011, Kiemle et al. [26] state that the random pressure error is less than 100 Pa in agreement with the Dee et al. [27] estimate around 70 Pa. The different values can be also explained by the regular improvement of forecast quality which continues up to now. The data available in the ECMWF website indicate an analysis bias of 3 Pa and a random error of less than 50 Pa for these last years. The estimate derived from the ensemble forecast is therefore representative of the state of the art, especially as our single orbit assessment necessarily underestimates the variability somewhat.
Given the advances in numerical weather prediction, the random error estimated by Bousquet et al. [1] at 7 ppb for the mole fraction of methane to dry air averaged over the atmospheric column can be reduced to 0.6 ppb and the bias from 0.7 ppb to 0.4 ppb, i.e., a gain of nearly 8% over the total random error and 3% over the total systematic error. However, we have to consider the differences between the actual observed profiles and those that can be estimated from the numerical forecasting data, which are not accounted for in this experiment. This is why we discuss the contribution of the various interpolations in the following.

Horizontal Interpolations
For our application, bilinear interpolations at a fixed vertical level of the NWP model are satisfactory. Although the quality of such a simple approach is often questioned, as there are many other ways to perform these interpolations, bilinear interpolations offer an excellent compromise between simplicity and accuracy. Alternatively, polynomial interpolations can be used to ensure multi-order derivability. However, they may be expensive to determine and subject to oscillations as their degree increases (Runge phenomenon). It is often preferable to use splines that are piecewise polynomials and to increase the number of polynomial pieces instead of increasing the degree of the polynomials used. Linear splines are obviously identical to the linear interpolation performed. Quadratic splines are rarely used in practice because with three degrees of freedom per interval, in dimension one, they require an additional constraint per interval which is not easy to fix objectively and the results are very dependent on this choice. Only cubic splines are used in certain circumstances. However, in their basic version, they require, in dimension two, sixteen points to determine the required first derivatives by finite difference, and they present problems with preserving extrema which can lead to negative values for humidity, for example. Nevertheless, it should be noted that where high accuracy is required, as in the case of using the semi-Lagrangian technique in a numerical model, there are several ways to save computational time and ensure extrema preservation [16]. However, there is no reason to use here these elaborate methods because of the inherent uncertainties in the quantities processed and because we have not identified any case in which a systematic error can occur.

Vertical Extrapolations
For the vertical extrapolations, systematic deviations may occur related to the method applied. Figure 7 shows that even for small differences in topography, different postprocessing of the meteorological data can lead to deviations of 0.1 ppb in the vertically averaged methane concentration retrieval. Although it will be possible, on the basis of Figure 7a, to filter out the most important uncertainties by using a threshold on the topography differences between the meteorological model and the DEM; we assume a standard deviation error of 0.1 ppb due to these vertical interpolations.

Time Interpolations
The main concern with time interpolation is to avoid systematic errors due to periodic forcings with periods shorter than one day which are poorly represented by linear interpolation between two instants. We have identified the diurnal and semi-diurnal atmospheric tides and using a three-hourly sampling is not a sufficient representation of these waves. In the processing of altimeter data for the calculation of the pressure correction, the mean atmospheric tides are subtracted from the pressure field before time interpolation and its component calculated exactly at the target time is then added [28]. Unfortunately, this approach is more difficult to implement in our case as we not only have to interpolate the pressure field but also the temperature and humidity fields which are 3D fields. We plan to investigate how to deal with this problem in the future, as further research on this issue is needed, with the challenge of minimising local biases on retrieved methane.

Conclusions
We would have liked to recommend for the MERLIN ground segment to use the ECMWF operational analysis data available every three hours and to interpolate them linearly in time. However, we cannot do this until a method is validated to deal with atmospheric tides because of the regional bias they may induce. In addition, we need to plan specific work on this issue. There is no specific problem expected with the use of bilinear interpolations on the horizontal surfaces of the model. For the vertical extrapolation, only the "standard gradient" method or the APACHE method should be used in the MERLIN data processing due to the possible large differences between the topography of the meteorological model and the actual topography. The APACHE extrapolation method appears to be physically more consistent, but using simulated data, it is not possible to demonstrate a significant contribution to MERLIN data processing. Work will need to be performed during the flight commissioning phase, especially in areas with highly variable topography to analyse methane variations according to the interpolation method used. We plan to study this point with airborne data from the CHARM-F instrument [29] deployed during the MAGICS campaign [30].
With the exception of the atmospheric tides problem, meteorological uncertainties to the values used in the ground segment are expected to contribute at most 0.4 ppb to the bias in methane retrieval from MERLIN data and 0.7 ppb to its standard deviation, while uncertainties related to instrumental noise for 50 km averages of the mole fraction of methane to dry air averaged over the atmospheric column correspond to a bias of 1.2 ppb and a standard deviation of 16 ppb [6]. With appropriate atmospheric tides processing, the MERLIN instrument will be therefore at the breakthrough performance level of 2 ppb for systematic errors and 18 ppb for random errors. with p the pressure, n the number of moles per m 3 , R = N A k B the universal gas constant (calculated from N A the Avogadro number and k B the Boltzmann's constant) and T the temperature; on the other hand, the hydrostatic relation: with M d the molar mass of the dry air and Φ the geopotential which is the energy level in the gravity field defined by g its intensity along the vertical (dΦ = gdz).
In order to facilitate their interpretation, the values of the temperature gradient with respect to the geopotential are generally expressed using as a unit of energy per unit of mass the geopotential metre defined by 1 mgp = g 0 J/kg where g 0 is a constant fixed to 9.80665.
With this unit, the values of the gradient Λ l with respect to the geopotential are indeed similar to those of a temperature gradient per geometric metre. We then have, in Table A1, for each layer: [Φ l−1 , Φ l ], its geopotential limits, and Λ l , its constant temperature gradient in K per kmgp:  [32]: On a vertical, we have Φ NWP (the ground geopotential), Ps NWP (the surface pressure) and {T (i) , q (i) } i=1,N (the temperature and humidity profiles) on the "full levels" {p (i) } i=1,N : {Φ (i+1/2) } i=0,N (the geopotential profile) is calculated on the half levels by integrating the hydrostatic relation from Φ NWP : where R is the universal gas constant, M, M d and M w the molar masses of the moist air, the dry air and the water vapour, respectively, and g 0 = 9.80665 ms −2 the gravity acceleration. First, ln(Ps DEM ) (the logarithm of the pressure at Φ DEM , the target altitude) is obtained by linear interpolation of the deviation of the logarithm of the pressure from the standard atmosphere, calculated at Φ (i) . With the following notation for the linear interpolation at the Y level of the variable defined by the X (i) at the Z (i) level: intlin Y,Z (i) where P std (Φ) is the standard pressure at Φ (see Appendix A). It is then possible to determine T DEM (the temperature at the target altitude). If Under the model topography, Λ (a constant temperature gradient with respect to the geopotential) is used to determine T DEM (the temperature) and Ps DEM (the pressure) from T * (the model surface temperature): for Φ DEM < Φ NWP : Thirdly, {T DEM (j) } j=1,N and {q DEM (j) } j=1,N (the temperature and humidity profiles) are computed at {p DEM (j) } j=1,N by linear interpolation according to pressure:

Appendix C. APACHE a Method for Maintaining the Boundary Layer
This method was developed at Météo-France in the 1990s from an idea by Jean-Francois Geleyn, but it has never been published except partially in Radmila Brožková's doctoral thesis, written in Czech, defended in 1995 at the University of Toulouse, in cotutelle with Charles University in Prague [19]. First, for the given value of Φ DEM , Ps DEM (the pressure corresponding to this geopotential) and {T DEM (j) , q DEM (j) } j=1,N (the vertical profiles of temperature and humidity) sampled over the {p DEM (j) } j=1,N are computed according to Appendix B, using the "standard gradient" option and the model data: Φ NWP , Ps NWP , {T(i), q(i)} i=1,N available at the {p (i) } i=1,N , as well as T * (the model surface temperature).
A second set of values is then determined by preserving for the layers defined from {p DEM (j+1/2) } j=0,N : the potential temperature gradient between a layer and the ground and the relative humidity of the initial layers. The idea is to preserve for the new layers the stability properties between the successive layers of the initial set. As the potential temperature is preserved during an adiabatic transformation, the adaptation of the initial profile to a new topography is similar to an adiabatic transformation. The relative humidity is preserved because the different physical processes are considered to lead to an equilibrium state with a particular relative humidity rather than a particular specific humidity, the value of which can vary strongly with the temperature of the environment. The following equations allow us to calculate Θ (the potential temperature) and Hu (the relative humidity): We have simplified the original approach by using only this saturation vapour pressure with respect to liquid water, instead of a linear combination, as a function of temperature, with the saturation vapour pressure with respect to ice.
Θs NWP and Θs DEM , respectively, the potential temperature of the model and target ground, are computed from Ps NWP and T * , and Ps DEM and T DEM , available according to the approach in Appendix B. Then, from the conservation of the potential temperature gradient between the two sets of layers, Θ DEM2 (i) (the potential temperature of the layer i) is determined: which allows to derive {T DEM2 (j) } j=1,N , a second temperature profile on the {p DEM (j) } j=1,N levels. In addition, at the same time, the conservation of the relative humidity between the two sets of layers gives the following equation for the specific humidity: in such a way as to keep the former in the free atmosphere (above P top the pressure of the top of the upper boundary layer) and to keep the latter close to the target surface (below P bot the pressure of the top of the lower boundary layer), while passing from one to the other in a smooth and regular manner. To achieve this, we set ∆ PBL = 175 hPa to define P top and P bot , and to avoid discontinuities, we ensure that the transition zone is not too small: P top = min(Ps NWP , Ps DEM ) − ∆ PBL P bot = max(max(Ps NWP , Ps DEM ) − ∆ PBL , min(Ps NWP , Ps DEM )) (A26) And the coefficients of the linear combination are defined as follows: : This advanced method leads to more realistic profiles but still has the disadvantage of changing the geopotential even far from the ground.