remote sensing Impact of Drought on Isoprene Fluxes Assessed Using Field Data, Satellite-Based GLEAM Soil Moisture and HCHO Observations from OMI

: Biogenic volatile organic compounds (BVOCs), primarily emitted by terrestrial vegetation, are highly reactive and have large effects on the oxidizing potential of the troposphere, air quality and climate. In terms of global emissions, isoprene is the most important BVOC. Droughts bring about changes in the surface emission of biogenic hydrocarbons mainly because plants suffer water stress. Past studies report that the current parameterization in the state-of-the-art Model of Emissions of Gases and Aerosols from Nature (MEGAN) v2.1, which is a function of the soil water content and the permanent wilting point, fails at representing the strong reduction in isoprene emissions observed in ﬁeld measurements conducted during a severe drought. Since the current algorithm was originally developed based on potted plants, in this study, we update the parameterization in the light of recent ecosystem-scale measurements of isoprene conducted during natural droughts in the central U.S. at the Missouri Ozarks AmeriFlux (MOFLUX) site. The updated parameterization results in stronger reductions in isoprene emissions. Evaluation using satellite formaldehyde (HCHO), a proxy for BVOC emissions, and a chemical-transport model, shows that the adjusted parameterization provides a better agreement between the modelled and observed HCHO temporal variability at local and regional scales in 2011–2012, even if it worsens the model agreement in a global, long-term evaluation. We discuss the limitations of the current parameterization, a function of highly uncertain soil properties such as porosity.


Introduction
Biogenic volatile organic compounds (BVOCs) are highly reactive and have large effects on the oxidizing potential of the troposphere, air quality and climate [1][2][3]. They are mainly emitted by terrestrial vegetation and amount to ca. 1000 Tg year −1 globally, with about half of the share consisting of one compound, isoprene [4][5][6][7][8]. Isoprene influences tropospheric hydroxyl radical levels through multiple processes [9][10][11][12][13][14], thereby altering the lifetime of methane. Under high-NO x conditions, it is a major precursor of tropospheric ozone [2,[15][16][17][18] and it also affects the growth of secondary organic aerosols [19][20][21][22], both of which are considered as important short-lived climate forcers. Extreme events such as droughts bring about changes in the surface emissions of biogenic hydrocarbons. Under mild to moderate drought, water stress can either have no significant impact on isoprene emission or can stimulate it, even in case of a significant inhibition of photosynthesis (e.g., [23][24][25]). In fact, there is evidence that isoprene improves plant thermotolerance [26][27][28]. This increase in isoprene emissions has been hypothesized as a possible result of the drought-induced stomatal closure, which causes a reduction in transpiration rates as the plants attempt to conserve water, and subsequently increases leaf temperature [29,30]. However, as water stress becomes more severe, the general consensus is that isoprene emission capacity of a plant is greatly reduced [23][24][25][30][31][32].
The effects of drought stress on isoprene emissions have been investigated in several studies, most of which were conducted in the laboratory [33][34][35][36][37][38]. Field observations at ecosystem-scale during natural droughts are scarce [30,32,39,40] and, to the authors' knowledge, only one site provides measurements of isoprene fluxes. During the 2011-2012 field campaigns in central U.S., at the Missouri Ozarks AmeriFlux (MOFLUX) site, isoprene and other BVOC fluxes were monitored for two consecutive growing seasons concomitant to two drought events [30,32]. Moreover, the effect of drought on isoprene may also be monitored at greater spatiotemporal scales using other proxies. Formaldehyde (HCHO) is a key intermediate in the oxidation of the large majority of hydrocarbons released by vegetation, anthropogenic activities and fire. BVOCs, among which isoprene is largely dominant, are the major precursors of HCHO over vegetated areas (e.g., [41][42][43]), and therefore, satellite HCHO is an excellent proxy for the emissions of BVOCs. The anthropogenic signal of HCHO is low and only detectable in highly polluted regions [44,45], and does not play a role in the isoprene-dominated MOFLUX site area. The contribution of fires to the HCHO abundances is generally important in the vicinity of fires or downwind, but only small agricultural fires occurred during the summers of 2011 and 2012 over Missouri. The BVOC emissions are therefore the main driver of the HCHO interannual variability over Eastern US [46]. Zheng et al. (2017) [47] have shown that remotely sensed sun-induced fluorescence (SIF) and HCHO columns qualitatively capture reductions in flux-derived gross primary productivity (GPP) and isoprene emission at the MOFLUX site, respectively, on weekly to monthly time scales, but with muted responses. They concluded that even fairly moderate reductions in satellite SIF and HCHO column could be indicative of severe soil drought.
The present study aims at assessing the response of isoprene emissions to drought modelled by the Model of Emissions of Gases and Aerosols from Nature (MEGAN; [5,48]). The current parameterization in MEGANv2.1, dependent on root-zone soil water content and permanent wilting point (denoted θ W ), was derived based on a study conducted on potted plants [33]. It was found unable to account for the response of isoprene emissions to the severe drought during the 2012 MOFLUX campaign, and the low value of θ W assumed for the site was designated the main culprit [32]. In addition to errors in the formulation and θ W , uncertainties have been pointed out in relation to the MEGANv2.1 drought parameterization. The reduction in global annual isoprene emissions entailed by this drought factor ranges between 7 and 50% for different model configurations [6,[48][49][50][51]. The range is even wider (10-70%) when varying the soil moisture input datasets ( [52,53]). In addition, MEGAN is typically run using reanalysis data, and thus does not benefit from the four-decades of progress in the field of microwave surface soil moisture retrievals [54].
The objective of this study is to update the drought parameterization in MEGANv2.1 [5] using (i) the 2011-2012 MOFLUX field measurements, and (ii) remotely sensing-based rootzone soil moisture from the Global Land Evaporation Amsterdam Model (GLEAM; [55,56]). GLEAM assimilates surface soil moisture retrievals from the ESA Climate Change Initiative soil moisture (ESA CCI SM; [57]) into a multilayer soil model driven by observationbased precipitation. GLEAM soil moisture recently proved to outperform the reanalysis estimates used in MEGAN in comparison to in situ measurements [58]. The influence on isoprene emissions of the updated MEGAN parameterization using GLEAM data is evaluated at different scales using the models (biogenic emission and chemistry-transport) Remote Sens. 2022, 14, 2021 3 of 23 and observations (in situ and remotely sensed) described in Section 2. The assessment of the re-adjusted algorithm at the site-, regional-and global-scales is presented in Section 3, with the regional study performed over the contiguous U.S. for 2011-2012. The evaluation method used here is based on spaceborne HCHO and a chemistry-transport model, and relies on the high HCHO yield from the oxidation of isoprene. This method has been widely used in past studies [42,43,[59][60][61][62]. Discussion and conclusions are drawn in Section 4.

Site Description and Eddy-Covariance Flux Measurements
The MOFLUX field site is located in central U.S., in the Ozarks region of Missouri state (38 • 44.65 N latitude, 92 • 12.00 W longitude, 219 m elevation; Figure 1). It is part of the AmeriFlux network and has provided continuous eddy covariance measurements of carbon dioxide, water and energy fluxes since 2004. The 32-m tower is surrounded by a continuous broadleaf deciduous forest within a 1-km radius. The tower extends about 10 m above the top of the canopy. The forest is oak dominated, with two thirds of the trees belonging to white and red oak species, and the remaining third being hickory, sugar maple, ash, and juniper. With this abundance of isoprene-emitting species of oaks, the Ozarks area has gained the nickname of 'isoprene-volcano', holding the record for the highest isoprene fluxes ever reported (53.3 mg m −2 h −1 ) during the 2011 campaign [30]. The relatively thin silt loam soils prevent the trees from reaching deep soil moisture during droughts, which exacerbates plant water stress [63]. Soil water retention curves ( [64]) for the Weller silt loam soils near the MOFLUX eddy covariance tower indicate a wilting point of 0.13 m 3 m −3 for the top 0-30 cm soil depth, and 0.23 m 3 m −3 for soil depths below 30 cm. These figures correspond to the water content at −1.5 MPa soil matric potential, a value that commonly defines θ W [65]. nalysis estimates used in MEGAN in comparison to in situ measurements [58]. The influence on isoprene emissions of the updated MEGAN parameterization using GLEAM data is evaluated at different scales using the models (biogenic emission and chemistrytransport) and observations (in situ and remotely sensed) described in Section 2. The assessment of the re-adjusted algorithm at the site-, regional-and global-scales is presented in Section 3, with the regional study performed over the contiguous U.S. for 2011-2012. The evaluation method used here is based on spaceborne HCHO and a chemistrytransport model, and relies on the high HCHO yield from the oxidation of isoprene. This method has been widely used in past studies [42,43,[59][60][61][62]. Discussion and conclusions are drawn in Section 4.

Site Description and Eddy-Covariance Flux Measurements
The MOFLUX field site is located in central U.S., in the Ozarks region of Missouri state (38°44.65′N latitude, 92°12.00′W longitude, 219 m elevation; Figure 1). It is part of the AmeriFlux network and has provided continuous eddy covariance measurements of carbon dioxide, water and energy fluxes since 2004. The 32-m tower is surrounded by a continuous broadleaf deciduous forest within a 1-km radius. The tower extends about 10 m above the top of the canopy. The forest is oak dominated, with two thirds of the trees belonging to white and red oak species, and the remaining third being hickory, sugar maple, ash, and juniper. With this abundance of isoprene-emitting species of oaks, the Ozarks area has gained the nickname of 'isoprene-volcano', holding the record for the highest isoprene fluxes ever reported (53.3 mg m −2 h −1 ) during the 2011 campaign [30]. The relatively thin silt loam soils prevent the trees from reaching deep soil moisture during droughts, which exacerbates plant water stress [63]. Soil water retention curves ( [64]) for the Weller silt loam soils near the MOFLUX eddy covariance tower indicate a wilting point of 0.13 m 3 m −3 for the top 0-30 cm soil depth, and 0.23 m 3 m −3 for soil depths below 30 cm. These figures correspond to the water content at −1.5 MPa soil matric potential, a value that commonly defines [65]. The area is classified as a humid subtropical climate in the Köppen classification and is prone to droughts between July and September [66,67]. In 2011 and 2012, the site has respectively suffered from mild and severe droughts for which two campaigns of VOC The area is classified as a humid subtropical climate in the Köppen classification and is prone to droughts between July and September [66,67]. In 2011 and 2012, the site has respectively suffered from mild and severe droughts for which two campaigns of VOC flux measurements were conducted. [30] used a chemiluminescence analyzer for the 2011 campaign spanning from 18 May to 2 September. The dataset is provided with a half-hour frequency. The 2012 campaign of [32] was extended to other compounds than isoprene. Eddy covariance fluxes were calculated half-hourly with a proton transfer reaction quadrupole mass spectrometer (PTR-Quad-MS) and covered a nearly six-months period from 2 May to 22 October. According to the U.S. Drought Monitor, the most intense period Remote Sens. 2022, 14, 2021 4 of 23 of drought occurred in August 2012 with exceptional drought conditions (i.e., level D4) affecting more than a third of Missouri land (https://www.drought.gov/states/missouri, accessed on 23 June 2021).

MEGAN-MOHYCAN: Biogenic VOC Emission Modelling
The biogenic emissions of isoprene, monoterpenes, and 2-methylbutenol are estimated using MEGAN [5,48] coupled to the MOdel for HYdrocarbon emissions by the CANopy (MOHYCAN; [49,68]), a multi-layer canopy environment model. The net emission rates F (µg m −2 h −1 ) into the above-canopy atmosphere are estimated as follows: where the emission factor ε (µg m −2 h −1 ) represents the emission rate at standard conditions as defined by [48]. Deviations from those conditions are accounted for by the activity factors encapsulated in γ 0 representing the response of biogenic emissions to their major identified environmental and phenological drivers. Activity factors account for leaf temperature and photosynthetic photon flux density (γ PT ), leaf area index (LAI), leaf age (γ A ), drought stress (γ SM ) and CO 2 inhibition (γ CO 2 ). The temperature and light response algorithm incorporates the influence of current and past conditions, as well as the attenuation of light by leaves within the canopy. The calculation of γ PT is performed by MOHYCAN, the canopy environment model, detailed by [49]. The adjustment factor C CE is specific for each canopy environment model and is set to 0.52 for MOHYCAN, so that γ 0 = 1 for standard canopy conditions. The response of biogenic isoprene emissions to drought in the MEGAN model is accounted for primarily in the γ SM activity factor through the direct impact of soil water stress on the emissions. Moreover, drought leads to enhanced stomatal resistances and therefore higher leaf temperatures accounted for in γ PT defined in MOHYCAN. For severe drought conditions, however, the impact of γ SM outweighs the impact of higher leaf temperatures, and can even lead to a complete shutdown of the emissions [33].
Currently, the MEGAN model provides the users with two algorithms of γ SM for the response of emissions to drought. In MEGANv2.1, the impact of drought is estimated based on the volumetric soil water content (denoted θ) and θ W . The factor γ SM accounts for direct deviations due to changes in θ as follows: The original value of ∆θ (0.06 m 3 m −3 ; [48]) was derived empirically based on the study of [33] conducted on potted plants, with one plant species (Populus deltoides Bartr.) and soil type. The value ∆θ was later lowered to a more conservative 0.04 m 3 m −3 [5], because emissions were being shut off in regions where there was no extreme drought, and will further be referred to as the default value for this parameter. The threshold point θ 1 is defined as θ W + ∆θ. As emphasized by [49], the soil water content should be paired consistently with the wilting point given the importance of the wilting point in the determination of soil moisture.
In the third version of the model (MEGANv3; [69]), isoprene emissions respond simultaneously to water stress and photosynthesis, as provided by the Community Land Model (CLM; [70]) in the Community Earth System Model (CESM; [71]), which is coupled to MEGAN. The parameterization of the drought stress activity factor γ SM in MEGANv3 (hereafter denoted γ * SM ) is based on CLM maximum rate of carboxylation (V CMAX ) and soil moisture stress function (β t ). A thorough description of these parameters is provided in [70] and the parameterization is fully discussed in [69]. The estimation of γ * SM follows: The empirical parameter α (=37) was calibrated by [69] based on field measurements at the MOFLUX site.
Two evaluation approaches are considered to assess the response of biogenic isoprene emissions to drought in the MEGAN biogenic emission model. Single-point estimations of isoprene fluxes using algorithms from Equations (2) and (3) are evaluated against measurements conducted in 2011 and 2012 presented hereinabove. Global and regional (over North America) approaches are also conducted through the assessment of tropospheric HCHO columns using satellite-based observations against simulations from a chemistry transport model that accounts for γ SM (Equation (2)) with different values of ∆θ.
The MEGAN-MOHYCAN model is set up for the MOFLUX site as follows. Singlepoint model simulations of daily isoprene fluxes for years 2011 and 2012 are performed using local meteorological conditions provided by the AmeriFlux website (https://ameriflux. lbl.gov/sites/siteinfo/US-MOz, accessed on 23 June 2021). Hourly measurements of air temperature and photosynthetic photon flux density from the tower are used. Daily LAI measurements for the 2012 campaign were obtained from [32]. Since direct local measurements of LAI were not available to us for 2011, we use the daily 2012 values, corrected for interannual variation based on monthly averaged Moderate Resolution Imaging Spectroradiometer (MODIS) LAI (MODIS 15A2H collection 6 available at https://lpdaac.usgs.gov, accessed on 15 June 2021) at 0.5 • resolution in the vicinity of the MOFLUX site. More precisely, the 2012 local values are multiplied by the ratio of MODIS LAI between the two years ( Figure S1). The canopy emission factor ε was set to 8.6 mg m −2 h −1 . The value is chosen such that the model without any drought stress correction (γ SM = 1) matches the average measured isoprene fluxes of the 2011 campaign, when the effect of drought was of lesser magnitude. As will be seen in Section 3, the γ SM parameterization has an almost negligible impact in 2011.
In the following, the simulation without stress (γ SM = 1) will be referred to as the baseline. The impact of drought stress is accounted for by multiplying the baseline emissions by either daily γ SM (Equation (2)) or γ * SM (Equation (3)). The latter was obtained using CLM5 (instead of CLM4.5 as in Jiang et al., 2018) and local environmental conditions observed at the MOFLUX site. The daily γ SM is obtained using daily satellite-based rootzone θ provided by the GLEAM model at 0.25 • × 0.25 • spatial resolution, along with the θ W distribution taken from the static soil properties defined in GLEAM (see Section 2.3). For reasons explained in Section 3.1, the empirical parameter ∆θ is re-evaluated and set to ∆θ = 0.12; the adjusted factor using GLEAM is referred to as γ OPT SM . The drought factor using the default ∆θ = 0.04 is denoted γ 0.04 SM .

Satellite-Based Soil Properties from GLEAM
GLEAM [55,56] is a set of algorithms to derive the different components of terrestrial evaporation and soil moisture (surface and root-zone) from satellite-based datasets of different hydro-climatic variables at a daily timescale. It employs the Priestley and Taylor equation [72] to estimate potential evaporation for different land cover types (short and tall vegetation, bare soil, ice and snow covers, and open water). A semi-empirical, land cover-dependent multiplicative stress factor is used to translate the potential evaporation estimates into actual evaporation. For tall and short vegetation, the stress factor is based on vegetation optical depth [73] and root-zone soil moisture. For bare soil, the stress factor is dependent on surface soil moisture only. Interception loss is estimated using the Gash analytical model [74]. To model surface and root-zone soil moisture, GLEAM applies a multi-layer soil water balance. The number of soil layers and the depth of the root-zone Remote Sens. 2022, 14, 2021 6 of 23 depend on land-cover type with a maximum depth of 10 cm for bare soil, 100 cm for low vegetation, and 250 cm for tall vegetation, respectively. The soil water available for evaporation, either via transpiration or bare soil evaporation, is estimated using static soil properties (e.g., θ W ) taken from the database of Global Gridded Surfaces of Selected Soil Characteristics generated by the International Geosphere-Biosphere Programme Data and Information System (IGBP-DIS, [75]). Finally, GLEAM incorporates a data assimilation module which ingests soil moisture observations via a Newtonian nudging scheme. In the version of GLEAM used in this study (GLEAMv3.5b), the latest satellite-based estimates of surface soil moisture from ESA-CCI [76] are assimilated, which include observations from the Soil Moisture Active Passive (SMAP) satellite. GLEAMv3.5b data at daily and 0.25 • × 0.25 • resolution are available for download at www.gleam.eu (accessed on 1 April 2022).

Remotely Sensed Formaldehyde Columns from OMI
We use formaldehyde observations retrieved from the Ozone Monitoring Instrument (OMI) spectrometer aboard the Aura mission launched in 2004 [77]. The OMI HCHO product used in this study was developed in the framework of the EU FP7 project QA4ECV (Quality Assurance for Essential Climate Variables, http://www.qa4ecv.eu, accessed on 1 June 2018) and is fully documented by [78,79]. The retrievals of HCHO slant columns are based on improved Differential Optical Absorption Spectroscopy (DOAS) algorithms that reduce the effect of interferences between species using the 328.5-359 nm spectral window with the HCHO absorption cross-section from [80]. Further, the background normalization of the slant columns is performed using the equatorial Pacific Ocean in order to compensate for possible systematic latitude-dependent offsets in spectral fitting. The conversion into a vertical column is performed with the air mass factor (AMF) obtained by combining an altitude-resolved AMF look-up table derived at 340 nm from the VLIDORTv2.6 radiative transfer model [81] and the daily a priori formaldehyde vertical profile shape calculated with the TM5-MP chemical transport model [82,83]. The scattering due to clouds is corrected using an independent pixel approximation [84], whereas the effects of non-absorbing aerosols are implicitly accounted for through the cloud correction scheme. The effect of absorbing aerosols is currently not accounted for in the retrieval algorithm. Those effects can be significant over regions dominated by biomass burning or over heavily industrialized regions.
The OMI QA4ECV dataset provides separate estimates for the systematic and random components of the total uncertainty and a thorough description of the error estimation is provided in [85]. The total uncertainty on the monthly HCHO vertical column generally ranges between 20% and 40%, depending on the observation conditions. For clear sky conditions, in absence of aerosols, the total error on the AMFs is estimated at 18% and the error increases with cloud fraction, up to 50% for a cloud fraction of 0.5. The omission of the aerosol correction may lead to a significant underestimation of the derived HCHO column by up to 40% over fire scenes [80]. To ensure the best quality, we filter out cloudy pixels (fractions higher than 20%) in addition to the quality assurance criterion (http://www. tropomi.eu/sites/default/files/files/publicSentinel-5P-Formaldehyde-Readme.pdf, accessed on 25 September 2021). This already limits the potential impact of thick plumes on the retrievals. Furthermore, we conducted an additional sensitivity analysis aiming to minimize the potential impact of aerosols on the model comparisons with OMI HCHO data. This was done by filtering out OMI pixels with high aerosol optical depth (AOD), as described in Section 2.6 below and Section S6 of the Supplementary Materials.
For the comparison with MOFLUX data, the OMI observations are averaged within a 20-km radius around the MOFLUX site. For the regional and global evaluation, observations at 0.5 • × 0.5 • and 2 • × 2.5 • are used, respectively.

Simulated HCHO Columns Using MAGRITTE and IMAGES CTMs
Simulations of HCHO columns are performed using the Model of Atmospheric composition at Global and Regional scales using Inversion Techniques for Trace gas Emissions (MAGRITTE) v1.1, a global and regional three-dimensional chemistry-transport model (CTM) of the troposphere [86,87]. Most model parameterizations are inherited from the global Intermediate Model of the Global and Annual Evolution of Species (IMAGESv2) model [46,61,[88][89][90]. MAGRITTE calculates the concentrations of 182 chemical compounds in the troposphere with a spin-up time of six months. The horizontal resolution of the model is either 2 • (latitude) × 2.5 • (longitude) at a global scale or 0.5 • × 0.5 • at a regional scale, while the vertical coordinate is a hybrid sigma-pressure system resolved in 40 unevenly spaced levels extending from the surface to the lower stratosphere (44 hPa). Meteorological fields are provided by European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalyses. The model calculates monthly averaged columns at the OMI overpass time (13:30), taking into account the sampling date statistics of the satellite. The chemical degradation mechanism of isoprene relies on the Leuven isoprene mechanism [11,91] and the Caltech oxidation mechanism [92], accounting for relevant theoretical and laboratory studies.
The bottom-up fluxes of HCHO precursors are prescribed as follows. The biomass burning inventory is obtained from the Global Fire Emission Database version 4 (GFED4s; [93]) on a daily basis (https://globalfiredata.org, accessed on 22 March 2013). The anthropogenic sources of non-methane VOC species are taken from the Emission Database for Global Atmospheric Research version 4.3.2 (EDGARv4.3.2, [94]) and the anthropogenic CO, NO x , OC, SO 2 , BC and NH 3 are obtained from the Hemispheric Transport of Air Pollution version 2 (HTAPv2; [95]) database for 2010. Due to their limited availability, the EDGARv4.3.2 emissions are set constant at their 2012 values after this year. Both inventories are available at https://edgar.jrc.ec.europa.eu (accessed on 28 August 2018). Over the U.S., the anthropogenic EDGAR emissions are scaled down to match emission estimates based on aircraft NO x concentrations and nitric acid deposition data, according to [96]. Biogenic emissions are prescribed as follows. As in [46], MEGAN-MOHYCAN global simulations provided fluxes of isoprene, monoterpene, and 2-methyl-3-buten-2-ol (MBO) using meteorological fields from ECMWF reanalysis (ERA-Interim; [97]) and gridded emission factors [5] based on detailed land cover description and species-specific emission factors. The global bottom-up inventory of biogenic isoprene emissions is publicly available from emissions.aeronomie.be/index.php/bottom-up/isoprenev2018 (accessed on 1 April 2022). Biogenic sources of ethylene, formaldehyde, and acetone, based on the MEGAN model are taken from the database of Emissions of atmospheric Compounds and Compilation of Ancillary Data (ECCAD, http://eccad.aeris-data.fr, accessed on 3 June 2019). Biogenic emissions of acetaldehyde and ethanol are parameterized as in [98]. Finally, the biogenic methanol emissions are provided by inverse modelling constrained by spaceborne methanol data [99]. Finally, monthly 3-dimensional methane concentration distributions are prescribed, based on an assimilation of surface data at 0.75 • resolution provided by the Copernicus Atmospheric Monitoring System (CAMS) (https://ads.atmosphere.copernicus.eu, accessed on 21 April 2022).
The HCHO column assessment is performed at three different scales: local, regional and global. Regional simulations of formaldehyde columns are conducted using the regional MAGRITTE model at a 0.5 • resolution for the U.S. (10-54 • N, 65-130 • W) for a period of 30 months starting on 1 July 2010 (including the spin-up time). For the local analysis, we extract HCHO columns at the 0.5 • grid cell of the regional model including the MOFLUX site. The lateral boundary conditions of the regional model are provided by a global model simulation using the same chemical mechanism and emission datasets. Three model simulations were performed: one simulation without water stress (referred to as MAGRITTE, γ SM = 1), and two simulations with the monthly γ SM factor calculated using Equation (2) with monthly root-zone θ from GLEAM with different values of the empirical parameter ∆θ, namely the default value 0.04 (referred to as MAGRITTE, γ 0.04 SM ) and the value of 0.12 obtained from the optimization of the model at the MOFLUX site (referred to as MAGRITTE, γ OPT SM ). The global, long term comparison relies on the monthly HCHO columns obtained with IMAGES, for which two simulations are performed over the 2005-2016 period at a horizontal resolution of 2 • × 2.5 • , using either the default value 0.04 (referred to as IMAGES, γ 0.04 SM ) or the value of 0.12 (referred to as IMAGES, γ OPT SM ).

Aerosol Optical Depth Dataset
For our alternative analysis based on OMI data excluding high-AOD pixels, we use gridded, total AOD at 469 nm from the CAMS global reanalysis EAC4 (ECMWF Atmospheric Composition Reanalysis 4; https://ads.atmosphere.copernicus.eu/#!/home, accessed on 22 February 2022). CAMS AOD is constrained by MODIS [100] and the Advanced Along-Track Scanning Radiometer (AATSR; [101]) aboard Envisat from 2003 to March 2012. Moreover, they are routinely evaluated against in situ data such as the AErosol RObotic NETwork (AERONET). A more complete description of the CAMS aerosol scheme can be found in [102].

Soil Water Content and Soil Moisture Stress Factor
The γ SM factor (Equation (2)) depends explicitly on the plant water status. Figure 2 shows the daily θ as measured at a depth of 100 cm and as represented by the root-zone (0-250 cm) GLEAM product at the grid cell encompassing the MOFLUX site in years 2011 and 2012. The exceptional conditions of 2012 are shown by the decrease in the observed θ, from 67% at the beginning of May to a minimum of 34% by the end of August. The drought condition was relieved with the return of precipitation illustrated by the sharp increase in the θ at the end of August 2012. The mild drought of 2011 had a lesser impact on θ, which decreased from 70% in July to 56% in August.
to as MAGRITTE, = 1), and two simulations with the monthly factor calculated using Equation (2) with monthly root-zone from GLEAM with different values of the empirical parameter Δ , namely the default value 0.04 (referred to as MAGRITTE, 0.04 ) and the value of 0.12 obtained from the optimization of the model at the MOFLUX site (referred to as MAGRITTE, ). The global, long term comparison relies on the monthly HCHO columns obtained with IMAGES, for which two simulations are performed over the 2005-2016 period at a horizontal resolution of 2° × 2.5°, using either the default value 0.04 (referred to as IMAGES, 0.04 ) or the value of 0.12 (referred to as IMAGES, ).

Aerosol Optical Depth Dataset
For our alternative analysis based on OMI data excluding high-AOD pixels, we use gridded, total AOD at 469 nm from the CAMS global reanalysis EAC4 (ECMWF Atmospheric Composition Reanalysis 4; https://ads.atmosphere.copernicus.eu/#!/home, accessed on 22 February 2022). CAMS AOD is constrained by MODIS [100] and the Advanced Along-Track Scanning Radiometer (AATSR; [101]) aboard Envisat from 2003 to March 2012. Moreover, they are routinely evaluated against in situ data such as the AErosol RObotic NETwork (AERONET). A more complete description of the CAMS aerosol scheme can be found in [102].

Soil Water Content and Soil Moisture Stress Factor
The factor (Equation (2)) depends explicitly on the plant water status. Figure    Following the γ SM parameterization of MEGANv2.1, both the locally measured θ (orange line on Figure 2) and the GLEAM θ (blue line) always exceed the default threshold value, θ 1 = θ W + 0.04, represented by dashed lines in Figure 2. The threshold is estimated at 0.27 m 3 m −3 based on local θ W , and at 0.21 m 3 m −3 in GLEAM based on the IGBP-DIS θ W value. This implies that the isoprene fluxes are not impacted by the drought stress activity factor (i.e., γ SM = 1) according to the MEGANv2.1 parameterization.
Therefore, the current parameterization fails to capture the drought effect due to the low value of the empirical parameter ∆θ, which was based on only one study [33]. Therefore the ∆θ value is re-calibrated using the MOFLUX measurements, in order to derive an optimal value of the parameter ∆θ that provides the best match with the MOFLUX isoprene dataset. To that aim, we vary the ∆θ parameter in the model and select the value which minimizes the root-mean square error (RMSE) between modelled and in situ daily fluxes in both years. Using the soil moisture and θ W from GLEAM, the optimized ∆θ is 0.12. Using this value, the GLEAM θ falls below the threshold θ 1 of 0.29 (= 0.12 + 0.17) during several weeks in 2011 (after 28 July) and the entire summer and early fall in 2012 (after 5 June).

Daily Time Series and Diel Cycle
During the mild drought year, the daily drought stress factors have virtually no impact on isoprene fluxes except when adopting the optimized factor γ OPT SM , which entails a flux reduction of about 10% during the month of August (Figure 3a). The bias-corrected MEGAN simulations explain 93% of the variability of daily isoprene fluxes during summer 2011 for the simulations using γ OPT SM or using γ SM = 1 (Figure 3b). Measurements are missing during 10-17 August due to an equipment failure. The reasons for the model overestimation in the second half of August are unclear, but could be partly due to an underestimation of the effect of leaf age as previously noted in MEGAN comparisons with flux measurements at a temperate forest [49].
In the severe stress year (2012), the model performance is lowest in the no-stress configuration (γ SM = 1) both in terms of correlation (r = 0.76) and root-mean square error (RMSE =~5 mg m −2 h −1 , Figure 3d). The γ * SM factor improves only slightly the comparison (r = 0.79 and RMSE =~3 mg m −2 h −1 ), whereas the application of γ OPT SM improves the fit in 2012 to r = 0.90 and RMSE =~2 mg m −2 h −1 . The stress factor γ OPT SM decreases the isoprene fluxes from early June to late October, with a maximum reduction of about 70% at the end of August, and a significant decrease of 40-50% in late summer (Figure 3c). In the measurement dataset, the highest daily flux occurs on June 24, at~18 mg m −2 h −1 . Without the stress factor, modelled isoprene fluxes peak 12 days later, on 6 July. Those fluxes overestimate the observations by 6.00 mg m −2 h −1 in July August (a percentage bias of about 98%). This bias is strongly reduced (to 0.42 mg m −2 h −1 ) due to the application of γ OPT SM i.e.,~7%. The drop in emissions is sharper with γ * SM , up to 90% in July August (Figure 3c), leading to a model underestimation reaching 3.13 mg m −2 h −1 (about −58%) for that same period (Figure 3d). The improvement of estimated isoprene fluxes using γ OPT SM is also seen in the drop in the RMSE from~5 mg m −2 h −1 to roughly 2 mg m −2 h −1 . Overall, the comparison shown in Figure 3 points to a better agreement with the use of γ OPT SM compared to γ * SM . Finally, the model succeeds at reproducing the normalized diel pattern of emission fluxes (Figure 4), independent of the application of the drought stress factor. Nevertheless, the diurnal peak occurs sooner in 2011, by one hour in June, and as much as 2 h in late summer. Remote Sens. 2022, 14, x FOR PEER REVIEW 10 of 23  [30] and in 2012 by [32].  [30] and in 2012 by [32]. Remote Sens. 2022, 14, x FOR PEER REVIEW 11 of 23

Local and Regional Simulations of Formaldehyde Columns
Here we examine the impact of drought stress on the temporal variability of isoprene fluxes and HCHO columns at the MOFLUX site. Both observed and calculated data are averaged biweekly (14 days) ( Figure 5) and monthly ( Table 1). As seen in Figure 3, isoprene fluxes exhibit great temporal variability at short timescales, but its evaluation using OMI HCHO observations cannot be performed due to the high noise and limited set of valid data (Section S2 of the Supplementary Materials).
Observed biweekly OMI HCHO columns averaged over a 20-km radius around the site (Figure 5b) follow essentially the same pattern as in situ measured isoprene fluxes (Figure 5a), although with a dampened amplitude. This dampening is expected due to the presence of other HCHO precursors less affected by drought (e.g., methane, methanol and anthropogenic VOCs) and due to complex transport and chemical effects, as discussed for instance by [41]. MAGRITTE simulations show the same biweekly variability during summer 2011, with HCHO columns reduced when using the drought factor . Whereas the observed columns are relatively constant in July and early August, MAGRITTE simulations exhibit a greater biweekly variability with a mid-July peak, reflecting the temporal variability of isoprene fluxes. In 2012, OMI observations display a strong decline of HCHO column following the peak in early July, which is not reproduced by MAGRITTE simulations. Overall, application of the drought factor reduces the positive biases in biweekly averages of HCHO columns estimated in the no-stress simulation as shown in Figure 5b. Note however that the OMI systematic errors are important, of the order of 50%, and encompass both MAGRITTE simulations, except in early August 2012.

Local and Regional Simulations of Formaldehyde Columns
Here we examine the impact of drought stress on the temporal variability of isoprene fluxes and HCHO columns at the MOFLUX site. Both observed and calculated data are averaged biweekly (14 days) ( Figure 5) and monthly ( Table 1). As seen in Figure 3, isoprene fluxes exhibit great temporal variability at short timescales, but its evaluation using OMI HCHO observations cannot be performed due to the high noise and limited set of valid data (Section S2 of the Supplementary Materials).
Observed biweekly OMI HCHO columns averaged over a 20-km radius around the site (Figure 5b) follow essentially the same pattern as in situ measured isoprene fluxes (Figure 5a), although with a dampened amplitude. This dampening is expected due to the presence of other HCHO precursors less affected by drought (e.g., methane, methanol and anthropogenic VOCs) and due to complex transport and chemical effects, as discussed for instance by [41]. MAGRITTE simulations show the same biweekly variability during summer 2011, with HCHO columns reduced when using the drought factor γ OPT SM . Whereas the observed columns are relatively constant in July and early August, MAGRITTE simulations exhibit a greater biweekly variability with a mid-July peak, reflecting the temporal variability of isoprene fluxes. In 2012, OMI observations display a strong decline of HCHO column following the peak in early July, which is not reproduced by MAGRITTE simulations. Overall, application of the drought factor γ OPT SM reduces the positive biases in biweekly averages of HCHO columns estimated in the no-stress simulation as shown in Figure 5b. Note however that the OMI systematic errors are important, of the order of 50%, and encompass both MAGRITTE simulations, except in early August 2012. Monthly isoprene fluxes were enhanced in June 2012 compared to the previous year ( Figure 5a and Table 1). More favourable (i.e., warm and sunny) environmental conditions at the onset of the severe drought explain this increase, as shown by the model ability to match the large emission enhancement in June (+47% in the observations and 34-60% in the model). In subsequent months, observed isoprene fluxes were about 40% and 50% lower than in the previous year, representing drastic decreases relative to the model expectation ignoring drought stress effects (+16% and −14% in July and August, respectively). The updated MEGAN model based on reproduces fairly well the seasonal variation (Figure 5a) as well as the inter-annual variability of summer monthly averages (Table 1).
A peak was observed in OMI HCHO in July 2011 at about 22 × 10 15 molecules cm −2 mainly due to high temperatures entailing fast VOC oxidation and strong biogenic emissions. About half this column was observed in August 2012 (12.5 × 10 15 molecules cm −2 ) when severe drought took its toll on biogenic emissions. While the small increase observed in June 2012 with respect to the previous year (+7%) is well simulated in all runs, the substantial drop in July (−19%) is very well represented in MAGRITTE, (−20%), whereas MAGRITTE, = 1 and MAGRITTE, 0.04 barely account for half of that reduction (Table 1). Yet in August, the relative impact of drought in MAGRITTE, is slightly overestimated (−28% vs. −23% in OMI). While those comparisons point to a better performance of MAGRITTE, at reproducing the interannual variability of monthly averaged HCHO columns, uncertainties in OMI relative changes (given between brackets in Table 1) do not allow for a robust conclusion to be drawn. In particular, in June, the uncertainty (±22%) exceeds the relative change (+8%).  Monthly isoprene fluxes were enhanced in June 2012 compared to the previous year ( Figure 5a and Table 1). More favourable (i.e., warm and sunny) environmental conditions at the onset of the severe drought explain this increase, as shown by the model ability to match the large emission enhancement in June (+47% in the observations and 34-60% in the model). In subsequent months, observed isoprene fluxes were about 40% and 50% lower than in the previous year, representing drastic decreases relative to the model expectation ignoring drought stress effects (+16% and −14% in July and August, respectively). The updated MEGAN model based on γ OPT SM reproduces fairly well the seasonal variation (Figure 5a) as well as the inter-annual variability of summer monthly averages (Table 1).
A peak was observed in OMI HCHO in July 2011 at about 22 × 10 15 molecules cm −2 mainly due to high temperatures entailing fast VOC oxidation and strong biogenic emissions. About half this column was observed in August 2012 (12.5 × 10 15 molecules cm −2 ) when severe drought took its toll on biogenic emissions. While the small increase observed in June 2012 with respect to the previous year (+7%) is well simulated in all runs, the substantial drop in July (−19%) is very well represented in MAGRITTE, γ OPT SM (−20%), whereas MAGRITTE, γ SM = 1 and MAGRITTE, γ 0.04 SM barely account for half of that reduction (Table 1). Yet in August, the relative impact of drought in MAGRITTE, γ OPT SM is slightly overestimated (−28% vs. −23% in OMI). While those comparisons point to a better performance of MAGRITTE, γ OPT SM at reproducing the interannual variability of monthly averaged HCHO columns, uncertainties in OMI relative changes (given between brackets in Table 1) do not allow for a robust conclusion to be drawn. In particular, in June, the uncertainty (±22%) exceeds the relative change (+8%).
Further, we assess the absolute difference in HCHO columns of 2012 and 2011 observed by OMI and modelled by MAGRITTE with γ SM = 1 and γ OPT SM ( Figure 6). Two regions depicted in Figure 6, where MAGRITTE predicts significant impacts of the drought stress effect, namely Texas (hereafter TX) and southeast U.S. (hereafter SEUS), are examined. Grid cells that are marked by a stippling are those for which the 2011 and 2012 columns are not statistically different (at the 68% confidence level), based on the random error estimates provided with the HCHO product.
A strong enhancement in 2012 columns over TX is observed by OMI, particularly in June but also in August. Since the area is primarily covered by prairie grasslands, HCHO columns are lower over central and west TX compared to the SEUS. The severe drought that hit Texas during the summer of 2011 (https://www.drought.gov/states/texas, accessed on 1 September 2021) led to enhanced temperatures and PAR ( Figure S3) and to lower humidity levels in comparison with 2012. The strong regional enhancement over TX observed in OMI, +21% of relative change (RC) on average (Table 2), is largely underestimated by the run ignoring drought stress (MAGRITTE, γ SM = 1) (+3%) because the faster VOC oxidation rate in 2012 (due to higher OH concentrations) and higher HCHO photolysis rate in 2011 (due to high radiation levels) were compensated by lower isoprene emissions in 2012 (due to lower temperatures and PAR). The model agreement with OMI greatly improves MAGRITTE, γ OPT SM run 3, when accounting for γ OPT SM . The activity factor decreased biogenic emissions over east Texas in June 2012, leading to enhanced differences (14%) in HCHO columns over the TX region. In subsequent months, the activity factor also improves the prediction of the column differences between 2011 and 2012, from −10% in MAGRITTE, γ SM = 1 to small positive changes. While the average change has been improved due to the inclusion of stress, its spatial distribution has not. For instance, MAGRITTE, γ SM = 1 run 1 is in better agreement with OMI observations than MAGRITTE, γ OPT SM in East Texas (Piney Woods), pointing to an overestimated drought stress effect in this region. This region exhibits relatively high HCHO columns due to the important source of isoprene oxidation stemming from the presence of a mixed evergreen-deciduous forestland and the 2012-2011 absolute difference in temperatures ( Figure S3) plays a dominant role in this interannual variation of isoprene emissions.
umns. Both MAGRITTE-calculated column differences lie within the error bar of the observations. Note that the filtering of OMI data based on AOD and GFED emissions (aimed to minimize the potential impact of biomass burning and aerosols) does not change the conclusions (Section S6 of the Supplement). The filtering leads to a decrease in the number of observations, most pronounced in the SEUS (Figure S5), which increases the noise level and leads to enhanced column differences ( Figure S6). On average over the TX and SEUS regions, however, the relative differences (Table S2) are very similar to those of Table 2.
The southeastern forests, a dense temperate trees coverage, extend from east Texas (Piney Woods) to North Carolina. HCHO columns in SEUS are primarily attributed to the oxidation of non-methane VOC due to the high biogenic isoprene and monoterpene emissions [41,42,62]. The region experienced lower temperatures and PAR in June and August 2012 ( Figure S3) which explains to a great extent the negative 2012-to-2011 HCHO differences observed by OMI and captured by the model, whereas opposite patterns are found in July. OMI observations in the SEUS show slight decreases in 2012 averaged columns in June (−3%) and July (−2%) followed by a larger drop of −14% in August ( Table 2). In June and August, the relative changes simulated by MAGRITTE, γ OPT SM in SEUS amount to −3% and −20%, respectively, i.e., in overall better agreement than MAGRITTE, γ 0.04 SM using the default parameter. The impact of the drought stress factor ( Figure S4) is mainly localized in Louisiana-Arkansas-Missouri and Georgia. This corresponds well with the 2012 patterns reported by the U.S. Drought Monitor (https://www.ncdc.noaa.gov/ sotc/drought/201208, accessed on 1 September 2021).
While average columns and relative changes point to an overall better agreement when using the optimized drought factor (with MAGRITTE, γ OPT SM , the uncertainties in OMI HCHO retrievals undermine the confidence in this result. Based on the random error estimates provided with the QA4ECV HCHO product, only about 50% of the differences shown on Figure 6 are statistically different from zero at 1σ (20% at 2σ i.e., 95% confidence level). Moreover, Table 2 shows that the uncertainty in the OMI relative changes based on the estimated random uncertainties of the OMI retrieval, is of the order of 20%, which is higher than (or of the same order as) the relative changes in OMI HCHO averaged columns. Both MAGRITTE-calculated column differences lie within the error bar of the observations. Note that the filtering of OMI data based on AOD and GFED emissions (aimed to minimize the potential impact of biomass burning and aerosols) does not change the conclusions (Section S6 of the Supplement). The filtering leads to a decrease in the number of observations, most pronounced in the SEUS ( Figure S5), which increases the noise level and leads to enhanced column differences ( Figure S6). On average over the TX and SEUS regions, however, the relative differences (Table S2) are very similar to those of Table 2.

Global Simulations of HCHO Columns
The interannual variability in seasonal averages from global simulations modelled by IMAGES over the course of 2005-2016 is assessed against OMI HCHO columns ( Figure S7). Figure 7 displays the difference in correlation coefficients between two runs, IMAGES, γ 0.04 SM and IMAGES, γ OPT SM . The seasonal averages are calculated over the months of June-August in the Northern Hemisphere (>0 • N), February-April in the 0 • -30 • S latitudinal band, and December-February below the 30th parallel, corresponding to months with low biomass burning activity.
Overall, the correlation coefficient in the interannual variability in HCHO columns is deteriorated in IMAGES, γ OPT SM compared to IMAGES, γ 0.04 SM , particularly, in southeastern U.S., the Gran Chaco (vicinity of Paraguay), eastern and southeastern Europe, and northern Australia. Those regions exhibit relatively high isoprene emissions (e.g., [6]) that lead to high, statistically significant correlations in HCHO in IMAGES, γ 0.04 SM ( Figure S7a). But stronger reductions in isoprene emissions in these areas due to the adjusted parameterization are not supported by OMI observations. Correlations are slightly improved in arid or semi-arid regions such as the Great Plains (central U.S.), shrubland of Western Australia, and the Argentinian Monte desert. Due to the sparse vegetation found in these regions, emissions of isoprene are of a lesser magnitude and the poor agreement between IMAGES, γ 0.04 SM and OMI observations is found in these vicinities ( Figure S7a) is likely due to the relatively low HCHO signal. As a result, improved correlations shown in Figure 7 remain statistically insignificant ( Figure S7b). The filtering of OMI data for high AOD and pyrogenic emissions have only a minimal impact on the conclusions of the global model analysis ( Figure S8). The filtering slightly decreases, as expected, the temporal correlation between the monthly HCHO columns from the model and the observation, due to degraded statistics, but the filtered dataset leads to the same patterns as seen in Figure 7.  Figure S7), respectively.

Discussion and Conclusions
Previous studies conducted on the MOFLUX site noted that the low value of wilting point used in the default parameterization of MEGANv2.1 at this site (0.084 m 3 m −3 , [103]) prevent the stress factor from capturing the effect of drought on isoprene emissions [30,32,52,68]. In this study, we focused on optimizing the empirical parameter Δ of the MEGANv2.1 soil moisture stress algorithm. The parameter Δ was originally set to 0.06 [48] and was later re-evaluated to a lower, more conservative, value of 0.04 m 3 m −3 [5] so that the emissions are shut off only in extreme drought conditions. Here we re-evaluated it (Δ = 0.12) based on isoprene flux measurements at the MOFLUX site conducted in mild and severe drought conditions, and GLEAM soil moisture. The adjusted parameterization leads to an improved prediction of daily isoprene emissions rates during the severe drought at the MOFLUX site (r = 0.90, RMSE = 2.02 mg m −2 h −1 ) compared to the nostress simulation (r = 0.76, RMSE = 4.98 mg m −2 h −1 ) and to the simulation adopting the MEGANv3 parameterization (r = 0.79, RMSE = 3.13 mg m −2 h −1 ). Analysis of OMI and modelled tropospheric HCHO column over the MOFLUX site supports this conclusion. At a regional scale, the relative changes in simulated formaldehyde between the two years show some improvement when including the adapted parameterization over certain regions in the U.S., but clearly not over all areas. Nonetheless, global long-term HCHO simulations using lead to heterogeneous results, with extended areas where the correlation with OMI observations worsens compared to the use of the default parameterization.
The study is subject to limitations and uncertainties. Foremost, the adjustment of the parameter Δ is influenced by many factors including the dataset, the soil depth of reference, and the MEGAN parametrized response functions to environmental conditions other than soil moisture. For instance, using in situ measurements of the at 100 cm and the local wilting point (0.23 m 3 m −3 ), the optimized Δ would be about twice higher (0.25 Figure 7. Difference in correlation coefficients of simulated HCHO columns with OMI HCHO columns at 2 • × 2.5 • spatial resolution: R 2 − R 1 , where R 1 and R 2 are the correlation coefficients from IMAGES, γ 0.04 SM and IMAGES, γ OPT SM simulations (displayed in Figure S7), respectively.

Discussion and Conclusions
Previous studies conducted on the MOFLUX site noted that the low value of wilting point used in the default parameterization of MEGANv2.1 at this site (0.084 m 3 m −3 , [103]) prevent the stress factor from capturing the effect of drought on isoprene emissions [30,32,52,68]. In this study, we focused on optimizing the empirical parameter ∆θ of the MEGANv2.1 soil moisture stress algorithm. The parameter ∆θ was originally set to 0.06 [48] and was later re-evaluated to a lower, more conservative, value of 0.04 m 3 m −3 [5] so that the emissions are shut off only in extreme drought conditions. Here we re-evaluated it (∆θ = 0.12) based on isoprene flux measurements at the MOFLUX site conducted in mild and severe drought conditions, and GLEAM soil moisture. The adjusted parameterization leads to an improved prediction of daily isoprene emissions rates during the severe drought at the MOFLUX site (r = 0.90, RMSE = 2.02 mg m −2 h −1 ) compared to the no-stress simulation (r = 0.76, RMSE = 4.98 mg m −2 h −1 ) and to the simulation adopting the MEGANv3 parameterization (r = 0.79, RMSE = 3.13 mg m −2 h −1 ). Analysis of OMI and modelled tropospheric HCHO column over the MOFLUX site supports this conclusion. At a regional scale, the relative changes in simulated formaldehyde between the two years show some improvement when including the adapted parameterization over certain regions in the U.S., but clearly not over all areas. Nonetheless, global long-term HCHO simulations using γ OPT SM lead to heterogeneous results, with extended areas where the correlation with OMI observations worsens compared to the use of the default parameterization.
The study is subject to limitations and uncertainties. Foremost, the adjustment of the parameter ∆θ is influenced by many factors including the θ dataset, the soil depth of reference, and the MEGAN parametrized response functions to environmental conditions other than soil moisture. For instance, using in situ measurements of the θ at 100 cm and the local wilting point (0.23 m 3 m −3 ), the optimized ∆θ would be about twice higher (0.25 m 3 m −3 ) than the value derived using GLEAM θ (0.12 m 3 m −3 ). In addition, the high sensitivity of isoprene emissions to soil properties using the MEGANv2.1 parameterization has been demonstrated in numerous studies and motivated the adoption of higher values of ∆θ despite the earlier rationale to only apply the isoprene reduction for extreme droughts, that led to its downward revision. Also, soil moisture should be seen as relative, i.e., framed by the soil properties the model uses. In particular, the high heterogeneity of porosity renders soil properties inaccurate since they are a function of porosity, such that the comparison of soil moisture in absolute terms is of limited relevance. Consistency in the use of soil characteristics dataset is also crucial as already pointed out in [49]. For year 2012, Seco et al. (2015) [32] used the MEGANv2.1 drought activity factor algorithm with the measured θ data at 10 cm jointly with the local wilting point (0.23 m 3 m −3 ) and found a substantial improvement in the overall model results against the MOFLUX isoprene flux measurements. However, the wilting point of 0.13 m 3 m −3 for the top 30 cm (Section 2.1) would be more appropriate when using the 10 cm θ data. Therefore, since the θ measured at 10 cm always exceeds the threshold of 0.13 + 0.04 = 0.17 m 3 m −3 , no drought stress is predicted by MEGANv2.1 based on 10 cm θ data. The level at which the θ is defined is important due to the general increase of θ with depth and to changes in the temporal evolution due to water percolation.
Regarding the assessment of drought impacts using formaldehyde columns, the following shortcomings and limitations are raised. OMI HCHO observations lack the precision required to conduct an evaluation on isoprene emission estimates at short temporal scales (daily and weekly), which would be most relevant for drought impact studies. With its much higher resolution and increased sensitivity to day-to-day variability [102], the TROPOMI sensor on Sentinel-5p offers promising perspectives for future research. As we have seen, the exclusion of OMI data characterized by high AOD or high pyrogenic emissions degrades the statistics while not fundamentally affecting the conclusions of the study. Modelled tropospheric HCHO columns over the MOFLUX site pointed to an overall better agreement with the adjusted parameterization and the radius of 20 km was chosen as the best for representing the site given that satellite overpasses within a smaller radius were too infrequent (<17/month for 10 km). At greater coverage (30 km), however, the MAGRITTE runs with the least soil moisture stress (as in MAGRITTE, γ SM = 1 and MAGRITTE, γ 0.04 SM ) are in best overall agreement with OMI observations (Table S1). At the global scale, high values of ∆θ are not supported by observations. Indeed, multi-year global simulations (2005-2016) with the CTM IMAGES indicated that the use of ∆θ = 0.12 m 3 m −3 does not improve significantly the comparison with OMI HCHO over dry regions, and it worsens the comparison over most forested areas. The best agreement with OMI in terms of interannual variability is found when the drought activity factor is not accounted for (not shown here). The reason for such disparity between the MOFLUX study and the large-scale evaluation using HCHO could lie in the spatial resolution. Disparities in soil water content between local and grid cell scales can be huge given the great heterogeneity in soil properties like porosity and subsequently for the 'absolute' values of soil moisture and the empirical parameter. Another issue with the current parameterization is that it is independent of various environmental factors other than the soil moisture characteristics. Furthermore, the response of isoprene emissions depends on the drought characteristics and plant species. The parametrization of drought stress should account for different effects across biomes and species [104][105][106]. For instance, Genard-Zielinski et al. (2018) [39] suggested an exponential reduction function of the SWC for drought-adapted species. In addition, longer and recurrent drought (e.g., 10 years) could have progressive and cumulative effects over time in forested ecosystems dominated by long-lived species [18]. Besides the complexity of the drought response, the lack of ecosystem-scale, long-term monitoring of soil moisture and isoprene emissions is another important issue. Direct comparison between MEGAN and flux measurements in various drought conditions has been only possible for the MOFLUX site. Other campaigns including isoprene concentration measurements in drought-impacted ecosystems ( [39,40]) provide qualitative insights but no direct constraint on the fluxes. A possible avenue for drought stress parameterizations could be the development of empirical relationships making use of remotely sensed photosynthetic parameters or vegetation indices such as SIF [47], Emissivity Difference Vegetation Index (EDVI; [107]) or Photochemical Reflectance Index (PRI; [108]). So far, however, progress has been limited in that direction.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/rs14092021/s1, Figure S1: Leaf area index at the MOFLUX site; Figure S2: Temporal frequency for the analysis of isoprene and HCHO columns at MOFLUX; Figure S3: Difference in maximal temperature and PAR over 2011 and 2012; Figure S4: Distribution of the MEGANv2.1. drought activity factor over eastern U.S.; Table S1: Sensitivity of OMI HCHO columns to the radius around the MOFLUX site; Figure S5, Figure S6 and Table S2: Screening of OMI columns for effects of fires and aerosols; Figure S7: Global correlation coefficients with OMI HCHO columns. Reference [109] is cited in the Supplementary Materials..