Inverse problem for retrieving greenhouse gas fluxes at the non-uniform underlying surface from measurements of their concentrations at several levels

The study focuses on the formulation, analysis, and solution of the remote sensing inverse problem to retrieve surface carbon dioxide (CO2) fluxes from measurements of CO2 concentrations at different levels within the atmospheric boundary layer. A three-dimensional hydrodynamic model of turbulent greenhouse gas (GHG) transport was used as a forward model to link the surface GHG fluxes to the drone observations of GHG concentrations. The 3D model provides a GHG concentration distribution by solving the diffusion-advection equation using information on wind speed, its direction, and turbulent exchange coefficients. The surface GHG fluxes are considered as a boundary condition. The spatial distributions of wind speed and turbulence coefficient “for a moment in time” are computed from the relaxation problem for the averaged Navier-Stokes and continuity equations, using a 1.5 order closure scheme (E-ω model). The inverse problem is to retrieve a surface GHG flux by minimizing the difference between the measured and modelled concentrations at several levels. The algorithm was applied to estimate CO2 fluxes over a non-uniform forest canopy at the Roshny-Chu experimental site in the foothills of the Greater Caucasus (Chechen Republic). To test the forward numerical problem, data on surface topography, vegetation height and density, spatial distribution of photosynthetically active solar radiation, as well as data on plant photosynthesis and soil CO2 fluxes were used.


Introduction
Due to the increasing anthropogenic influence on the atmosphere, the problem of investigating and assessing greenhouse gas (GHG) fluxes between the land surface and the atmosphere has become very challenging [1][2][3].Among the various experimental methods used to measure GHG fluxes in the field, eddy covariance and chamber methods are the most widely used [4][5].However, due to numerous methodological limitations, these techniques can only be used at selected experimental sites.In IOP Publishing doi:10.1088/1742-6596/2701/1/012141 2 particular, the eddy covariance method requires horizontal homogeneity of the underlying surface, fully developed turbulence, negligible divergence and convergence of the airflow, small fluctuations in air density, etc. [6][7].It cannot be used for representative flux measurements in areas with nonuniform vegetation distribution and complex topography.The chamber method can be used to measure soil and leaf GHG fluxes at experimental sites with highly heterogeneous topography and mosaic vegetation, but it does not allow flux averaging over a large area without a large number of simultaneous measurements on different experimental plots [8].Considering the potential limitations of different in situ methods for regional flux estimation, remote sensing in this case can be a very useful tool for flux estimation over non-uniform land surfaces.Over the past decades, numerous inverse modelling techniques have been applied to derive GHG fluxes, especially CO 2 , from remote sensing data [9][10][11][12][13][14].However, the novel techniques and models for GHG assessment using remote sensing information are still very much needed.
In our study, we propose a novel approach for regional estimation of CO 2 fluxes from UAV-based (based on unmanned aerial vehicles) measurements of GHG concentrations at two levels above the ground surface and a three-dimensional (3D) hydrodynamic model.The model allows the calculation of the stationary GHG concentration distribution over a non-uniform surface with known GHG fluxes at the lower boundary of the computational domain.The unknown GHG flux is determined by solving the inverse problem of minimising the difference between measured and modelled GHG concentrations at two levels above the ground surface.

Formulation of the full forward problem for GHG concentration calculation
To describe the spatial distribution of GHG concentrations within the atmospheric boundary layer over a non-uniform surface, a 3D hydrodynamic model of turbulent GHG transfer in the atmospheric boundary layer has been implemented [15].It is based on the E-ω closure scheme [16][17][18] of the averaged Navier-Stokes equation, the continuity equation and the 'diffusion-advection' equation for the transfer of GHGs [19][20][21][22].In this study we use a modified version of a previously developed 3D hydrodynamic model [15,23] to describe the airflow distribution over the area with complex topography and vegetation.
The following forward (direct) model is used to link surface GHG fluxes to their concentrations in the atmospheric boundary layer: where Π is a rectangle ( ) ( , is the function that describes the lower boundary of the computational domain and characterises the surface topography, n  is an upward unit normal vector to the surface ( ) , H is the upper boundary of the modelling domain, corresponding to the height of the atmospheric boundary layer (ABL), V  is a known steadystate distribution of the mean air flow velocity taking into account the heterogeneities of the underlying surface, K is a turbulent coefficient, Sc is the Schmidt number for corresponding GHG, 0 C corresponds to the mean concentration of GHGs within the ABL, ( ) is the projection of the flux onto the normal vector n  , the terms b F and ph F describe the GHG sources and sinks.
The main parameterisations of b F and ph F functions used to describe the sources and sinks of CO 2 within a plant canopy have been discussed in Mukhartova et al [15].The 'instantaneous' distributions of the mean wind velocity components and the turbulent coefficient were obtained using the relaxation algorithm for the averaged Navier-Stokes equation and the continuity equation using the "1.5-order closure" scheme (Eω model), analogous to the algorithm used in [15]: where f is the Coriolis parameter, g u and g v are the geostrophic wind components, P δ is an overpressure due to the interaction of the air flow with topography and vegetation elements, E is a turbulent kinetic energy (TKE), ( ) , , is a function describing the spatial distribution of the projected plant area density, d c is a dimensionless drag coefficient.
The main difference between this model and the model used in [15] is the presence of the term modelling interaction with topography elements.The topography is modelled as a low permeable obstacle with a high drag coefficient 0 G using the function , where 0 = z corresponds to the grid point with the lowest elevation of the underlying surface.
The turbulence coefficient K is expressed in a 1.5-order closure scheme from the turbulent kinetic energy E and its dissipation rate ε as ε , and the functions E and E ε ω = satisfy the equations of the "diffusion -reaction -advection" type [15][16][17][18], also taking into account the surface topography as an obstacle with a high drag coefficient.
The boundary conditions for the system of equations ( 2) -( 5), as well as the equations for E and ω are described in Mukhartova et al [15].The initial vertical wind distribution is consistent with the logarithmic wind profile and has been modified to account for the surface topography.Let the absolute value of the mean wind speed at the upper boundary of the modelling domain be known and equal to max V .The friction velocity * u can then be approximated by the expression: where κ is the von Kármán constant and 0 z is the surface roughness (for grassy vegetation).The initial wind velocity distribution is described as: For the initial parameterisation of the functions E and ω, we use the classical expressions that are valid over a homogeneous underlying surface, replacing the constant friction velocity * u in these expressions by the function ( 6): Turbulent GHG fluxes can be estimated from a steady-state spatial concentration distribution as follows For advective GHG fluxes we use the parameterization according to [24][25][26]: The use of the full forward model to estimate CO 2 fluxes at some level above the plant canopy is preferable if the rates of soil respiration and canopy CO 2 uptake or release (due to plant photosynthesis and respiration) are known.
In our study, we used the full forward model to find the "theoretical" flux distribution at some level above the canopy and to model the "measured" concentrations at some reference levels in order to verify the representativeness of the inverse modelling results.

Reduced forward problem for GHG concentration calculation
To solve the inverse problem, detailed parameterization of GHG sources and sinks within the canopy (as used in the full forward model) greatly complicates the modelling procedure.As a simplification of the forward model for further implementation in the inverse problem, GHG concentrations can be considered above a given level flux h located higher than the vegetation canopy.The reduced forward problem is: where n  is the unit upward vector normal to the surface is a normal component of known flax at the level flux h .When solving the inverse problem, it's necessary to retrieve the function ( ) y x q n , if the concentration fields at several levels above the level flux h are known.

The inverse problem for flux estimation at a reference height above the ground
If the GHG concentration is known from measurements taken at least at two levels above the plant canopy ( ) ( ) , then the inverse problem of ( ) y x q n , retrieval can be formulated as follows: we have to find a function ( )  (9), where ( ) y x q n , is used as input data.

Numerical solution
Stable implicit finite difference schemes with process splitting are implemented for the numerical solution of the initial boundary problems for the mean wind velocity components, TKE and its dissipation rate, as well as for the GHG concentrations [15].The full forward problem is solved in two steps.In the first step, the "instantaneous" distribution of wind speed and turbulence coefficient is computed, taking into account the heterogeneity of the underlying surface.In the second step, the diffusion-advection equation for the GHG concentrations is solved using the obtained wind speed and K distributions.
In our study, the grid spacing along the x-and y-directions is assumed to be 40 m.In the vertical direction, 524 layers of a quasi-uniform grid with non-linear step growth with height are used.The minimum vertical step is 0.5 m near the surface and the maximum step is 19 m near the top of the computational domain (500 m above the grid point with the lowest elevation of the underlying surface).
The Nelder-Mead multidimensional unconstrained nonlinear minimisation algorithm [27] is used to retrieve ( ) y x q n , at the grid points.

Modelling results
An experimental forest site located at the bottom of the Great Caucasus Mountains (43°03'47.0"N,45°23'13.3"E)near the village of Roshny-Chu in the Chechen Republic, Russia was selected to test and verify the developed model algorithm.The Roshny-Chu area is characterised by a transition type from the low mountain landscapes of the mountain forest zone to the mountain-valley forest-meadowshrub landscape.The experimental plot area is approximately 243 ha and is characterised by a complex topography with elevations ranging from 400 to 600 m (Figure 1).The vegetation is represented by mixed deciduous forests, mainly beech and hornbeam, with admixtures of oak, elm, black alder, alder and maple.In the understorey there are a number of species of wild fruit trees.To parameterise the developed 3D model and solve the full forward problem, experimental data on soil CO 2 efflux and leaf photosynthesis and respiration of different tree species were used.All these parameters were measured during intensive field campaigns in the summer of 2022.Maps of seasonal variation in leaf area index (LAI) were derived using Landsat8 normalised difference vegetation index (NDVI).Numerical experiments were conducted for the prevailing southeast wind direction in the study area for the summer season.The spatial distribution of incoming solar radiation was derived for clear sky and cloudless weather conditions, taking into account the complex surface topography.In the further analysis of the inverse problem, the calculated CO 2 flux at a height of 30 m above the ground surface is considered as "theoretical".The CO 2 concentration fields, which in our case are considered as "measured" at some reference levels above the plant canopy due to the lack of direct UAV measurements, were also determined using the full forward modelling.For this purpose, the CO 2 concentrations are modelled in the nodes of the sparse grid at the reference levels, and possible instrumental errors are added.
Figures 2 and 3 show the spatial distributions of wind speed, K coefficient and theoretical CO 2 flux at 30 m above the ground.To analyse the representativeness of the reduced forward problem, we compared the CO 2 concentration distributions calculated from the model ( 9) using the theoretical flux rate ( ) ( )  above the topographic elevations.The maximum overestimation of the concentration does not exceed 9 ppm.The differences between the concentrations derived by the full forward and the reduced forward models decrease with increasing altitude.However, the sensitivity of the models to the properties of the underlying surface and the accuracy of the lower boundary flux estimation also decrease with altitude.Therefore, the reference heights for UAV concentration measurements should not be positioned too high above the ground surface.For the solution of the inverse problem, the height flux h was chosen to be 30 m, and the heights 1 h ∆ and 2 h ∆ were chosen to be 45 m and 60 m above the ground surface, respectively.
In the inverse problem, the following approximate expression was used as an initial parameterisation of the CO 2 flux: A pre-filtering of the "measured" concentration fields to reduce the error was also performed for the numerical solution.
Figure 6 shows the theoretical distribution of the CO 2 flux, its simplified approximation by equation (11) and the result of solving the inverse problem.

Conclusions
The previously developed three-dimensional hydrodynamic model of the turbulent transfer of GHGs within the atmospheric boundary layer was modified to describe CO 2 fluxes over a complex terrain.We investigated the possible application of a reduced forward model to calculate the GHG concentration above the underlying surface with known CO 2 flux distributions at the upper boundary of the plant canopy.The good qualitative and quantitative agreement of the CO 2 concentrations derived using the reduced forward and full forward models in the case of known fluxes was shown.The inverse problem of retrieving the GHG flux from the known concentration distributions at two reference levels above the ground surface was performed.The CO 2 fluxes retrieved by solving the inverse problem show a good quantitative agreement with the theoretical one and reflect a strong complexity of the spatial flux distribution influenced by surface topography and wind direction.Thus, the numerical experiments carried out showed the ability of the developed approach to provide an adequate and representative flux estimation.However, the effects of measurement errors, the heights of selected reference levels above the ground surface, as well as wind direction on the accuracy of flux estimation are still open and require further multifaceted theoretical and experimental studies.

Figure 1 .
Figure 1.Surface topography (a) and satellite image (b) of the Roshny-Chu site.

Figure 2 .
Figure 2. Spatial distribution of horizontal (a) and vertical (b) components of wind speed at 30m above ground level.

Figure 3 .
Figure 3. Spatial distribution of the turbulence coefficient K (a) and the normal part n q of the total CO 2 flux (b) at 30 m above the ground level.
flux h as input data with the corresponding CO 2 concentration distributions obtained from simulations of the full forward model (1).The results of this comparison are shown in Figures 4 and 5.

Figure 4 .Figure 5 .
Figure 4. Distributions of CO 2 concentrations in a cross-section 2 y L y = : full forward model (a); reduced forward model (b).A thin solid line corresponds to the level flux h , a dashed line corresponds to the level

Figure 6 .
Figure 6.Theoretical CO 2 flux distribution (a) and initial CO 2 flux approximation (b), as well as the CO 2 fluxes obtained by solving the inverse problem (c).