Articles | Volume 14, issue 8
Research article
13 Aug 2021
Research article |  | 13 Aug 2021

Boundary layer water vapour statistics from high-spatial-resolution spaceborne imaging spectroscopy

Mark T. Richardson, David R. Thompson, Marcin J. Kurowski, and Matthew D. Lebsock

Daytime clear-sky total column water vapour (TCWV) is commonly retrieved from visible and shortwave infrared reflectance (VSWIR) measurements, and modern missions such as the upcoming Earth Surface Mineral Dust Source Investigation (EMIT) offer unprecedented horizontal resolution of order 30–80 m. We provide evidence that for convective planetary boundary layers (PBLs), spatial variability in TCWV corresponds to variability in PBL water vapour. Using an observing system simulation experiment (OSSE) applied to large eddy simulation (LES) output, we show that EMIT can retrieve horizontal variability in PBL water vapour, provided that the domain surface is uniformly composed of either vegetated surfaces or mineral surfaces. Random retrieval errors are easily quantified and removed, but biases from −7 % to +34 % remain in retrieved spatial standard deviation and are primarily related to the retrieval's assumed atmospheric profiles. Future retrieval development could greatly mitigate these errors. Finally, we account for changing solar zenith angle (SZA) from 15 to 60 and show that the non-vertical solar path destroys the correspondence between footprint-retrieved TCWV and the true TCWV directly above that footprint. Even at the 250 m horizontal resolution regularly obtained by current sensors, the derived maps correspond poorly to true TCWV at the pixel scale, with r2<0.6 at SZA=30. However, the derived histograms of TCWV in an area are closely related to the true histograms of TCWV at the nominal footprint resolution. Upcoming VSWIR instruments, primarily targeting surface properties, can therefore offer new information on PBL water vapour spatial statistics to the atmospheric community.

1 Introduction

Thermodynamic information about the planetary boundary layer (PBL), including information about water vapour (qv), is a targeted observable recommended by NASA's Decadal Survey (National Academies of Science, Engineering, and Medicine, 2018). PBL qv estimates would go beyond the current total column water vapour (TCWV) and free-tropopause products to provide new information about the vertical moisture structure for weather and climate applications. The Decadal Survey explicitly recognised the PBL's importance since it “literally couples the surface of the Earth to the atmosphere above”, and among other important factors, gradients of moisture between the surface and PBL and between the PBL and free troposphere are strong controls on vertical atmospheric heat and moisture transport. The formation of boundary layer clouds was also highlighted due to their importance for Earth's energy balance. A critical measurement gap in the current observations of PBL thermodynamics is the inability to quantify mesoscale variations in PBL qv. Mesoscale aggregation in PBL water vapour appears to play an important role in determining the timing of deep convective events (Stirling and Petch, 2004; Wulfmeyer et al., 2006). Furthermore, in situ observations suggest that the majority of the variation in the TCWV prior to convective initiation can be explained by variability within the PBL (Couvreux et al., 2009). The mesoscale spatial variability of qv is not resolved by current global weather or climate models, but instead it must be parameterised. Modern approaches to parameterise PBL variability include eddy-diffusivity/mass-flux approaches (Suselj et al., 2019) and higher-order closure approaches that include prognostic equations for higher-order moments such as the variance (Golaz et al., 2002; Larson et al., 2002). However, we lack observations at a global scale to evaluate the small-scale variability produced by these models. This paper will address the feasibility of addressing this measurement gap using upcoming observations from very high-spatial-resolution visible and shortwave infrared reflectance (VSWIR) observations from space.

This study is primarily motivated by the ongoing development of spaceborne hyperspectral VSWIR measurement capacity at fine horizontal resolution. We focus on the EMIT mission, planned to launch to the International Space Station (ISS) in 2022 with an average footprint size (Δx) of 60 m (Green and Thompson, 2020). However, similar or improved capacity is anticipated in response to NASA's Surface Biology and Geology (SBG) designated observable, with the Hyperspectral Infrared Imager (HysPIRI, Lee et al., 2015) concept offering Δx of 30–60 m, and for ESA's Copernicus Hyperspectral Imaging Mission for the Environment (CHIME), also known as Sentinel 10, for which the prime contractor was selected in July 2020 and whose Mission Requirements Document refers repeatedly to Δx<30 m (Rast et al., 2019).

Of present missions, this analysis may be applicable to the Italian PRecursore IperSpettrale della Missione Applicativa (PRISMA, Candela et al., 2016), which provides similar spectral range and sampling to EMIT at Δx=30 m. Some of the conclusions will also apply to other recent instruments, such as Sentinel-2's Multi-Spectral Imager (MSI, Drusch et al., 2012), which offers Δx=20 m, albeit with far fewer channels, or the DLR Earth Sensing Imaging Spectrometer (DESIS, Krutz et al., 2019), which provides hyperspectral measurements over a smaller wavelength range. These modern and upcoming instruments offer Δx that are substantially smaller than past VSWIR instruments that retrieve TCWV, such as ESA's Medium Resolution Imaging Spectrometer (MERIS) on Envisat, whose smallest provided Δx is approximately 0.25 km×0.30km, which allowed the identification of horizontal convective rolls during a high-pressure event over Germany (Carbajal Henken et al., 2015) but cannot resolve the smaller scales of variability. Recently, Thompson et al. (2021) used VSWIR measurements from the Airborne Visible Infrared Imaging Spectrometer-Next Generation (AVIRIS-NG) to capture information about PBL qv variability at spatial scales <1 km, which cannot be determined with footprint sizes similar to MERIS.

EMIT-like instruments could allow retrieval of bulk PBL qv, which we henceforth refer to as the partial column water vapour in the PBL (PCWVPBL) via two demonstrated approaches. The first approach uses VSWIR measurements alone, and the second combines separate above-PBL water vapour (PCWVupper) and TCWV to obtain PCWVPBL=TCWV-PCWVupper. A third approach, which has not been demonstrated operationally to our knowledge, is to perform joint retrievals using both VSWIR and vertically resolved sounding measurements.

The direct VSWIR-only method can be seen in Trent et al. (2018), who estimated PCWVPBL from the Greenhouse Gases Observing Satellite (GOSAT, Kuze et al., 2009), while the second is explored in Millán et al. (2016), who paired TCWV from passive microwave measurements with PCWVupper above horizontally uniform clouds from Moderate Resolution Imaging Spectroradiometer (MODIS) near-infrared retrievals. The resultant PCWVPBL values showed good agreement with radiosondes and ERA-Interim reanalysis, and a promising candidate approach is to use VSWIR TCWV in place of the microwave measurements.

The physical principle of VSWIR TCWV retrievals is differential optical absorption spectroscopy (DOAS). More TCWV leads to increasing depth of H2O absorption features relative to other wavelengths. This applies to TCWVVSWIR from missions including MERIS (Bennartz and Fischer, 2001; Guanter et al., 2008), MODIS (Diedrich et al., 2015; Gao and Kaufman, 2003), TROPOMI (Borger et al., 2020; Schneider et al., 2020), SCIAMACHY (Noël et al., 2004), GOME (Noël et al., 1999), GOME-2 (Grossi et al., 2015) and OCO-2 (Nelson et al., 2016).

These instruments vary in spectral range and sampling, but all must contend with the measured spectra responding to properties other than TCWV. The retrievals only operate for daytime cloud-free scenes and commonly only over land, since water surfaces are dark such that insufficient light reaches the sensor to allow for a TCWV retrieval, with exceptions for sun glint as exploited in the aforementioned AVIRIS-NG study (Thompson et al., 2021). Thompson et al. selected these AVIRIS-NG flights because DOAS techniques respond to the total light path absorption including the slanted sunlight path from the top of atmosphere (TOA) to the surface. This horizontally smears the effective footprint size, with larger smearing for larger solar zenith angle (SZA). As footprints become smaller, the proportional effect of this smearing may become more important, and so here we apply solar ray tracing to determine whether observations with a nominal Δx of 20–50 m obtain useful information about the spatial statistics of PCWVPBL at that spatial resolution. We use two performance metrics: (i) the correlation between retrieved TCWV and true TCWV, which was used as input for our forward simulations, and (ii) the spatial standard deviation σx of retrieved TCWV within a snapshot relative to the large eddy simulation (LES) output PCWVPBLσx, which we refer to as the true σx.

We employ a new type of Observing System Simulation Experiment framework and perform simulated VSWIR retrievals of TCWV from high-spatial-resolution LES output to determine whether horizontal spatial variability in PBL qv can be obtained from retrieved TCWV, and conclusions are limited to daytime non-cloudy conditions. The purpose of this is a detailed sensitivity study using retrieval code and tools already developed for EMIT. We consider Δx≥40 m since this is appropriate for EMIT and several LES cases in our archive were run at that resolution.

Here we test the use of the iterative optimal estimation code Imaging Spectrometer Optical Fitting (Isofit) for a spaceborne application, specifically target TCWV and address the following questions.

  1. In LES, how does horizontal variability in TCWV relate to PCWVPBL?

  2. What uncertainties are introduced into the retrieval by EMIT instrumental error, non-uniform AOD and different surface types, and can these errors be anticipated and quantified from observations alone?

  3. What is the correlation coefficient between retrieved and true TCWV, and can the spatial standard deviation be estimated? How does this depend on LESs of different convective PBL types?

  4. How does the solar path across different SZAs affect these conclusions?

This scope excludes important factors such as topography, inter-channel correlated errors in the instrument, imperfect cloud masking and cloud 3D radiative effects, and our paper is structured to address these questions in turn, with each analysis section containing its own methodology, results and discussion. Section 2 explores the raw LES output to address question 1, Sect. 2 describes the synthetic retrievals and analysis methodology to address questions 2–3, Sect. 4 adds solar path analysis to address question 4, and Sect. 5 discusses and concludes.

2 Large eddy simulations

2.1 Model setup, scenarios and snapshot selection

We use output from five LES simulations named RICO, ARM, ARM_lsconv, BOMEX and DRY, which are summarised with references in Table 1. They all represent convective boundary layers characterised by either low-altitude or no cloud cover. The 23 separate snapshots are identified by timestamp; e.g. ARM_18000s is 5 h into the ARM simulation. Simulation Δx sets the implied measurement horizontal resolution and varies from 20 to 50 m.

The simulations are performed with two different models: EULAG (Prusa et al., 2008; ARM and ARM_lsconv) and JPL-UCONN LES (Matheou and Chung, 2014; RICO, BOMEX, DRY). Each simulation applies periodic lateral boundary conditions and a horizontally homogeneous initial state. For the RICO case, interactive sensible and latent heat surface fluxes over constant-temperature ocean are used, while the other cases are driven by prescribed (either constant for DRY and BOMEX or time-dependent for ARM and ARM_lsconv) surface fluxes. All other setup details are explained in the Table 1 references: these references also show how the ARM, BOMEX and RICO LES simulations, which were based on detailed field campaigns, accurately reproduce the main features observed during those campaigns. Each 3D LES snapshot is merged with 1D MERRA-2 reanalysis profiles aloft to produce a full-depth atmospheric column. Reanalysis data are chosen for the dates and locations of the field campaigns the LESs refer to. The DRY and ARM_lsconv cases share the same upper-atmospheric profiles as ARM. In all cases except for DRY, Table 1 rows (vii)–(ix) show that the LESs capture >85 % of total TCWV. For retrieval purposes we ignore the LES surface type and apply an assumed surface reflectance spectrum below the LES profiles.

Table 1Summary of LES properties. Where ranges are provided, these are the full range of clear-sky mean values from the snapshots used for each case. Row (vi) is the fraction of columns whose integrated liquid water path >1.3×10-3 mm and differs from mean cloud fraction in Fig. 1 due to overlap. The TCWV in row (vi) is derived from the combined LES and reanalysis profile and separated into the LES and reanalysis partial column water vapour amounts in rows (vii) and (viii).

Download Print Version | Download XLSX

2.2 Profiles and PBL height

Definitions of PBL height zPBL vary widely. We found similar results from four standard thermodynamic calculations (von Engeln and Teixeira, 2013), so henceforth we define zPBL as the altitude of max(dθ/dz), where θ is the all-sky mean potential temperature. Mean all-sky profiles of T and q, horizontal standard deviation in q(σq), and cloud fraction are shown in Fig. 1. Changes in σq are the largest differences between time steps but are small (<10 %) relative to the mean, so measuring this variability will require precise observations. Also, σq is negligible in the layers in the free troposphere that lie above the PBL but are resolved by the LES, implying that the LES domains capture qv variability. We later support this claim using real-world airborne lidar retrievals.

Figure 1Output all-sky profiles for the (a–d) ARM, (e–h) ARM_lsconv, (i–l) RICO, (m–p) BOMEX and (q–t) DRY LES. In each panel the separate coloured lines represent different timesteps, the black horizontal line is the top of the LES and the dashed blue horizontal line is the PBL height calculated from the first shown timestep, whose lines are in the same blue. The column beginning with (a) is the mean T profile, that with (b) the mean q profile, that with (c) the profile of the spatial standard deviation in q and that with (d) the cloud fraction. Note that due to overlap, the fraction of cloudy columns listed in Table 1 is higher than the peak mean profile cloud fraction.


2.3 Water vapour spatial variability statistics and the relationship between TCWV and PCWVPBL

Figure 1 displays all-sky conditions, but our retrievals only target clear sky, thereby missing a moister tail to the distribution (Supplement Fig. 1). Within-cloud retrievals would require alternative measurement approaches, such as differential absorption radar (Roy et al., 2018, 2020), and the restriction to clear-sky scenes is a limitation that also applies to current thermal infrared and lidar retrievals.

We assess TCWV-PCWVPBL spatial variability by calculating clear-sky PCWV up to capping altitudes from 0.5 to 5 km and then correlating these with TCWV. Figure 2 confirms that >90 % of horizontal variance in LES TCWV at these scales is explained by PCWVPBL. It is reasonable to ask whether this finding that the PBL variance dominates the TCWV variance is representative of the real atmosphere. Indeed, the LES results are supported by the same statistics calculated from High Altitude Lidar Observatory (HALO) flights over the Pacific Ocean in April 2019 (Bedka et al., 2021), as presented in Thompson et al. (2021) and shown in Fig. 2f. In these calculations TCWV is only calculated up to 8 km due to flight altitude, but these real-world data include free-tropospheric moisture variability and furthermore will have lower r values due to the presence of random retrieval error. The horizontal resolution is ∼3 km versus the 20–50 m of LES, and the HALO sampling is sparse and often separated by hundreds of kilometres due to clouds. Nevertheless, the HALO flights show that horizontal TCWV variability can be well captured within 3 km altitude in real scenes and provide evidence that the LES domains capture horizontal variability in qv.

Figure 2Correlation coefficient between clear-sky partial column water vapour (PCWV) integrated up to given capping altitudes, and the TCWV. (a)(e) contain the snapshots of each individual LES run and (f) reproduces the values calculated from High Altitude Lidar Observatory flights over the Pacific (Bedka et al., 2021) as presented in Thompson et al. (2021). The LES profiles also have a horizontal bar appended at the derived PBL top height. The flight data differ from the LES outputs in that horizontal resolution is approximately 3 km along-track, they are dispersed over thousands of kilometres, and the TCWV is only up to 8 km due to the flight altitudes.


The TCWV-PCWVPBL fit coefficients for ARM, ARM_lsconv, BOMEX and RICO range from 0.99 to 1.04 mm mm−1; i.e. a 1 mm change in PCWVPBL means a 0.99–1.04 mm change in TCWV. This confirms that almost all horizontal qv variability occurs within the mean PBL height. For the DRY case, coefficients range from 1.06 to 1.12 mm mm−1. These coefficients mean that PCWVupper spatially correlates with PCWVPBL, which could be explained by moister plumes rising and having higher local zPBL than the domain-mean value used in the calculation. In summary, we have answered question 1 from Sect. 1 and can expect spatial variability in retrieved TCWV for these cases to represent real variability in PCWVPBL, and so we use TCWV and PCWVPBL interchangeably from now on.

3 Simulated EMIT retrievals of TCWV in LES

This experiment requires a large number of inversions over a wide spatial field. Simulating synthetic spectra and performing a retrieval for every grid point proved to be prohibitively computationally expensive. Consequently, we develop an emulator to statistically reproduce the result of the full inversion but with dramatically better efficiency. Retrievals will include a range of surfaces in a subset of the snapshots (to identify sensitivity to surface type) and then a fixed surface type across all snapshots (to identify sensitivity to atmospheric conditions). Sensitivity tests will be performed on individual subsets of snapshots as required, and a correction method for identifying the random component of retrieval error will be introduced. Section 3.1 describes the relevant methods, Sect. 3.2 gives the results and Sect. 3.3 discusses limitations.

3.1 Retrieval methodology

3.1.1 MODTRAN6.0 forward model, EMIT instrument characteristics and Isofit retrievals

We use the same retrieval code as in Thompson et al. (2021), Imaging Spectrometer Optimal Fitting (Isofit), for our synthetic retrievals (, last access: 2 June 2021). This iterative optimal estimation code simultaneously retrieves surface reflectance, aerosol optical depth (AOD) and TCWV, differing from older techniques that retrieve properties sequentially (e.g. Guanter et al., 2008 for MERIS). Isofit is described and shown to have a closed error budget in Thompson et al. (2018) and has been applied to observations from several airborne campaigns (Thompson et al., 2019, 2020, 2021). Conceptually it targets surface reflectance ρs, and the estimation of TCWV is seen as part of an atmospheric correction.

Forward simulations use MODTRAN6.0 (Berk et al., 2014, 2015), which provides a plane-parallel solution to the radiative transfer equation. Atmospheric reflectance and transmittance vectors ρa, t and spherical sky albedo s are calculated at wavenumber separation Δk=0.1cm−1 (Δλ≈0.002 nm) before being convolved with the EMIT spectral response function, and the instrument is assumed to be nadir-viewing from 100 km altitude. With no substantial atmosphere above 100 km, this gives the same results as the ISS altitude near 400 km, where EMIT will be hosted. A correlated-k method and the HITRAN database (Rothman et al., 2009) are used for gaseous absorption, while scattering is handled by DISORT (Laszlo et al., 2016; Stamnes et al., 1988). The EMIT instrument properties are derived from the current mission instrument model, which accounts for all signal-independent noise terms like electronic noise, and photon shot noise calculated using predicted efficiencies of the instrument mirrors, lens, grating, and focal plane array. The spectral range is 380–2500 nm, with Δλchannel=10 nm and full width at half maximum averaging ΔλFWHM≈11 nm.

For forward simulations, merged LES-reanalysis T and q profiles are interpolated onto a profile with eight points from 0 to 6 km, and then vertical resolution slowly degrades over 6–100 km. Interpolated TCWV differs from the LES reanalysis, but we assume that conclusions regarding derived sensitivities and errors will not be strongly affected.

The forward radiance vector I is calculated using a standard Lambertian approximation (e.g. as in Vermote et al., 1997):

(1) I = I 0 μ 0 π ρ a + t ρ s 1 - s ρ s ,

where I0 is the downward TOA solar radiance, μ0 the cosine of the solar zenith angle, ρs is the surface reflectance and represents channel-by-channel multiplication. The ρs elements represent the hemispheric-directional distribution function (Schaepman-Strub et al., 2006). The atmospheric coefficient vectors t, ρa and s represent the transmittance of the solar-reflected optical path, the path reflectance, and the spherical sky albedo, respectively. These coefficients are obtained from simulations over a black surface. Using Eq. (1) in forward simulations results in negligible differences to retrieved TCWV compared with inserting the surface directly into MODTRAN forward simulations (Supplement Fig. 2). Use of Eq. (1) means that just one MODTRAN simulation is needed per column rather than one for each combination of column and surface type. The pseudo-observation, Iobs, is I with random uncorrelated noise added, generated using the EMIT noise model.

The Iobs are input as observations to Isofit, while its state vector x contains surface reflectance in each channel, TCWV and aerosol optical depth at λ=550 nm (AOD), i.e. x=[ρs  AOD  TCWV]. We mask the most strongly absorbing channels due to lack of any surface information, so the retrieval uses 176 EMIT channels, and therefore x has 176+2=178 elements.

The ρs elements are constrained via a covariance matrix whose mean is derived from a library of real surfaces, thereby capturing realistic spectral shapes. We retrieve absolute ρs, rather than the normalised value discussed in Thompson et al. (2018), and the prior is loosely constrained, however, ensuring that most information comes from the measurements.

Isofit uses Eq. (1) with a lookup table (LUT) for its forward model, populating ρa, t, and s for selected AOD and TCWV and then linearly interpolating in TCWV, AOD space to estimate Iobs given x. The LUT uses the default midlatitude summer profile and scales its q(z) and aerosol extinction profile to match desired AOD (from 0.05 to 0.30) and TCWV (from 5 to 53.5 mm). The Isofit default configuration uses the U.S. Standard Atmosphere 1976 (Sissenwine et al., 1976), but MODTRAN applies a relative humidity limit, and the U.S. Standard Atmosphere 1976 is cool enough that MODTRAN automatically restricts its moisture content, such that the TCWV cannot reach the values seen in any LES case except for DRY. The midlatitude summer TCWV limit is just over 53.5 mm, so that defines our LUT maximum.

Our prior and first guess TCWV is 40 mm with a standard deviation of 7.5 mm, although observationally a heuristic band ratio is commonly used to provide a first guess and a locally appropriate prior would be selected. This choice of prior does not change our derived spatial statistics (Supplement Fig. 3), although it results in a small shift of mean retrieved TCWV and reflectance (e.g. posterior TCWV shifts by 0.15 mm when the prior is shifted by 32.5 mm).

3.1.2 Profile subsets, emulator development and sensitivity tests

All retrievals use radiances simulated at SZA=45, using the profiles associated with an individual footprint and assuming a plane-parallel atmosphere. We define “clear sky” as where cloud water path <1×10-3 mm, approximately τ<0.3 in a typical sub-adiabatic cloud (e.g. Szczodrak et al., 2001). Clear-sky columns are ranked by TCWV, and 101 columns equally spaced in terms of this ranking are taken (Supplement Fig. 4 justifies N=101).

All snapshots in a given LES case are combined and Isofit-retrieved TCWVret is used to fit an emulator in combination with the forward-model TCWV via

(2) TCWV ret = a 1 TCWV + a 2 + ϵ ,

where a1 and a2 are the trend and intercept parameters from an optimised-least-squares fit and ϵ is random zero-centred Gaussian noise with standard deviation from the emulator fit residuals. Tests with SZA from 14 to 60 show no significant differences in a1 with SZA, while the standard deviation of ϵ increases by up to 25 % at SZA=60 relative to SZA=45 (Supplement Fig. 5, Supplement Table 1). Section 3.1.4 shows how we are able to identify and remove the effect of ϵ on derived statistics, so given that a1 did not change with SZA in these tests, we anticipate that our conclusions will largely apply to SZA up to and including 60.

Forward-simulation AOD varied from 0.1 to 0.2, and most footprints were assigned AOD=0.2. Supplement Figs. 6–7 show weak sensitivity of retrieved TCWV to AOD. The analysis is separated into two parts: Sect. 3.2.1 shows results for sensitivity of TCWVret to changes in the surface spectrum within selected ARM snapshots and Sect. 3.2.2 shows changes in retrieved TCWV over a single surface type for all snapshots.

3.1.3 Development and fitting of the retrieval emulator

For each emulator we use all snapshots within an LES run to fit Eq. (2) (separate snapshot fits do not affect conclusions, Supplement Fig. 8), and full-snapshot fields of TCWVret are then emulated using Eq. (2) with LES TCWV as input. The surface analysis uses the first three ARM snapshots and seven surface spectra from the Isofit surface model clusters, three of which are typical of vegetation and the others of mineral surfaces. The database used to generate the surface model includes artificial surfaces, which are largely captured by the “mineral” spectra. An additional test was run with ARM_18000s profiles over uniform Lambertian surfaces with ρs=0.1–0.5 in increments of 0.1. The atmospheric analysis uses the MODTRAN cropland and ocean ρs spectra for all 23 snapshots, although poor performance over dark surfaces means that the main emulator results are reported only for the land-surface retrievals.

Figure 3a shows typical spectra simulated over several surfaces: notably, the MODTRAN ρs spectra have sharp changes that are not included in the Isofit surface model and therefore provide a challenging test of the retrieval code's ability to retrieve TCWV outside of the surface conditions for which it was developed.

Figure 3Examples of (a) simulated spectra and (b) used surface reflectances in the forward model (solid lines) and those retrieved by Isofit using EMIT instrument characteristics (dashed lines). Each colour refers to a surface type as listed in the panel (a) legend.


With regards to the emulator parameters, non-unity a1 represents biases in the local retrieval sensitivity dTCWVret/dTCWV. Possible causes will be discussed in Sect. 3.2.3, but this is the main concern for retrieval of local variability statistics because the retrieved standard deviation will be scaled by a1, and this scaling will be undetectable in the absence of independent validation data. Changes in a1 also change the derived spatial r2, since a1>1 increases retrieved σx variance and will increase r2. The parameter a2 is related to a combination of the mean bias and the magnitude of a1 within a snapshot and may depend on factors such as surface type or biases in the LUT-assumed T and q profiles as seen for MERIS retrievals in Lindstrot et al. (2012). For our spatial statistics, a2 has no effect since it is subtracted during calculation. The parameter ε represents non-systematic errors within a scene.

Importantly, σε is not the typical error seen in validation or inter-comparison exercises (Diedrich et al., 2015; Nelson et al., 2016; Pérez-Ramírez et al., 2014), since in these studies the varying biases between products in different conditions will add to the reported errors and make them larger than the σε appropriate for our retrievals.

3.1.4 Estimating random error from retrieved fields

Random retrieval error ε with standard deviation σε adds variance and therefore reduces r2 while adding a high-bias term to estimated σx. Knowing σε would allow removal of its bias contribution to σx, and clearly interpretation of spatial variability at a footprint level requires that σε is small relative to σx. TCWV variability between columns separated by 50 m in the horizontal is far smaller than at larger separations. We will exploit this to estimate the spatially constant σε using an approach based around the second-order structure function S2. Here we describe the recipe and mathematical justification; see Supplement Figs. 9–10 for a step-by-step illustration. For a TCWV field,

(3) S 2 ( Δ r ) = E [ ( TCWV ( x + Δ r ) - TCWV ( x ) ) 2 ] .

This is the variance between pairwise footprints separated by the distance Δr, and retrieved S2 includes contributions from the spatial variance characteristic at that separation, σx2(Δr), and the observational uncertainty σϵ2. The subtraction removes the field mean TCWV, and each of the terms TCWV(x) and TCWV(xr) will contribute σx2(Δr)+σϵ2 to the variance. We treat these as independent, so their variances add to give the retrieved S2,ret:

(4) S 2 , ret ( Δ r ) = 2 σ x 2 ( Δ r ) + 2 σ ϵ 2 .

For ARM_18000s, σx(Δr=50 m) is 0.03 mm, compared with the full-snapshot σx of 0.29 mm. We exploit the smallness of σx at small Δr by smoothing the field in one direction with no overlap between smoothed footprints and then calculating the structure function at Δr=1 footprint (20–50 m, depending on the LES) perpendicular to the smoothing direction. For n-footprint smoothing, the independent component of variance shrinks by 1/n, which we attribute to σϵ2. The steps are the following.

  • (i)

    Select a direction and evaluate S2r) in that direction for Δr=1 footprint separation.

  • (ii)

    Smooth the field in the direction perpendicular to Δr by averaging over nfoot=2 footprints.

  • (iii)

    Recalculate S2(Δr,nfoot=2), treat the calculated value (i) as S2(Δr,nfoot=1), regress S2r,nfoot) against 1/nfoot, and take the best-fit trend to be equivalent to 2σϵ2.

By smoothing in one direction and then calculating orthogonally, the separation distance Δr does not grow with smoothing, and so we maintain the advantages of the small σx(Δr=20–50 m). To estimate TCWV σε with EMIT-like Δx, this method outperforms a standard spatial smoothing filter approach (Supplement Fig. 11).

3.1.5 Calculating spatial statistics and relationship with spatial smoothing

We calculate the spatial standard deviation σx of clear-sky TCWV and TCWVret for each snapshot. The random error σε is then estimated following Sect. 3.2.3 and subtracted in quadrature:

(5) σ x ,ret,corr = σ x ,ret 2 - σ ϵ ,ret 2 ,

where the subscript “ret” means retrieved and “corr” means corrected.

The other target statistic is r2 between TCWV and TCWVret; we calculate this directly and also estimate it via

(6) r est 2 = σ x ,ret 2 - σ ϵ ,ret 2 σ x ,ret 2 .

Where emulator trend a1=1, estimated error from Sect. 3.2.3 is accurate, and there are no spatially varying errors, Eq. (6) should reproduce retrieval r2. However, a1≠1 means each σx,ret2 term will be multiplied by a12, resulting in an erroneous r2 estimate. User requirements for r2 will depend on application; we arbitrarily select r2=0.9 as a target and compare Eq. (6) estimates with the true field values. True r2 is unknowable without perfect knowledge of the TCWV field, but operational estimation using Eq. (6) would allow users to determine whether their requirements are likely to be met.

If r2 is too low for the desired application, then averaging over footprints may address this; although it results in loss of fine spatial information, it may be necessary to suppress errors or may be enforced by effective horizontal smearing where SZA>0.

We show the results of sequentially smoothing the TCWV and TCWVret field on both σx and r2 and smooth in both horizontal directions, for example going from 50 m×50 m to 100 m×100 m. Smoothed footprints do not overlap and so are independent, and the smoothing is done on TCWVret rather than on the radiance field. This avoids the requirement for additional forward-model runs and furthermore allows predictions of how r2 changes with resolution by applying Eq. (6) with a minor modification:

(7) r est 2 = σ x ,ret 2 - σ ϵ ,ret 2 n σ x ,ret 2 ,

where n is the number of footprints over which TCWVret has been smoothed, e.g. for the 50 m×50 m to 100 m×100 m transition n=4. In this case, σε,ret must be calculated at the native resolution and therefore exploits the smaller TCWV variance at Δr∼50 m rather than the higher variance in a smoothed field with larger Δr.

3.2 Simulated retrieval results

3.2.1 TCWV retrievals over different surfaces

We first remind readers that “retrieval error” here only includes errors present in these synthetic retrievals and excludes several real-world sources, such as how the true atmosphere is not plane-parallel as assumed in our radiative transfer. Retrieved surface ρs spectra and TCWVret versus forward-model TCWV are shown in Fig. 4. Surface ρs are retrieved well, with mean bias magnitude equivalent to 0.2 %–1.6 % of true ρs (e.g. for Lambertian ρs=0.1, the mean bias is 0.00021), and standard deviation across all channels is 2 %–4 % of true ρs. The largest contribution to errors is from spikes near λ∼2.06µm. Inspection found that the MODTRAN CO2 concentration changes between default profiles versus prescribed T and q profiles. In future an up-to-date CO2 mixing ratio will be assigned, but the higher LUT value (361 ppmv) versus the forward-model value (323 ppmv) results in the retrieval overly brightening the surface in the strong CO2 band near λ∼2.06µm.

Figure 4(a)(c) Retrieved reflectance spectra for (a) vegetation, (b) mineral and (c) spectrally uniform surfaces. Lines show the mean of all simulated retrievals and shading shows ±1σ. (d)(f) Retrieved TCWV as a function of true TCWV for the same. The vegetation and mineral cases use three snapshots (N=303) and the Lambertian surfaces just use ARM_18000s (N=101).


Comparing Fig. 4d–f, TCWVret over mineral surfaces is a mean 0.44 mm higher than over vegetation. From panel f, some of this difference is likely related to mean surface brightness: darker surfaces give higher TCWVret. The other differences in TCWVret between surfaces must be due to spectral shape, but it appears that surface-induced errors are small when considering only mixed vegetation or mixed mineral surfaces. Regardless of the surface, a bias of order ∼1 mm remains, which is similar to the largest difference introduced by surface type and may be related to other retrieval errors such as inappropriate atmospheric profile shapes assumed in the LUT. However, the derived spatial statistics we are interested in here are not affected by any mean bias.

Figure 5 shows example scenes with different surface types. The true TCWV standard deviation σx is 0.28 mm (panel a), while over the uniform surfaces the retrieval gives 0.33 mm (panels b and c), with the larger value due mainly to the σε contribution. Over the striped surfaces it is 0.40 mm (panel d) due to the additional variance from combining surface types. However, if the top or bottom half of panel d is selected, then both return σx of 0.33 mm, i.e. the same as over a fixed vegetation or mineral surface. Statistics should not be taken over scenes with both vegetation and mineral surfaces, but the Isofit surface classification, which is output by the retrieval, should be used to identify areas of sufficiently similar surface type for calculation of TCWV spatial statistics. The rest of the analysis assumes the MODTRAN cropland default surface.

Figure 5ARM_18000s (a) true TCWV, (b) retrieved TCWV over a uniform vegetated surface, (c) retrieved TCWV over a uniform mineral surface, (d) retrieved TCWV over stripes of uniform surface types as labelled in the figure, (e) difference induced in retrieved TCWV by surface type relative to mixed vegetation as (d) minus (b). (f) Difference relative to a mixed mineral surface as (d) minus (c). Clouds are masked in all cases. Panel (f) represents (d) minus (c).


3.2.2 TCWV retrievals over vegetation surfaces in all LES snapshots

Figure 6 shows TCWV retrievals over the MODTRAN cropland and ocean surfaces. The poor performance over ocean absent sun glint justifies our land-only investigation. Over land the mean bias ranges from −3.0 % (DRY) to +1.8 % (BOMEX), while the within-scene σε is from 0.52 % (ARM_lsconv) to 0.67 % (BOMEX). As discussed in Sect. 3.1.3, VSWIR TCWV validation studies typically report error metrics larger than our σε, but their values include inter-product differences in bias, which are potentially far larger. Bias is indeed sensitive to the assumed meteorological profiles, since re-running the ARM_18000s retrievals using a LUT developed with the MODTRAN default “tropical” atmospheric profile shifts the mean bias from 0.33±0.04 mm to 0.14±0.04 mm (mean ±2σ).

For the purpose of spatial variability in TCWV at scales of tens of kilometres, the distinction between large-area and small-area retrieval errors is important. Generally speaking, the error in an individual column TCWV retrieval is of order 2 %–3 % since that includes the bias term, but for estimates of sub-10 km spatial variability, the within-LES 0.5 %–0.7 % is the error of interest.

3.2.3 Emulator parameters

Emulator parameters with ±2σ confidence intervals are listed in Table 2, and significant (p<0.05) non-unity trends can be seen most clearly for BOMEX (green) and RICO (purple) in Fig. 6a; the retrieved properties are more variable than reality, with trends of 1.34 and 1.22 mm mm−1, respectively. Meanwhile, the ARM and ARM_lsconv trends are both <1mmmm-1. Therefore σx calculated for BOMEX will be 34 % too high and for ARM 6 % too low.

Table 2Emulator parameters relating true TCWV to TCWVret. In Eq. (2) the trend is a1, the intercept is a2, and residual σ is the standard deviation used in generating the samples of ϵ. Values are shown ±2σ.

Download Print Version | Download XLSX

We argue that the most likely causes of emulated trend bias are related to the vertical T and q profile. Firstly, dI/dq is non-linear and varies with atmospheric conditions due to line broadening and interaction with aerosol layers. The a1 fit parameter may therefore be sensitive to differences between true profiles and those assumed in the retrieval LUT. Secondly, the LUT uniformly scales q(z) profiles, whereas the horizontal variability in q tends to peak at specific altitudes (Fig. 1).

Two tests provide some evidence for this. Firstly, when using different standard atmospheres to generate lookup tables for the DRY case, the retrieval gradient changes by 5 %, and secondly, when re-running all BOMEX retrievals with forward radiances generated using the same q profile shape that has been scaled to match the original range of TCWV, the retrieval gradient changes by 9 % (Supplement Fig. 12). These results suggest that retrievals could be improved by using more accurate meteorological profiles in the LUT development and by using a more appropriate scaling for q as a function of z in the LUT.

3.2.4 Snapshot statistics and estimation of random error

Figure 7a shows how σx of TCWVret is overestimated in every snapshot (circles). Figure 7b shows that the estimated retrieval error σε agrees excellently with the truth, and after removing σε, the triangles in Fig. 7a show the consistent overestimate is removed. Random error, such as that introduced by some instrumental uncertainties, can be precisely identified and removed from the spatial variance calculations. To estimate σx, the largest error source we consider is due to emulator slope. Other potential sources would be due to surface variation, which can be mitigated by selecting regions of similar surface classification as suggested in Sect. 4.1, and due to spatially varying errors, such as inter-pixel calibration biases or those induced by unmodelled temperature gradients across the sensors. Separate approaches are required to account for these issues.

Figure 6Retrieved TCWV as a function of the truth for all snapshots in each LES case over (a) cropland and (b) ocean. Note that the TCWV values differ from those derived from the LES due to differences in the MODTRAN layer interpolation and calculations.


Figure 7(a) Estimated clear-sky horizontal standard deviation as a function of the true value for each snapshot for raw retrievals (circles) and retrievals after removal of the random component of retrieval error, e.g. that induced by instrumental noise (triangles). (b) The estimate of retrieval error as in Sect. 3.1.4 as a function of the true error in each case.


Next, we consider the r2 coefficient between TCWV and TCWVret, with an illustration in Fig. 8, where the RICO_14400s TCWVret fields are shown at the native resolution and after smoothing down to 80 m. The random retrieval error is visible as speckling (Fig. 8a and b) and clearly reduces following smoothing. The 2D histograms (Fig. 8c and d) demonstrate the increase in r2 from 0.82 to 0.95 following a coarsening of the Δx resolution from 40 to 80 m.

Figure 8(a) Retrieved TCWV at 40 m resolution, (b) retrieved TCWV at 80 m resolution, (c) 2D histogram of retrieved TCWV as a function of the truth at 40 m resolution, and (d) 2D histogram of the same at 80 m resolution. The squared Pearson correlation coefficient, r2, is written in the upper left corner of (c, d).


Figure 9 summarises the true and estimated statistical values as horizontal resolution is sequentially degraded. Comparison of Fig. 9a and b reveals that there is only a small decrease in σx as resolution coarsens up to hundreds of metres, and the biases between estimated and true values follow emulator a1 trends as expected, with ARM and ARM_lsconv too low and DRY, RICO and BOMEX too high.

Figure 9Changes in the true and retrieved statistical properties for LES as a function of spatial resolution Δx. (a) The true standard deviation calculated directly from the LES output, (b) retrieved standard deviation after removing the estimated retrieval error as in Sect. 3.1.4, (c) r2 between true TCWV and TCWVret, and (d) estimated r2 using Eq. (7).


Regarding r2 in Fig. 9c and d, Eq. (7) reliably predicts true r2, so a user could determine the spatial resolution required to achieve a desired r2. In all snapshots r2>0.90 at 150 m resolution, and in 21 of 23 cases this is achieved at 100 m. Therefore, with the errors accounted for here, the EMIT instrument could capture 90 % of spatial variability in PCWVPBL at 100 m resolution in the PBL conditions examined here, a factor of 7.5 improvement in the MERIS full-resolution retrievals. However, this conclusion does not account for the spatial smearing caused by SZA.

3.3 Discussion of retrieval results and limitations

This section has addressed questions (2) and (3) from Sect. 1 and shown that random errors introduced by EMIT's instrumental error can be accurately identified and removed. Provided that an observed domain consists of mixed vegetation or mixed mineral surfaces, then our derived error in σx using EMIT is from −7 % to +34 %. Isofit returns surface type, meaning that such domains can be identified from retrievals.

Computational limitations forced adoption of an emulator approach, which provides a useful framework to assess error sources. Firstly, this framework shows that the errors of interest for retrieval of spatial statistics of PCWVPBL are the gradient a1, equivalent to dTCWVret/dTCWV, and random error σε. We show that σε can be estimated and removed and that the main error is that in a1, most likely driven by the retrieval's atmospheric profile assumptions, which can be addressed in future development. Our method to derive σε also allows users to predict spatial correlation; in particular, we found that an r2>0.9 requirement requires smoothing to 100–150 m resolution. This is a factor of 3–8 improvement in sampling relative to MERIS full resolution.

Limitations include the use of the same radiative transfer code for forward and inverse simulations, so spectroscopic errors were ignored, as were errors in cloud and shadow masking, those caused by topography, or errors that correlate between footprints.

Spectroscopy errors can be estimated (Thompson et al., 2020) and should shrink in future with developments, with ongoing research in water vapour absorption spectroscopy (Elsey et al., 2020; Lechevallier et al., 2018; Menang et al., 2021) and a history of targeted development of spectroscopy to improve retrievals, such as for OCO-2 (Drouin et al., 2016; O'Dell et al., 2018; Payne et al., 2020). The surface remote-sensing community has tools for addressing topography (Kobayashi and Sanga-Ngoie, 2008; Teillet et al., 1982), and there are also approaches to dealing with nearby clouds to minimise the effect of imperfect cloud edge identification, shadowing and 3D cloud-radiative effects (Massie et al., 2021). Nevertheless, these are all topics that are worth evaluating for Isofit-like TCWV retrievals.

We note that our σε is smaller than the errors reported in product intercomparison studies, but those studies implicitly capture variance due to differing mean biases, i.e. the a2 term in our emulator, which is larger than the other terms. An evaluation of our retrieved σx would require independent validated sources such as passive microwave or differential absorption lidar data with Δx≤50 m that are collocated with VSWIR TCWVret. Reported comparisons are typically of bias and root-mean-squared error (RMSE) of satellite VSWIR retrievals relative to surface-based or other satellite products and are calculated from datasets across a range of times and sometimes places. Furthermore, the comparison data generally have larger Δx and may not be perfectly collocated in time and space, introducing additional variance that contributes to reported RMSE. Typical published analyses include within their RMSE uncertainties due to differences in space and time of measurements and any differences between the a2 terms between the VSWIR and validation dataset retrievals. Therefore, these reported errors cannot be compared with our values, which are calculated within individual LES runs. We can, however, report that our errors are similar to Thompson et al. (2021)'s airborne Isofit retrieval statistics against nearby AERONET surface stations, which reported an RMSE of 2.8 mm. Flight C data from Fig. 9 of that paper show a spatial standard deviation of 0.19 mm when smoothed to Δx=48 m, which is within the LES-simulated range, and σε is estimated at 0.18 mm, although that is not comparable to our values since it uses AVIRIS-NG rather than EMIT and is over ocean sun glint rather than land.

Reported RMSEs over land for other VSWIR instruments include 0.9–1.3 mm for OCO-2 (Nelson et al., 2016), 1.4–3.7 mm for MERIS (Lindstrot et al., 2012), 0.9–2.0 mm for MODIS (Diedrich et al., 2015), 1.3–3.3 mm for OLCI (Preusker et al., 2021) and up to 2.4 mm for Sentinel-2 (Obregón et al., 2019). The range of TCWVret simulated in Fig. 4 is therefore consistent with typical errors reported for other instruments. Interestingly, Obregón et al. (2019) report a gradient of 0.9 between Sentinel-2 and AERONET TCWVret. This is derived from data across multiple sites and times and so cannot be compared to our gradients derived from individual LES cases but indicates that different retrievals may indeed have relationships between TCWV and TCWVret which are not 1:1, and thus our non-unity a1 parameters, which scale derived σx, are credible.

Figure 10ARM_lsconv_36000s-integrated water path calculated along (a) SZA=15 and (b) SZA=60 with the upward path directly up at zenith angle 0; values labelled TCWV in colour bar for simplicity. Panel (c) shows the difference for each footprint by subtracting the true TCWV at SZA=0 from panel (a), and panel (d) shows the same for subtracting the SZA=0 value from the SZA=60 value. The “cloud” mask in each case is now extended to include cloud shadows, and the illumination comes from the top of each panel; i.e. sunlight travelling down through the atmosphere has a component in the negative y direction.


4 Effect of SZA variation on retrieved properties

4.1 Calculation of TCWVret accounting for the light path at different solar zenith angles

Along-path-integrated water vapour (IWV) for SZA ranging from 0 to 60 inclusive in increments of 15 is calculated using ray tracing. The sunlight's horizontal component is in the negative y direction, viewing zenith angle is 0, and the ray is traced from the top of atmosphere to the centre of each surface footprint. Each partial grid cell encountered has its q weighted by the pressure-corrected path through that cell before obtaining IWV. The cloud mask is extended by a “shadow mask” where cloud LWP>1×10-3 mm along the solar direct ray path. This IWV is referred to as a TCWV for consistency with standard retrieval terminology, even though it is not directly a measure of the column over the footprint. The Sect. 3 analysis is then repeated using the same emulators developed using radiative transfer with SZA=45, a plane-parallel assumption and footprint column profiles. Different SZAs may change the sensitivities somewhat, but we do not expect results substantially outside the range of those considered here.

4.2 Effect of SZA variation on retrieved properties

Figure 10a and d show apparent TCWV in ARM_lsconv_36000s (i.e. when convection is most developed) at SZA=15 and 60, and in Fig. 10c and d the clear vertical pattern of positive followed by negative biases relative to true TCWV is clear, with greater magnitude and larger regions of continuous positive or negative bias at higher SZA.

Figure 112D histograms between clear-sky TCWV (true value integrated only in column over footprint) and the retrieved values at the corresponding footprint with SZA of (a) 15, (b) 30, (c) 45, and (d) 60. The r2 coefficient is in each panel, and the footprint resolution is the native output of Δx=50 m.


Figure 11 shows that this spatial smearing destroys the correspondence between footprint and path TCWV, with r2 around 0.1 with SZA as small as 30. This can be compensated only somewhat by spatial smoothing, as Fig. 12 shows that even footprints degraded to 300 m are affected by SZA. The calculated σx at SZA=0 match those from Fig. 9, with biases from the emulator slope parameter in Table 2. Larger SZA in these cases increases the magnitude of this bias, but the difference in σx as SZA changes from 15 to 60 is smaller than the RICO or BOMEX emulator-trend-induced biases. The retrieved σx with footprint size tracks reality, suggesting that the horizontal distribution statistics might still be captured even at large SZA. Furthermore, the statistical error estimation from Sect. 3.1.4 has effectively identical performance regardless of SZA (not shown).

However, Fig. 12b shows that a VSWIR-retrieved map TCWVret does not accurately represent the actual spatial variability in TCWV and by extension PCWVPBL even for SZA=15, and this is a fundamental limitation caused by the solar path through the atmosphere. In fact, the TCWVret map corresponds better to the TCWV map at the horizontal location where the downward solar path enters the PBL, but improvement in r2 is limited (Supplement Fig. 13).

Figure 12Clear-sky TCWV horizontal spatial statistics calculated for ARM_18000s (blue) and RICO_14400s (orange) as a function of the horizontal footprint size. (a) Standard deviation σx as in Fig. 9 and including the random error correction from Sect. 3.1.4. (b) Correlation coefficient between column true TCWV and that retrieved for the same footprint as SZA changes. Each line style represents a different SZA as labelled in the legend of (a).


Figure 13 shows that while the retrieved TCWV distributions are biased, as previously discussed, SZA increases cause only minor visible changes in distribution shape. This indicates that important statistics of the TCWV (and therefore PCWVPBL) field can be obtained at the native footprint resolution, despite the poor correspondence of any individual footprint to the column located at that position. The primary advantages of finer spatial resolution are that (i) it allows better calculation of σε than at coarser resolution using Sect. 3.1.4's method, due to the smaller Δr between footprints and, (ii) when calculating statistics such as standard deviation on local scales, statistical errors are reduced by the larger number of footprints. For example, Δx=50 m represents approximately 25 times more measurements than MODIS or MERIS. If standard deviation were desired for a 1 km×1 km region, N=16 from 250 m footprints results in a sampling error of ±17.7 % versus ±3.5 % for N=400 from 50 m footprints.

Figure 13Histograms of footprint-estimated clear-sky TCWV, with the truth shown in grey shading. The retrieval estimates are all scaled to remove the variance due to estimated random error. (a) Variation with SZA calculated at footprint size Δx=50 m and (b) variation with footprint size at SZA=0. In both panels the blue histograms are the same.


5 Discussion and conclusions

Modern and upcoming VSWIR instruments promise unprecedented horizontal resolution for the study of surface properties, with emphases ranging from mineral regions that are the source of dust (EMIT) to routine observation of agriculture and biodiversity (CHIME) to the full spectrum of study under the NASA 2017 Decadal Survey's Surface Biology and Geology (SBG) designated observable.

This study suggests potential synergies with the Decadal Survey's PBL targeted observable by showing that PCWVPBL variability at high spatial resolution can be inferred using the TCWVret that will be obtained from EMIT observations. While these measurements lack the vertical resolution that is necessary to advance PBL science, they provide a unique constraint on the mesoscale moisture variability and aggregation within the convective PBL. This analysis is restricted to daytime convective PBLs over land surfaces, which excludes deep convection but still represents a large fraction of meteorological conditions in the tropical to mid latitudes. Importantly, these are the precise conditions in which it is suspected that PBL moisture aggregation influences the timing of deep convective events. Furthermore, given the large number of scenes in which we expect to be able to derive these spatial statistics, these observations could prove useful for constraining the manner in which small-scale variability is parameterised in shallow convection or unified parameterisation schemes. The Isofit development team has curated additional spectra for a universal prior that includes cryosphere surfaces, but additional work may be necessary to evaluate TCWV over snow, since there is a snow absorption feature near λ=1µm whose depth depends on snow grain size (Painter et al., 2007) and which overlaps qv absorption. This may introduce surface–atmosphere covariance that affects the retrieval.

NASA's 2017 Decadal Survey encourages multi-instrument applications, and the VSWIR retrievals discussed here could be combined with radio occultation, thermal infrared (TIR) or passive microwave sounders, which have far larger horizontal resolution but obtain vertical profiles. Early explorations of joint VSWIR-TIR retrievals are promising, suggesting that the sensors provide complementary information on both atmospheric and surface properties. VSWIR could provide a prior constraint on TCWV in a collocated TIR retrieval, or the TIR-retrieved PCWVupper could be subtracted from VSWIR TCWV to estimate PCWVPBL, with VSWIR also providing the horizontal statistics of clear-sky PCWVPBL variability within the TIR footprint. Another opportunity is to use coincident TIR-retrieved profiles of T and q to either build a more appropriate LUT for the VSWIR retrieval or to select from among pre-computed LUTs.

In Isofit, the atmospheric component contributes a bias to dTCWVret/dTCWV and may be the largest source of our errors in σx, which range from −7 % to +34 % of true σx. Development allowing the use of prescribed profiles and the ability to assign variability in q to lower altitudes rather than uniform scaling at all altitudes should reduce these errors, as accounting for temperature reduced biases in MERIS TCWVret (Lindstrot et al., 2012).

This study also showed how SZA as small as 15 significantly degrades the accuracy of retrieved spatial patterns in TCWV, even at coarser resolutions similar to current sensors such as MERIS. However, the TCWV distribution was far less sensitive to SZA. While our results should strongly affect the interpretation of retrieved maps of TCWV from instruments like MERIS, they suggest that moments of the PCWVPBL distribution can be obtained at unprecedented horizontal resolution, which may be of use to developers of modern PBL schemes that use or assume such moments. We note that the LES TCWV distributions and their variation with spatial scale may not be realistic, since they tend to be overly dissipative on scales ≤6 grid cells (Bryan et al., 2003), but it is not clear that these biases affect our conclusion regarding the ability to obtain distributional statistics that represent horizontal variability at scales as small as 50 m.

Future work could address uncertainties that are ignored here, such as topography or cloud 3D radiative effects via 3D radiative transfer simulations which avoid several of our assumptions, such as a plane-parallel atmosphere. A particular limitation is that this analysis did not consider vertical structure or PBL height beyond using that derived from the LES mean profiles. In reality there may be errors in locally estimated PBL height, or that obtained from other sensors may be inconsistent with the max(dθ/dz) value used here, and targeted research on this topic would be worthwhile. Observational evaluation of these uncertainties could be performed using collocated airborne measurements of column water vapour from VSWIR and other instruments such as differential absorption lidar or passive microwave imagers, provided they can obtain sufficiently high spatial resolution. Finally, this work could be extended to other sensors, such as MSI on Sentinel-2, which is not hyperspectral but provides an exceptionally fine Δx of approximately 20 m. Additional high-resolution analysis may be required for this, since Fig. 9a and b imply increases in retrieved σx at Δx=20 m for the two simulations that were run at that resolution.

Despite these caveats, we have shown ways in which atmospheric correction outputs of surface property retrievals for EMIT can provide unique information on fine-scale PBL water vapour variability and also identified specific development tasks to improve the quality of its atmospheric outputs. With current tools it therefore seems likely that missions such as EMIT and CHIME, which are primarily designated as targeting surface observables, can provide unique information to the atmospheric science community.

Code availability

The Isofit retrieval package is available on GitHub (, last access: 2 June 2021) (, Brodrick et al., 2021) and MODTRAN from Spectral Sciences (, licence required, last access: 19 March 2020, Berk et al., 2014, 2015​​​​​​​).

Data availability

The surface models are either default MODTRAN or available from the Isofit GitHub under data/reflectance/surface_model_ucsb. The instrument noise model is from the Isofit GitHub under data/sbg_noise_coeffs.txt (, Brodrick et al., 2021). The LES output was generated using published large eddy simulation models, and the cases are described in the references in row (xi) of Table 1.


The supplement related to this article is available online at:

Author contributions

The authors contributed as follows: MTR (conceptualisation, investigation, methodology, visualisation, writing), DRT (methodology, writing and resources, namely Isofit), MJK (writing and resources, namely LES output), and MDL (conceptualisation, methodology, writing).

Competing interests

The authors declare that they have no conflict of interest.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The authors also thank Amin Nehrir and Brian Carroll for assistance with the HALO lidar data.

Financial support

This research has been supported by the National Aeronautics and Space Administration (grant no. NNN12AA01C).

Review statement

This paper was edited by Alexander Kokhanovsky and reviewed by two anonymous referees.


Bedka, K. M., Nehrir, A. R., Kavaya, M., Barton-Grimley, R., Beaubien, M., Carroll, B., Collins, J., Cooney, J., Emmitt, G. D., Greco, S., Kooi, S., Lee, T., Liu, Z., Rodier, S., and Skofronick-Jackson, G.: Airborne lidar observations of wind, water vapor, and aerosol profiles during the NASA Aeolus calibration and validation (Cal/Val) test flight campaign, Atmos. Meas. Tech., 14, 4305–4334,, 2021. 

Bennartz, R. and Fischer, J.: Retrieval of columnar water vapour over land from backscattered solar radiation using the Medium Resolution Imaging Spectrometer, Remote Sens. Environ., 78, 274–283,, 2001. 

Berk, A., Conforti, P., Kennett, R., Perkins, T., Hawes, F., and van den Bosch, J.: MODTRAN6: a major upgrade of the MODTRAN radiative transfer code, Velez-Reyes, M. and Kruse, F. A. (Eds.), Proceedings Volume 9088, Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XX, SPIE Defense + Security, 2014, Baltimore, Maryland, United States, p. 90880H,, 2014 (data available at:, last access: 19 March 2020). 

Berk, A., Conforti, P., and Hawes, F.: An accelerated line-by-line option for MODTRAN combining on-the-fly generation of line center absorption within 0.1 cm−1 bins and pre-computed line tails, edited by: Velez-Reyes, M. and Kruse, F. A., Proceedings Volume 9472, Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XXI, SPIE Defense + Security, 2015, Baltimore, Maryland, United States, p. 947217,, 2015 (data available at:, last access: 19 March 2020). 

Borger, C., Beirle, S., Dörner, S., Sihler, H., and Wagner, T.: Total column water vapour retrieval from S-5P/TROPOMI in the visible blue spectral range, Atmos. Meas. Tech., 13, 2751–2783,, 2020. 

Brodrick, P., Erickson, A., Fahlen, J., Olson, W., Thompson, D. R., Shiklomanov, A., Serbin, S. P., Carmon, N., and McGibbney, L. J.: Isofit 2.8.0, Zenodo [data set],, 2021. 

Brown, A. R., Cederwall, R. T., Chlond, A., Duynkerke, P. G., Golaz, J.-C., Khairoutdinov, M., Lewellen, D. C., Lock, A. P., MacVean, M. K., Moeng, C.-H., Neggers, R. A. J., Siebesma, A. P., and Stevens, B.: Large-eddy simulation of the diurnal cycle of shallow cumulus convection over land, Q. J. Roy. Meteor. Soc., 128, 1075–1093,, 2002. 

Bryan, G. H., Wyngaard, J. C., and Fritsch, J. M.: Resolution Requirements for the Simulation of Deep Moist Convection, Mon. Weather Rev., 131, 2394–2416,<2394:RRFTSO>2.0.CO;2, 2003. 

Candela, L., Formaro, R., Guarini, R., Loizzo, R., Longo, F., and Varacalli, G.: The PRISMA mission, in: 2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), IEEE, Beijing, China, 10–15 July 2016, pp. 253–256, 2016. 

Carbajal Henken, C. K., Diedrich, H., Preusker, R., and Fischer, J.: MERIS full-resolution total column water vapor: Observing horizontal convective rolls, Geophys. Res. Lett., 42, 10074–10081,, 2015. 

Couvreux, F., Guichard, F., Austin, P. H., and Chen, F.: Nature of the Mesoscale Boundary Layer Height and Water Vapor Variability Observed 14 June 2002 during the IHOP_2002 Campaign, Mon. Weather Rev., 137, 414–432,, 2009. 

Diedrich, H., Preusker, R., Lindstrot, R., and Fischer, J.: Retrieval of daytime total columnar water vapour from MODIS measurements over land surfaces, Atmos. Meas. Tech., 8, 823–836,, 2015. 

Drouin, B. J., Benner, D. C., Brown, L. R., Cich, M. J., Crawford, T. J., Devi, V. M., Guillaume, A., Hodges, J. T., Mlawer, E. J., Robichaud, D. J., Oyafuso, F., Payne, V. H., Sung, K., Wishnow, E. H., and Yu, S.: Multispectrum analysis of the oxygen A-band, J. Quant. Spectrosc. Ra., 186, 118–138,, 2016. 

Drusch, M., Del Bello, U., Carlier, S., Colin, O., Fernandez, V., Gascon, F., Hoersch, B., Isola, C., Laberinti, P., Martimort, P., Meygret, A., Spoto, F., Sy, O., Marchese, F., and Bargellini, P.: Sentinel-2: ESA's Optical High-Resolution Mission for GMES Operational Services, Remote Sens. Environ., 120, 25–36,, 2012. 

Elsey, J., Coleman, M. D., Gardiner, T. D., Menang, K. P., and Shine, K. P.: Atmospheric observations of the water vapour continuum in the near-infrared windows between 2500 and 6600 cm−1, Atmos. Meas. Tech., 13, 2335–2361,, 2020. 

Gao, B.-C. and Kaufman, Y. J.: Water vapor retrievals using Moderate Resolution Imaging Spectroradiometer (MODIS) near-infrared channels, J. Geophys. Res.-Atmos., 108, 4389,, 2003. 

Golaz, J.-C., Larson, V. E., and Cotton, W. R.: A PDF-Based Model for Boundary Layer Clouds. Part I: Method and Model Description, J. Atmos. Sci., 59, 3540–3551,<3540:APBMFB>2.0.CO;2, 2002. 

Green, R. O. and Thompson, D. R.: An Earth Science Imaging Spectroscopy Mission: The Earth Surface Mineral Dust Source Investigation (EMIT), in: IGARSS 2020–2020 IEEE International Geoscience and Remote Sensing Symposium, IEEE, Waikoloa, Hawaii, USA, 26 September–2 October 2020, pp. 6262–6265,, 2020. 

Grossi, M., Valks, P., Loyola, D., Aberle, B., Slijkhuis, S., Wagner, T., Beirle, S., and Lang, R.: Total column water vapour measurements from GOME-2 MetOp-A and MetOp-B, Atmos. Meas. Tech., 8, 1111–1133,, 2015. 

Guanter, L., Gómez-Chova, L., and Moreno, J.: Coupled retrieval of aerosol optical thickness, columnar water vapor and surface reflectance maps from ENVISAT/MERIS data over land, Remote Sens. Environ., 112, 2898–2913,, 2008. 

Kobayashi, S. and Sanga-Ngoie, K.: The integrated radiometric correction of optical remote sensing imageries, Int. J. Remote Sens., 29, 5957–5985,, 2008. 

Krutz, D., Müller, R., Knodt, U., Günther, B., Walter, I., Sebastian, I., Säuberlich, T., Reulke, R., Carmona, E., Eckardt, A., Venus, H., Fischer, C., Zender, B., Arloth, S., Lieder, M., Neidhardt, M., Grote, U., Schrandt, F., Gelmi, S., and Wojtkowiak, A.: The Instrument Design of the DLR Earth Sensing Imaging Spectrometer (DESIS), Sensors, 19, 1622,, 2019. 

Kurowski, M. J., Grabowski, W. W., Suselj, K., and Teixeira, J.: The Strong Impact of Weak Horizontal Convergence on Continental Shallow Convection, J. Atmos. Sci., 77, 3119–3137,, 2020. 

Kuze, A., Suto, H., Nakajima, M., and Hamazaki, T.: Thermal and near infrared sensor for carbon observation Fourier-transform spectrometer on the Greenhouse Gases Observing Satellite for greenhouse gases monitoring, Appl. Opt., 48, 6716,, 2009. 

Larson, V. E., Golaz, J.-C., and Cotton, W. R.: Small-Scale and Mesoscale Variability in Cloudy Boundary Layers: Joint Probability Density Functions, J. Atmos. Sci., 59, 3519–3539,<3519:SSAMVI>2.0.CO;2, 2002. 

Laszlo, I., Stamnes, K., Wiscombe, W. J., and Tsay, S.-C.: The Discrete Ordinate Algorithm, DISORT for Radiative Transfer, in Light Scattering Reviews, vol. 11, pp. 3–65, Springer, Berlin, Heidelberg, 2016. 

Lechevallier, L., Vasilchenko, S., Grilli, R., Mondelain, D., Romanini, D., and Campargue, A.: The water vapour self-continuum absorption in the infrared atmospheric windows: new laser measurements near 3.3 and 2.0 µm, Atmos. Meas. Tech., 11, 2159–2171,, 2018. 

Lee, C. M., Cable, M. L., Hook, S. J., Green, R. O., Ustin, S. L., Mandl, D. J., and Middleton, E. M.: An introduction to the NASA Hyperspectral InfraRed Imager (HyspIRI) mission and preparatory activities, Remote Sens. Environ., 167, 6–19,, 2015. 

Lindstrot, R., Preusker, R., Diedrich, H., Doppler, L., Bennartz, R., and Fischer, J.: 1D-Var retrieval of daytime total columnar water vapour from MERIS measurements, Atmos. Meas. Tech., 5, 631–646,, 2012. 

Massie, S. T., Cronk, H., Merrelli, A., O'Dell, C., Schmidt, K. S., Chen, H., and Baker, D.: Analysis of 3D cloud effects in OCO-2 XCO2 retrievals, Atmos. Meas. Tech., 14, 1475–1499,, 2021. 

Matheou, G. and Chung, D.: Large-Eddy Simulation of Stratified Turbulence. Part II: Application of the Stretched-Vortex Model to the Atmospheric Boundary Layer, J. Atmos. Sci., 71, 4439–4460,, 2014. 

Menang, K. P., Gbode, I. E., and Adeyeri, O. E.: The effect of the differences in near-infrared water vapour continuum models on the absorption of solar radiation, Meteorol. Atmos. Phys.,, 2021. 

Millán, L., Lebsock, M., Fishbein, E., Kalmus, P., and Teixeira, J.: Quantifying Marine Boundary Layer Water Vapor beneath Low Clouds with Near-Infrared and Microwave Imagery, J. Appl. Meteorol. Clim., 55, 213–225,, 2016. 

National Academies of Science, Engineering, and Medicine: Thriving on Our Changing Planet, National Academies Press, Washington, DC, 2018. 

Nelson, R. R., Crisp, D., Ott, L. E., and O'Dell, C. W.: High-accuracy measurements of total column water vapor from the Orbiting Carbon Observatory-2, Geophys. Res. Lett., 43, 12261–12269,, 2016. 

Noël, S., Buchwitz, M., Bovensmann, H., Hoogen, R., and Burrows, J. P.: Atmospheric water vapor amounts retrieved from GOME satellite data, Geophys. Res. Lett., 26, 1841–1844,, 1999. 

Noël, S., Buchwitz, M., and Burrows, J. P.: First retrieval of global water vapour column amounts from SCIAMACHY measurements, Atmos. Chem. Phys., 4, 111–125,, 2004. 

Obregón, M. Á., Rodrigues, G., Costa, M. J., Potes, M., and Silva, A. M.: Validation of ESA Sentinel-2 L2 A Aerosol Optical Thickness and Columnar Water Vapour during 2017–2018, Remote Sens., 11, 1649,, 2019. 

O'Dell, C. W., Eldering, A., Wennberg, P. O., Crisp, D., Gunson, M. R., Fisher, B., Frankenberg, C., Kiel, M., Lindqvist, H., Mandrake, L., Merrelli, A., Natraj, V., Nelson, R. R., Osterman, G. B., Payne, V. H., Taylor, T. E., Wunch, D., Drouin, B. J., Oyafuso, F., Chang, A., McDuffie, J., Smyth, M., Baker, D. F., Basu, S., Chevallier, F., Crowell, S. M. R., Feng, L., Palmer, P. I., Dubey, M., García, O. E., Griffith, D. W. T., Hase, F., Iraci, L. T., Kivi, R., Morino, I., Notholt, J., Ohyama, H., Petri, C., Roehl, C. M., Sha, M. K., Strong, K., Sussmann, R., Te, Y., Uchino, O., and Velazco, V. A.: Improved retrievals of carbon dioxide from Orbiting Carbon Observatory-2 with the version 8 ACOS algorithm, Atmos. Meas. Tech., 11, 6539–6576,, 2018. 

Painter, T. H., Molotch, N. P., Cassidy, M., Flanner, M., and Steffen, K.: Contact spectroscopy for determination of stratigraphy of snow optical grain size, J. Glaciol., 53, 121–127,, 2007. 

Payne, V. H., Drouin, B. J., Oyafuso, F., Kuai, L., Fisher, B. M., Sung, K., Nemchick, D., Crawford, T. J., Smyth, M., Crisp, D., Adkins, E., Hodges, J. T., Long, D. A., Mlawer, E. J., Merrelli, A., Lunny, E., and O'Dell, C. W.: Absorption coefficient (ABSCO) tables for the Orbiting Carbon Observatories: Version 5.1, J. Quant. Spectrosc. Ra., 255, 107217,, 2020. 

Pérez-Ramírez, D., Whiteman, D. N., Smirnov, A., Lyamani, H., Holben, B. N., Pinker, R., Andrade, M., and Alados-Arboledas, L.: Evaluation of AERONET precipitable water vapor versus microwave radiometry, GPS, and radiosondes at ARM sites, J. Geophys. Res.-Atmos., 119, 9596–9613,, 2014. 

Preusker, R., Carbajal Henken, C., and Fischer, J.: Retrieval of Daytime Total Column Water Vapour from OLCI Measurements over Land Surfaces, Remote Sens., 13, 932,, 2021. 

Prusa, J. M., Smolarkiewicz, P. K., and Wyszogrodzki, A. A.: EULAG, a computational model for multiscale flows, Comput. Fluids, 37, 1193–1207,, 2008. 

Rast, M., Ananasso, C., Bach, H., Ben-Dor, E., Chabrillat, S., Colombo, R., Del Bello, U., Feret, J., Giardino, C., Green, R., Guanter, L., Marsh, S., Nieke, J., Ong, C. C. H., Rum, G., Schaepman, M., Schlerf, M., Skidmore, A., and Strobl, P.: Copernicus hyperspectral imaging mission for the environment: Mission requirements document, v. 2.1, ESA/ESTEC, Noordwijk, the Netherlands, 2019. 

Rothman, L. S., Gordon, I. E., Barbe, A., Benner, D. C., Bernath, P. F., Birk, M., Boudon, V., Brown, L. R., Campargue, A., Champion, J.-P., Chance, K., Coudert, L. H., Dana, V., Devi, V. M., Fally, S., Flaud, J.-M., Gamache, R. R., Goldman, A., Jacquemart, D., Kleiner, I., Lacome, N., Lafferty, W. J., Mandin, J.-Y., Massie, S. T., Mikhailenko, S. N., Miller, C. E., Moazzen-Ahmadi, N., Naumenko, O. V., Nikitin, A. V., Orphal, J., Perevalov, V. I., Perrin, A., Predoi-Cross, A., Rinsland, C. P., Rotger, M., Šimečková, M., Smith, M. A. H., Sung, K., Tashkun, S. A., Tennyson, J., Toth, R. A., Vandaele, A. C., and Vander Auwera, J.: The HITRAN 2008 molecular spectroscopic database, J. Quant. Spectrosc. Ra., 110, 533–572,, 2009. 

Roy, R. J., Lebsock, M., Millán, L., Dengler, R., Rodriguez Monje, R., Siles, J. V., and Cooper, K. B.: Boundary-layer water vapor profiling using differential absorption radar, Atmos. Meas. Tech., 11, 6511–6523,, 2018. 

Roy, R. J., Lebsock, M., Millán, L., and Cooper, K. B.: Validation of a G-Band Differential Absorption Cloud Radar for Humidity Remote Sensing, J. Atmos. Ocean. Techn., 37, 1085–1102,, 2020. 

Schaepman-Strub, G., Schaepman, M. E., Painter, T. H., Dangel, S., and Martonchik, J. V.: Reflectance quantities in optical remote sensing–definitions and case studies, Remote Sens. Environ., 103, 27–42,, 2006. 

Schneider, A., Borsdorff, T., aan de Brugh, J., Aemisegger, F., Feist, D. G., Kivi, R., Hase, F., Schneider, M., and Landgraf, J.: First data set of H2O/HDO columns from the Tropospheric Monitoring Instrument (TROPOMI), Atmos. Meas. Tech., 13, 85–100,, 2020. 

Siebesma, A. P., Bretherton, C. S., Brown, A., Chlond, A., Cuxart, J., Duynkerke, P. G., Jiang, H., Khairoutdinov, M., Lewellen, D., Moeng, C.-H., Sanchez, E., Stevens, B., and Stevens, D. E.: A Large Eddy Simulation Intercomparison Study of Shallow Cumulus Convection, J. Atmos. Sci., 60, 1201–1219,<1201:ALESIS>2.0.CO;2, 2003. 

Sissenwine, N., Dubin, M., and Teweles, S.: US Standard Atmosphere, National Oceanographic and Atmospheric Administration, Washington, DC, 1976. 

Stamnes, K., Tsay, S.-C., Wiscombe, W., and Jayaweera, K.: Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media, Appl. Opt., 27, 2502,, 1988. 

Stirling, A. J. and Petch, J. C.: The impacts of spatial variability on the development of convection, Q. J. Roy. Meteor. Soc., 130, 3189–3206,, 2004. 

Suselj, K., Kurowski, M. J., and Teixeira, J.: A Unified Eddy-Diffusivity/Mass-Flux Approach for Modeling Atmospheric Convection, J. Atmos. Sci., 76, 2505–2537,, 2019. 

Szczodrak, M., Austin, P. H., and Krummel, P. B.: Variability of Optical Depth and Effective Radius in Marine Stratocumulus Clouds, J. Atmos. Sci., 58, 2912–2926,<2912:VOODAE>2.0.CO;2, 2001. 

Teillet, P. M., Guindon, B., and Goodenough, D. G.: On the Slope-Aspect Correction of Multispectral Scanner Data, Can. J. Remote Sens., 8, 84–106,, 1982. 

Thompson, D. R., Natraj, V., Green, R. O., Helmlinger, M. C., Gao, B.-C., and Eastwood, M. L.: Optimal estimation for imaging spectrometer atmospheric correction, Remote Sens. Environ., 216, 355–373,, 2018. 

Thompson, D. R., Cawse-Nicholson, K., Erickson, Z., Fichot, C. G., Frankenberg, C., Gao, B.-C., Gierach, M. M., Green, R. O., Jensen, D., Natraj, V., and Thompson, A.: A unified approach to estimate land and water reflectances with uncertainties for coastal imaging spectroscopy, Remote Sens. Environ., 231, 111198,, 2019. 

Thompson, D. R., Braverman, A., Brodrick, P. G., Candela, A., Carmon, N., Clark, R. N., Connelly, D., Green, R. O., Kokaly, R. F., Li, L., Mahowald, N., Miller, R. L., Okin, G. S., Painter, T. H., Swayze, G. A., Turmon, M., Susilouto, J., and Wettergreen, D. S.: Quantifying uncertainty for remote spectroscopy of surface composition, Remote Sens. Environ., 247, 111898,, 2020. 

Thompson, D. R., Kahn, B. H., Brodrick, P. G., Lebsock, M. D., Richardson, M., and Green, R. O.: Spectroscopic imaging of sub-kilometer spatial structure in lower-tropospheric water vapor, Atmos. Meas. Tech., 14, 2827–2840,, 2021. 

Trent, T., Boesch, H., Somkuti, P., and Scott, N.: Observing Water Vapour in the Planetary Boundary Layer from the Short-Wave Infrared, Remote Sens., 10, 1469,, 2018.  

vanZanten, M. C., Stevens, B., Nuijens, L., Siebesma, A. P., Ackerman, A. S., Burnet, F., Cheng, A., Couvreux, F., Jiang, H., Khairoutdinov, M., Kogan, Y., Lewellen, D. C., Mechem, D., Nakamura, K., Noda, A., Shipway, B. J., Slawinska, J., Wang, S., and Wyszogrodzki, A.: Controls on precipitation and cloudiness in simulations of trade-wind cumulus as observed during RICO, J. Adv. Model. Earth Sy., 3, M06001,, 2011. 

Vermote, E. F., El Saleous, N., Justice, C. O., Kaufman, Y. J., Privette, J. L., Remer, L., Roger, J. C., and Tanré, D.: Atmospheric correction of visible to middle-infrared EOS-MODIS data over land surfaces: Background, operational algorithm and validation, J. Geophys. Res.-Atmos., 102, 17131–17141,, 1997. 

von Engeln, A. and Teixeira, J.: A Planetary Boundary Layer Height Climatology Derived from ECMWF Reanalysis Data, J. Climate, 26, 6575–6590,, 2013. 

Wulfmeyer, V., Bauer, H.-S., Grzeschik, M., Behrendt, A., Vandenberghe, F., Browell, E. V., Ismail, S., and Ferrare, R. A.: Four-Dimensional Variational Assimilation of Water Vapor Differential Absorption Lidar Data: The First Case Study within IHOP_2002, Mon. Weather Rev., 134, 209–230,, 2006. 

Short summary
Modern and upcoming hyperspectral imagers will take images with spatial resolutions as fine as 20 m. They can retrieve column water vapour, and we show evidence that from these column measurements you can get statistics of planetary boundary layer (PBL) water vapour. This is important information for climate models that need to account for sub-grid mixing of water vapour near the surface in their PBL schemes.