Spectroscopic imaging of sub-kilometer spatial structure in lower-tropospheric water vapor

Abstract

The subgrid spatial variability of water vapor is an important geophysical parameter for modeling tropical convention and cloud processes in atmospheric models. This study maps sub-kilometer spatial structures in total atmospheric column water vapor with visible to shortwave infrared (VSWIR) imaging spectroscopy. We describe our inversion approach and validate its accuracy with coincident measurements by airborne imaging spectrometers and the AERONET ground-based observation network. Next, data from NASA's AVIRIS-NG spectrometer enable the highest-resolution measurement to date of water vapor's spatial variability and scaling properties. We find second-order structure function scaling exponents consistent with prior studies of convective atmospheres. Airborne lidar data show that this total column measurement provides information about variability in the lower troposphere. We conclude by discussing the implications of these measurements and paths toward future campaigns to build upon these results.

1 Introduction

The complex spatial distribution of atmospheric water vapor surrounding clouds and precipitation structures has important consequences for parameterizing moist processes in atmospheric models. At the scale of general circulation models (GCMs), water vapor plays an important role in tropical moist convection and its associated precipitation . The mean and variability of precipitation rate in the tropics are strongly dependent on the atmospheric water vapor , a fact which has implications for parameterizing convection. Another ubiquitous property of convection is its tendency to aggregate . There is evidence the degree of aggregation will change as the climate warms, potentially changing the cloud feedback . Models and observations suggest that the tendency of convection to aggregate depends on the degree of spatial variance in the water vapor field. Over land surfaces with heterogeneous surface conditions, the variability in atmospheric water vapor can be larger and is seen as a critical component of the timing of deep convection . These variations in water vapor over convective continental environments are primarily driven by variability below 2 km altitude and within the planetary boundary layer (PBL) . Accurate water vapor parameterization is also important for cloud-resolving or convection-permitting models operating at kilometer scales and large-eddy simulations at sub-kilometer resolution. Across all scales, water vapor variability, and its coupling to cloud types and multi-scale organization, is key for advancing the parameterization and simulation of cloud processes.

To the end of advancing remote observations of atmospheric water vapor, this paper focuses on a specific measurement that is independently useful and also typifies the more general challenge of observing variability. Structure functions measure the change in the vapor field as a function of distance, quantifying its spatial texture across scales. They are used to analyze atmospheric components including temperature, winds, and trace gas concentrations, at scales from tens to hundreds of kilometers . Water vapor structure functions can indicate the recent history of the air mass. They can distinguish convective and non-convective systems and suggest precipitation rates embedded within a moist column of air with variable levels of column water vapor . For a one-dimensional field $f\left(i\right)$ indexed by location $i$, the $n$th-order structure function $S_n\left(r\right)$ is

1$$S_n\left(r\right)=E\left[\right|f(i+r)-f\left(i\right){|}^n],$$ where $r$ is a separation distance between pairs of points, and $E$ is the expectation over locations. More generally, $S_n\left(r\right)$ can represent variability along one direction of a multidimensional field. At least one study of water vapor data has found minor differences in satellite cross-track and along-track directions . However, it is more typical to assume that water vapor scaling has no preferred horizontal orientation, and that the structure function is rotationally symmetric. $S_n\left(r\right)$ is estimated using the mean of observed water vapor values at different spatial offsets. Over a restricted range of distances, structure functions can be described with a power law: 2$$S_n\left(r\right)\propto r^{{\mathit{\zeta }}_n},$$ where ${\mathit{\zeta }}_n$ is the scaling exponent of the order of $n$. The scaling exponent of the order of 2 is related to the commonly used Fourier power spectrum exponent $\mathit{\beta }$: 3$$\mathit{\beta }=-({\mathit{\zeta }}_{\mathrm{2}}+\mathrm{1}).$$ These values can diagnose specific atmospheric transport processes. Following Kolmogorov theory, a passive tracer in turbulence has a theoretical second-order scaling exponent ${\mathit{\zeta }}_{\mathrm{2}}$ of $\mathrm{2}/\mathrm{3}$, or equivalently a Fourier power spectrum exponent $\mathit{\beta }$ of $-\mathrm{5}/\mathrm{3}$ .

Previous studies have used a range of instruments to measure water vapor structure functions. In situ aircraft sensors measure a one-dimensional time series along the aircraft flight track. measured separation distances from 150 to 1500 km using this technique. They found ${\mathit{\zeta }}_{\mathrm{2}}$ ranging from the theoretical value of $\mathrm{2}/\mathrm{3}$ to unity. Later found ${\mathit{\zeta }}_{\mathrm{2}}$ consistent with $\mathrm{2}/\mathrm{3}$ over distances of several kilometers. These surveys have recently been augmented by airborne lidar measurements, which provide vertical profiles along the flight track. A series of lidar observations by measured air masses at scales down to 2 km. They reported ${\mathit{\zeta }}_{\mathrm{2}}=\mathrm{0.6}-\mathrm{0.75}$ in convective air masses and ${\mathit{\zeta }}_{\mathrm{2}}=\mathrm{1.0}-\mathrm{1.2}$ in non-convective air masses. Convective environments had shallower scaling indicating higher spatial variability at small scales. They hypothesized that this flatter power spectrum slope, which was also common in boundary layer air masses, was related to the small-scale injection of water vapor anomalies by convective eddies. However, their study did not conclusively identify the cause. More recently, compared airborne lidar measurements with simulations at scales greater than 11 km. This study confirmed that power-law exponents were strongly related to altitude and the presence or absence of convection. They showed ${\mathit{\zeta }}_{\mathrm{2}}$ above 1 for non-convective air masses in the free troposphere, near 0.6 in the boundary layer, and as low as 0.2 in convective air masses.

Contrasting with localized, high-resolution aircraft data, orbital data such as the Atmospheric Infrared Sounder (AIRS) on EOS Aqua have provided more comprehensive power-law exponents across the mesoscale and synoptic scales but without comparing them to GCMs. developed a global climatology of scaling properties between 150–1200 km based on vertically resolved temperature ($T$) and specific humidity ($q$) observations by AIRS. extended these results and compared them to climate GCMs, MERRA reanalysis, and VOCALS-REx in situ observations of vertically resolved $T$ and $q$ within and above shallow cumulus off of the coast of Peru. They found strong evidence of ${\mathit{\zeta }}_{\mathrm{2}}=\mathrm{2}/\mathrm{3}$ scaling in the boundary layer and tropopause at all latitudes, with steeper scaling in the mid-troposphere and at low latitudes. This was consistent with the view of water vapor as a passive scalar in turbulent flow and implied more small-scale variability than was predicted by contemporaneous GCMs. This increased scaling of the tropical free troposphere was also consistent with theoretical predictions of steeper scaling of column water vapor within and near strongly precipitating convection . Interestingly, their analysis of aircraft data from the VOCALS-REx experiment suggested a steepening of the curve at the smallest scales below 10 km . This result suggested a possible change in the scaling properties at the finest resolutions.

These studies contribute to a growing body of literature on water vapor scaling. However, important gaps remain. High-spatial-resolution data are sparse, and most studies explore spatial scales above 1 km. To our knowledge no study has yet corroborated the steepening phenomenon at scales less than 10 km or probed the structure functions at scales less than 2 km. Most aircraft data consist of one-dimensional time series, rather than the two-dimensional maps available from instruments like AIRS. Measurements confined to the flight trajectory provide fewer samples, increasing uncertainty in the derived exponents .

A new generation of orbital instrumentation may shed new light on fine-scale water vapor. Visible to shortwave infrared (VSWIR) imaging spectrometers, such as NASA's upcoming EMIT mission or the Surface Biology and Geology (SBG) investigation , are highly sensitive to the water vapor column absorption . They typically have a spectral resolution of 5–10 nm and span the 380–2500 nm interval, a range which overlaps significant water absorption features. A typical push-broom instrument could have as many as 1200 cross-track measurements with ground sampling distances of 30–60 m. Most such investigations target surface properties. However, a byproduct of this analysis will be accurate column-averaged water vapor measurements at high spatial resolution over wide areas. The spatial resolution and accuracy of the derived water vapor column measurements will be unprecedented, providing opportunities to probe the horizontal variability of water vapor on global scales. Before these missions launch, archives of airborne precursor data provide an opportunity to validate the technique and begin the investigation.

This paper demonstrates direct mapping of the sub-kilometer spatial structure in column water vapor using an airborne VSWIR imaging spectrometer. We first describe our model for vapor absorption features in the near to shortwave infrared. We estimate clear-sky vapor concentrations by inverting a combined model of the atmosphere and surface reflectance after filtering out cloud-affected footprints. We validate the approach with NASA's “Classic” Airborne Visible Infrared Imaging Spectrometer, AVIRIS-C , evaluating overflights of the AERONET observation network . We analyze imaging spectroscopy data acquired during the “Next Generation” Airborne Visible Infrared Imaging Spectrometer (AVIRIS-NG) India campaign of 2018. This campaign includes several scenes with highly favorable solar geometry, providing a uniquely high-resolution measurement of the H${}_{\mathrm{2}}$O vapor column. This enables estimation of the H${}_{\mathrm{2}}$O spatial structure functions at sub-kilometer scales. We find confirming evidence of ${\mathit{\zeta }}_{\mathrm{2}}=\mathrm{2}/\mathrm{3}$ scaling in some but not all atmospheres, with curve steepening that continues down to at least 100 m. Finally, we discuss the relationship between the total column measurement and variability in the lower troposphere. We assess water vapor profiles from airborne lidar campaigns . These data show that lower-tropospheric variability consistently dominates the total column, making the VSWIR measurement informative about lower-atmospheric water vapor. We conclude by discussing the implications of these measurements and future campaigns that build upon the results.

2 Methods

2.1 Atmospheric model

Our technique estimates water vapor independently for each spatial location by inverting an atmospheric radiative transfer model. We define a state vector $\mathit{x}$ containing all the free parameters in the system. It includes the surface reflectance vector $\mathit{\rho }$, the column water vapor concentration $q_{\mathrm{v}}$, and the aerosol optical depth at 550 nm. Aerosol optical properties were derived from a canonical sulfate aerosol type and validated in prior studies of the India campaign . A forward model $\mathit{F}\mathbf{\left(}\mathit{x}\mathbf{\right)}$ maps this state onto an observed radiance at the sensor, ${\mathit{L}}_o$:

4$${\mathit{L}}_o=\mathit{F}\mathbf{\left(}\mathit{x}\mathbf{\right)}+\mathit{ϵ}.$$ Boldface represents a vector or matrix; e.g., ${\mathit{L}}_o$ has one element for each spectrometer channel. The random variable $\mathit{ϵ}$ represents instrument noise, distributed according to a zero-mean Gaussian with covariance ${\mathbf{\Sigma }}_e$. We can decompose ${\mathit{L}}_o$ into different photon paths, with $\odot$ representing element-wise multiplication: 5$${\mathit{L}}_o={\mathit{L}}_{\mathrm{atm}}+\mathit{t}\odot {\mathit{L}}_{\mathrm{dn}}\odot {\displaystyle \frac{\mathit{\rho }}{\mathrm{1}-\mathit{S}\odot \mathit{\rho }}}+\mathit{ϵ}.$$ Here ${\mathit{L}}_{\mathrm{atm}}$ represents the path radiance caused by molecular and particle scattering; these photons never reach the surface. The second term represents all the photons that interact with the surface at least once. ${\mathit{L}}_{\mathrm{dn}}$ is the downwelling illumination at the bottom of the atmosphere. $\mathit{\rho }$ is the spectral surface reflectance. $\mathit{t}$ is the atmospheric transmittance from the surface to the sensor along its line of sight. $\mathit{S}$ is the spherical sky albedo observed from the ground. The relation in Eq. () holds for a locally homogeneous and Lambertian surface, but small departures from these conditions are not catastrophic, and in any case the assumptions hold sufficiently well for the scenes in this study.

Our atmospheric model calculates the optical coefficients ${\mathit{L}}_{\mathrm{atm}}$, $\mathit{t}$, and $\mathit{S}$ using the MODTRAN 6.0 software package . Specifically, we use the discrete ordinate radiative transfer (DISORT) code with eight-stream multiple-scattering calculations. The absorption model uses a correlated-$k$ representation with 0.1 cm${}^{-\mathrm{1}}$ bins. Vertical profiles are assigned from a 20-layer stratified atmosphere. For computational efficiency, we do not run the complete radiative transfer model (RTM) for each evaluation of $\mathit{F}\mathbf{\left(}\mathit{x}\mathbf{\right)}$. Instead, we calculate optical coefficients in advance to fill a lookup table for each component in the right side of Eq. (), and then we interpolate within this table to determine the precise radiance for any given state vector. Figure  illustrates the atmospheric transmittance in the measurement interval from 380–2500 nm, including several prominent absorption features related to H${}_{\mathrm{2}}$O rovibrational overtones. Absorption features at 940 and 1140 nm carry most of the water column information. Stronger features at 1380 and 1880 nm are saturated at atmospheric path lengths and consequently less useful.

Figure 1

Atmospheric transmissivity from 380–2500 nm reveals multiple H${}_{\mathrm{2}}$O absorption features, indicated with red arrows.

[Figure omitted. See PDF]

The remote H${}_{\mathrm{2}}$O measurements aim to quantify the integrated mass of water vapor directly above each pixel in the scene. However, this is not always possible for a solar-reflected signal. Neglecting scattering, which is small in the near infrared, the remote water vapor observation measures absorption along a two-part optical path from the sun to the ground to the sensor. The sun is seldom directly overhead, and its downwelling illumination enters the atmosphere at some horizontal offset from the reflection point (Fig. ). Consequently, a spectrum's spatial sensitivity can extend far beyond the target pixel in the direction of the sun. Projected onto the ground, the vapor absorption path forms a long, thin footprint that extends tens or hundreds of meters from the target pixel in the sunward direction. To quantify this effect, we define the water vapor spatial sensitivity as the extinction-weighted distance between downwelling and upwelling beams. We use $g\left(h\right)$ to represent the relative sensitivity to water vapor at $h$, some horizontal offset distance in the solar direction. For a nadir-pointed observation, assuming that the vertical profile of water vapor is locally constant, the relationship is 6$$g\left(h\right)=\left\{\begin{array}{cc}\underset{z\in [z_{\mathrm{s}},z_o]}{\int }\mathit{\kappa }\left(z\right)\mathrm{d}z& \text{if}h=\mathrm{0}\\ & \\ \mathit{\kappa }(\frac{h}{\tan {\mathit{\theta }}_{\mathrm{s}}}+z_{\mathrm{s}})\frac{\mathrm{1}}{\cos {\mathit{\theta }}_{\mathrm{s}}}& \text{if}h>\mathrm{0}\end{array},\right.$$ where $\mathit{\kappa }\left(z\right)$ is the water vapor number density at altitude $z$. The variables $z_{\mathrm{s}}$ and $z_o$ represent the surface elevation and sensor altitude, respectively. The spatial sensitivity is strongly dependent on solar zenith angle (${\mathit{\theta }}_{\mathrm{s}}$). It is also dependent to some extent on the view angle; if the observation is off-nadir, $g\left(h\right)$ includes the sensitivity of both upward and downward paths. We define the average spatial offset of a measurement, ${\mathit{\mu }}_{\mathrm{h}}$, as the mean horizontal position of the water vapor along the two-part optical path: 7$${\mathit{\mu }}_{\mathrm{h}}={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}hg\left(h\right)\mathrm{d}h}{{\int }_{\mathrm{0}}^{\mathrm{\infty }}g\left(h\right)\mathrm{d}h}}.$$ In addition to shifting the center of the spatial response, larger solar zenith angles can coarsen the effective spatial resolution by stretching the water vapor sensitivity footprint in the sunward direction. We define the effective spatial resolution, $R_{\mathrm{e}}$, as the symmetric horizontal distance from the center point which encloses 68.2 % of the spatial response. This would be equivalent to the area inside a single standard deviation if the response function were Gaussian. 8$$\mathrm{0.682}={\displaystyle \frac{{\int }_{{\mathit{\mu }}_{\mathrm{h}}-R_{\mathrm{e}}}^{{\mathit{\mu }}_{\mathrm{h}}+R_{\mathrm{e}}}g\left(h\right)\mathrm{d}h}{{\int }_{\mathrm{0}}^{\mathrm{\infty }}g\left(h\right)\mathrm{d}h}}$$ The left panel of Fig.  shows the relative sensitivity to water vapor at different horizontal offsets for a nadir-viewing measurement of a tropical atmosphere at 4 km acquisition altitude, similar to the observing geometry of the following experiments. We show the response areas for three solar zenith angles. As the solar angles increase, the response function extends farther from the target pixel, with increasingly “thick tails” caused by water absorption along the slanted downwelling path. Figure  shows the resulting resolutions as a function of solar zenith angle at different viewing geometries and atmospheric profiles, at sensor altitudes of 4 km and low Earth orbit (LEO). It shows that the spatial sensitivity is only weakly dependent on the sensor altitude and somewhat dependent on the viewing angle and the vertical profile of water vapor in the atmosphere. The solar zenith angle is the dominant influence on spatial resolution, motivating a careful selection of flight lines for the following experiments. The spatial footprint projected on the ground is not radially symmetric; it is long and thin, but retains the native spatial resolution along its short axis. Consequently, for isotropic structure functions the effective spatial resolution is the worst case, and sensitivity improves as one calculates structure functions in directions orthogonal to the sun.

Figure 2

The spatial sensitivity of solar-reflected measurements depends on water vapor absorption along the optical path from the sun to the ground to the sensor. Aircraft image credit: NASA.

[Figure omitted. See PDF]

Figure 3

The relative sensitivity to water vapor at 4 km acquisition altitude in a tropical atmosphere as a function of horizontal distance from nadir, partitioned by solar zenith angle ${\mathit{\theta }}_{\mathrm{s}}$.

[Figure omitted. See PDF]

Figure 4

The effective spatial resolution as a function of solar zenith angle ${\mathit{\theta }}_{\mathrm{s}}$, at both 4 km and low-Earth-orbit (LEO) acquisition altitudes. The red line shows a sub-arctic winter atmospheric profile indicating the range of spatial resolution across two very different atmospheric conditions. In reality, the sun never reaches low zenith angles in polar regions of the globe.

[Figure omitted. See PDF]

2.2 Inversion methodology

To estimate $q_{\mathrm{v}}$, we invert a combined model of surface and atmosphere. We use a Bayesian maximum a posteriori formalism that is common among atmospheric sounding missions. Recent work extended this approach to the visible–shortwave infrared spectral interval . The solar-reflected regime is strongly influenced by variability in surface reflectance. Consequently, we fit atmospheric parameters simultaneously with a flexible surface model. As described in Sect. , our state vector $\mathit{x}$ includes surface reflectance in every channel, the column water vapor concentration $q_{\mathrm{v}}$, and the aerosol optical depth at 550 nm. A forward model $\mathit{F}\mathbf{\left(}\mathit{x}\mathbf{\right)}$ transforms this state vector into a simulated radiance at the sensor following Eq. ().

Our inversion determines the most likely state vector to explain the observation ${\mathit{L}}_o$. It includes background knowledge with a multivariate Gaussian prior over state vector elements, with mean ${\mathit{\mu }}_a$ and covariance ${\mathbf{\Sigma }}_a$. We fit this distribution as in , with a library of diverse reflectance spectra. Shrinkage regularization ensures that the inversion can represent spectra not spanned by original library subspace. The inversion balances this prior against the measurement noise in Eq. (), which is determined from a component-wise instrument performance model . The optimal state vector minimizes the following cost function:

9$$\begin{array}{rl}\mathit{\chi }\left(\mathit{x}\right)& =\left(\mathit{F}\mathbf{\right(}\mathit{x}\mathbf{)}-{\mathit{L}}_o{)}^T{\mathbf{\Sigma }}_e^{-\mathrm{1}}(\mathit{F}\mathbf{\left(}\mathit{x}\mathbf{\right)}-{\mathit{L}}_o)\\ & +(\mathit{x}-{\mathit{\mu }}_a{)}^T{\mathbf{\Sigma }}_a^{-\mathrm{1}}(\mathit{x}-{\mathit{\mu }}_a),\end{array}$$ where ${\mathbf{\Sigma }}_e$ is the instrument noise. This cost is proportional to the negative logarithm of the posterior probability, a product of multivariate Gaussian prior and likelihood terms. We solve it with a trust region gradient descent optimization, a nonlinear optimization approach that respects positivity constraints on the free parameters . In principle, any gradient-based optimization would suffice with an appropriate starting point. The inversion typically converges in 10–20 iterations.

2.3 Postprocessing and structure functions

Before analyzing the spatial structure of the resulting H${}_{\mathrm{2}}$O maps, we perform several postprocessing steps to improve the map consistency. Even a highly accurate retrieval is likely to suffer some biases due to the influence of the surface type or magnitude. For example, very slight errors in atmospheric path radiance estimates are proportionally larger for dark targets, which can induce a spurious dependence between albedo and water vapor. Surface reflectance features that overlap water vapor absorptions can also influence retrievals. Finally, minor differences in the radiometric response or linearity of different cross-track elements can create striping artifacts in derived products. We address all of these issues with a single empirical correction. We first define a grid of locations $i\in \mathcal{L}$ where the retrieval is performed. We assert that the retrieved water vapor at each location $i$, written ${\widehat{q}}_i$, is the combination of an underlying “true” water vapor signal $q_i$, zero-mean measurement noise ${\mathit{ϵ}}_q$, and non-stochastic interference by surface and systematic instrument effects. We define a feature vector ${\mathit{\psi }}_i$ to include the surface reflectance in each channel and a sparse position vector $\mathit{p}$ encoding the associated cross-track position on the focal plane array (FPA) in a stacked binary representation. The position encoding enables the model to represent radiometric sensitivity of different FPA elements. The column vector ${\mathit{\psi }}_i$ thus combines ${\mathit{\rho }}_i$ and ${\mathit{p}}_i$. For simplicity, we assume the surface interference is a linear combination of these features, weighted by coefficients $\mathit{\varphi }$:

10$${\widehat{q}}_i=q_i+{\mathit{\varphi }}^T{\mathit{\psi }}_i+{\mathit{ϵ}}_q.$$ To estimate $\mathit{\varphi }$, we treat the true water vapor component as a random variable. Over large spatial scales, biases due to diverse content will tend to average out, and the sample mean of retrieved values can be used to approximate the true mean. Rearranging algebraically, we define the local vapor anomaly ${\mathcal{A}}_i$ as 11$${\mathcal{A}}_i={\widehat{q}}_i-\left[{\displaystyle \frac{\mathrm{1}}{\left|\mathcal{L}\right|}}\sum _{j\in \mathcal{L}}{\widehat{q}}_j\right]={\mathit{\varphi }}^T{\mathit{\psi }}_i+{\mathit{ϵ}}_{\mathrm{w}},$$ where ${\mathit{ϵ}}_{\mathrm{w}}$ now incorporates variability due to retrieval noise as well as variability in the true water vapor field. Since ${\mathit{ϵ}}_w$ is zero mean, Eq. () reduces to a straightforward linear regression problem. For a data matrix $\mathit{U}$ with rows made up of all reflectance spectra, and a vector $\mathcal{A}$ of anomaly values, we estimate $\mathit{\varphi }$ with ordinary least-squares regression. We subtract the predicted error from the original estimate to yield the bias-corrected estimate $q_i^{\prime }$: 12$$q_i^{\prime }={\widehat{q}}_i-{\widehat{\mathit{\varphi }}}^T{\mathit{\psi }}_i\text{for}\widehat{\mathit{\varphi }}=\mathit{U}({\mathit{U}}^T\mathit{U}{)}^{-\mathrm{1}}{\mathit{U}}^T\mathcal{A}.$$ To apply the bias correction, we segment long flight lines into segments of no longer than 2000 pixels each, and we apply the interference correction to each segment independently.

A final postprocessing operation smooths each water vapor image with a Gaussian spatial filter. This dramatically reduces the retrieval noise, making it possible to resolve much finer structures. This noise reduction is also beneficial for estimating structure functions, which are bounded artificially on the low end by the average squared noise in neighboring pixels. We use leave-one-out cross-validation to select an optimal blurring kernel width. To score a candidate width, we compare every point in the scene to the prediction made by applying the blurring kernel to that location, excluding the point under test with appropriate renormalization. The optimal kernel standard deviations range from four to six in different scenes; we use the lower value to preserve fine spatial structure at sub-100 m scales.

After these postprocessing steps, we calculate empirical structure functions for each pixel shift distance $d$, computing the squared differences between image locations shifted in the along-track direction. Recall that the long flight lines are corrected in parts no longer than 2000 pixels in length, revealing structure functions up to 4 km. To avoid minor offsets in the vapor field on each side of a border, we do not permit shifts across multiple segments. We aggregate the statistics of all segments to form a structure function estimate for each flight line. We manually mask any visible artifacts in water vapor images to exclude them from the calculation. These artifacts are primarily due to the presence of clouds or their shadows. Clouds can also disrupt their local light fields with scattered illumination, an effect visible in water vapor maps as obvious halos around clouds and their shadows. To mitigate this effect, we dilate the masks by 200 m horizontally.

3 Results: AERONET overflights

We first assess the absolute accuracy of water vapor absorption measurements using airborne overflights of the AERONET robotic observation network . The AERONET sun photometers view the sun directly, estimating water vapor by solar extinction . Consequently, they do not measure the same optical path as the down-looking sensor. Nevertheless, these coincident overflights validate the approach and provide a ceiling for unmodeled uncertainties. The airborne instrument is AVIRIS-C . It flew on an ER-2 aircraft at approximately 20 km altitude. Our dataset spans 6 years of operations over California, from 2013–2019, during which it overflew active terrestrial AERONET sites on over 100 occasions. Our atmospheric model for these observations used an Air Force Geophysics Laboratory (AFGL) midlatitude summer profile . Aerosols were light throughout the overflights, so we omitted them from the state vector.

We use several filtering methods to remove uncertain data points. Many overflights are contaminated by clouds. Our comparison excludes obvious clouds recognized by an estimated surface reflectance at 450 nm above 0.2. However, even when opaque clouds do not directly cover the AERONET location, high-altitude cirrus, cloud shadows, or scattered irradiance can contaminate a retrieval. Additionally, high spatiotemporal variability in water vapor can cause a discrepancy; it would aggravate the difference in optical paths, as well as the imperfect temporal coincidence between the AERONET acquisitions and the overflight. Considering the strong influence of ${\mathit{\theta }}_{\mathrm{s}}$ demonstrated in Fig. , and that ${\mathit{\theta }}_{\mathrm{s}}$ during typical flights often reaches 45${}^{\circ }$ or greater, the ground-based measurement could see a very different vapor field if humidity were not horizontally homogeneous. To address these issues, we remove any data points where the in situ standard deviation in H${}_{\mathrm{2}}$O is larger than 0.1 g cm${}^{-\mathrm{2}}$, estimated by comparing the 10 temporally closest acquisitions. After this filtering step, 64 data points remain. Some discrepancies in the optical paths remain, which become larger for column water vapor in the free troposphere than in the planetary boundary layer. As neither measurement resolves the detailed vertical profile, it is not possible to distinguish these cases without additional information.

Figure  compares the AERONET and remote retrievals for all sites. Error bars indicate the variability in water vapor over the 10 closest time steps (AERONET) or the 10 $\times$ 10 enclosing rectangle of pixels (AVIRIS-C). Each AERONET time step is approximately 5 min, so our window provides a measure of variability in the hour around the flight. These overflights represent a wide range of atmospheric conditions and solar angles. They span column water vapor concentrations from approximately 0.3 to 2.5 g cm${}^{-\mathrm{2}}$. The airborne and up-looking measurements show strong agreement, with a correlation coefficient $R$ of 0.95. A few AVIRIS-C observations with very high variability are likely clouds or cloud shadows that survived the filtering process. Nevertheless, the result is broadly consistent with AERONET accuracy of 12 %–15 % for column water vapor, a claim which has been validated independently by comparison to microwave radiometers . This provides confidence in the accuracy of the water vapor retrieval. There is a small bias between the two datasets. However, even if we attributed it entirely to the airborne data, such an offset would not influence the structure function measurement.

Figure 5

Coincident AVIRIS-C vs. AERONET water vapor measurements, filtering out AERONET instances with high temporal variability ($\mathrm{1}\mathit{\sigma }>\mathrm{0.1}$ g cm${}^{-\mathrm{2}}$ over 10 temporally closest measurements).

[Figure omitted. See PDF]

4 Results: structure function measurement

We apply the H${}_{\mathrm{2}}$O retrieval methodology to four flight lines from a 2018 AVIRIS-NG campaign in India. We selected these flight lines for favorable solar-observing geometry and because they represented two distinct conditions observed on different days. All flight lines were at 4 km altitude over water, providing a uniform, topographically flat surface for vapor retrievals. On 12 May, two flight lines encountered scattered low clouds. On 14 May, two more flight lines encountered clear skies. We restricted our study to flights with solar zenith angles less than 10${}^{\circ }$, a condition which occurred on 2 flight days. This provided effective spatial resolutions from 80 to 250 m at nadir. In comparison, the native spatial sampling of AVIRIS-NG at these altitudes was 4 m. In addition to improving spatial resolution, the extreme solar angles also produced significant sunglint, revealing near-infrared water features that would otherwise be masked by the absorption of the water surface. Information on each flight line appears in Table . We applied the standard AVIRIS-NG radiance calibration procedure , including corrections to the spectral response function . We then performed water vapor retrievals on all flight lines using a tropical atmospheric profile . AVIRIS-NG had a higher intrinsic signal-to-noise ratio (SNR) than AVIRIS Classic, and the solar angles provided more signal than a typical AERONET overflight from the validation experiment. Consequently, we expected a more sensitive retrieval in the India datasets.

Table 1

Flight lines used for structure function analysis.

Label	Flight line ID	Date	${\mathit{\theta }}_{\mathbf{s}}$	Nadir ${\mathit{R}}_{\mathbf{e}}$ (m)	Latitude (${}^{\circ }$ N)	Longitude (${}^{\circ }$ W)	Length (m)
A	ang20180512t052609	12 May 2018	9.7	250	21.649	87.775	24000
B	ang20180512t053942	12 May 2018	6.9	177	21.639	87.789	24000
C	ang20180514t055115	14 May 2018	4.1	105	21.524	88.346	8000
D	ang20180514t060206	14 May 2018	3.1	79	21.526	88.325	8000

To provide context for interpreting the vapor fields, we analyzed MERRA-2 reanalysis data for each day. The atmospheric conditions were generally similar across the overflights, with light trade winds and the lack of an obvious inversion to stratify boundary layer processes. There were also some differences between the days. Wind velocity changed slightly but was not obviously tied to any change in atmospheric turbulence. The relative humidity was generally higher on 14 May. The lapse rate was slightly more variable: 8 K km${}^{-\mathrm{1}}$ at 760 hPa as opposed to 7 K km${}^{-\mathrm{1}}$ on 12 May. These changes would be consistent with a slightly more turbulent atmosphere and a shallower scaling exponent, though there was no obvious step change in atmospheric stability. Figure  shows the temperature, lapse rate, and specific and relative humidity for each day.

Figure 6

MERRA-2 reanalysis predictions for the 2 flight days. Counterclockwise from upper left: temperature, specific humidity, lapse rate, and relative humidity.

[Figure omitted. See PDF]

Figure  shows a typical spectrum fit. The top panel shows the radiance measurement. The converged model matches the measurement closely, with the largest discrepancies in the shortest wavelengths. The middle panel shows the residual error. Some structure at the sub-2 % level is consistent with minor model discrepancies or calibration errors. These mostly affect the shortest wavelengths outside the dominant water absorption intervals. The reflectance spectrum in the bottom panel shows a good-quality retrieval of a water surface. There is considerable sunglint, a spectrally flat additive signal that uniformly elevates the spectrum. This raises the near- to shortwave-infrared reflectance, which would normally be zero over liquid water, to values of 12 %–13 %. Elevated reflectance in the presence of sunglint enables high-accuracy retrievals of $q_{\mathrm{v}}$ over open water. This is a reliable consequence of small solar zenith angles under a wide range of wind conditions.

Figure 7

(a) Example radiance spectrum. (b) Model fit residual. (c) Estimated reflectance with 95 % confidence bounds.

[Figure omitted. See PDF]

Figure  shows the column water vapor maps for individual segments of the flight lines. In the top panels, dark blue areas indicate manual masks applied to exclude cloudy regions. After the interference correction, the vapor maps are generally free of artifacts except for a minor band at the zero phase angle. This artifact, indicated by a white arrow, may be related to issues with the sun glint bidirectional reflectance distribution or the aircraft shadow. Therefore, it was excluded from the statistics. The smallest-scale variability is almost certainly related to noise in the retrieval algorithm. Larger-scale coherent structures show that meaningful differences in water vapor content are resolved at the level of 1 %–3 %, over scales comparable to the predicted spatial resolutions of 80–250 m. The center and bottom panels show example maps for flight lines C and D on 14 May, which were entirely free of clouds. In aggregate, the flight lines encountered water vapor values ranging from about 3.23 to 3.42 g cm${}^{-\mathrm{2}}$. The smoothing operation dramatically increases the contrast of spatial features (Fig. ).

Figure 8

H${}_{\mathrm{2}}$O vapor maps. In flight lines A and B, dark blue areas are masked to avoid clouds. A white arrow in flight line B indicates an artifact that occurs where the solar phase angle is zero. A white box in flight line C indicates the area shown by Fig. .

[Figure omitted. See PDF]

Figure 9

(a) Initial water vapor map from flight line C. (b) Noise reduction after kernel smoothing reveals fine-scale structure at sub-kilometer scales.

[Figure omitted. See PDF]

Figure  shows structure functions for the four flight lines with second-order scaling exponents for representative intervals. Grey lines indicate the best-fitting power-law exponents over the 500–1000 m interval. We also plot the canonical $\mathrm{2}/\mathrm{3}$ slope for reference. The flight lines' configurations of vapor and clouds were all unique, but their scaling functions from the common flight days show similar profiles. This intra-day similarity demonstrates the repeatability of the measurement in a consistent air mass. For flight lines A and B, ${\mathit{\zeta }}_{\mathrm{2}}$ is $\mathrm{0.64}$ and $\mathrm{0.63}$, respectively. These values are very close to the classic $\mathrm{2}/\mathrm{3}$ value of a passive scalar in turbulence. The slope continues to scales well under 1 km, with a steepening of the curve below about 500 m. This corroborates the VOCALS-REx data analyzed by , although those observations were taken in a stable regime within and above stratocumulus clouds off of the coast of South America. In this case, the observed steepening continues to scales approaching the effective spatial resolution of approximately 250 m. The right panel shows flights C and D on 14 May. These flight lines have shallower profiles in the 500–1000 m range, with scaling exponents of $\mathrm{0.41}$ and $\mathrm{0.29}$. The curves are less consistent across the two overflights, which might be due to small sample sizes. Flight lines C and D contained just one segment of 2000 m, compared with six segments in each of flight lines A and B.

Figure 10

Structure functions for the maps in Fig. . (a) Flight lines A and B. (b) Flight lines C and D.

[Figure omitted. See PDF]

To confirm that the results were robust to slight differences in viewing angles, we recalculated the same result using the central half of the flight line data. This excluded all off-nadir view angles greater than approximately 9${}^{\circ }$. The results did not change significantly; flight lines A, B, C, and D showed ${\mathit{\zeta }}_{\mathrm{2}}$= 0.55, 0.58, 0.41, and 0.29, respectively. These differences would not be large enough to change our interpretation. In all cases, the maps comprise millions of independent measurements, making it possible to resolve structures of magnitudes smaller than the single-pixel noise. All acquisitions, including the extreme shallow scaling of flight lines C and D, are in the range of exponents predicted for convective air masses by . This is consistent with the presence of convective clouds in the flight lines.

5 Discussion

The airborne experiments demonstrate the ability of VSWIR spectroscopy to measure total column water vapor with high precision but cannot directly apportion this observed variability to different vertical layers. The state of the lower troposphere generally dominates the overall water vapor content, but it does not necessarily follow that the troposphere determines the observed variability over short spatial scales. Our derived ${\mathit{\zeta }}_{\mathrm{2}}$ coefficients refer to the horizontal variability of water vapor integrated along the full solar-reflected optical path, and variability at different altitudes could contribute to estimated ${\mathit{\zeta }}_{\mathrm{2}}$. For example, one might observe a superposition of different atmospheric regimes at different layers, such as convective and non-convective, or boundary layer and free troposphere, which would blur the variability in the total column. This could be significant for attributing the variability or for comparing the observations with more localized measurements. Given that this study is one of the first that explores water vapor variability at such high horizontal resolution, there is a dearth of independent data against which to interpret vertical sensitivity.

We can gain insight by framing the relationship between full column and stratified variability as a statistical question. Other vertically resolved measurements give evidence that the derived exponents reflect structure in the lower troposphere and planetary boundary layer. Preliminary support comes from reprocessed data from the five 2019 Cal/Val campaign flights of the High Altitude Lidar Observatory (HALO) over the eastern Pacific . HALO’s differential absorption lidar (DIAL) water vapor retrievals were reprocessed for our analysis with 3 km along-track averaging, 0.5 km vertical averaging for altitudes 4 km or lower, and 1 km vertical averaging for altitudes higher than 4 km. This along-track averaging provided an optimal balance of spatial resolution and sensitivity. We integrated the column water vapor vertically to obtain partial column water vapor (PCWV) and calculated the squared Pearson correlation coefficient ($r^{\mathrm{2}}$) between the PCWV up to different altitudes, and that up to 7 km, which approximated the total column water vapor (TCWV). Picking a higher altitude would have severely restricted the number of valid lidar footprints due to the flight altitudes, and there was also very little vapor content above 7 km. To create a dataset which was most similar to the AVIRIS-NG flight data, we extracted HALO data over ocean in the least-cloudy parts of each flight, with valid water vapor retrievals from the lowest bin up to 7 km. Measurements occurred approximately every 200 m and were averaged over 3 km. We required at least 90 % of the 200 m footprints to be non-cloudy within a given 3 km along-track measurement, and then we identified the “least cloudy” larger areas by moving a 750 km window along the flight line and requiring $>\mathrm{20}$ % of the 3 km spacing profiles be clear. The resulting partial column water vapor quantities appear in Fig. .

Figure 11

Partial column water vapor for five HALO flights calculated at different capping altitudes.

[Figure omitted. See PDF]

Figure  shows that in every flight, over 70 % of the variance in TCWV is explained by altitudes less than 2 km, and more than 90 % of the variance is explained by variance at altitudes less than 3 km. While the horizontal resolution and meteorology of HALO and AVIRIS-NG flights are different, this provides quantitative evidence that the ${\mathit{\zeta }}_{\mathrm{2}}$ values derived here most likely refer to bulk variability within the PBL. This adds to mesoscale studies finding that processes in the PBL largely control the total atmospheric water vapor content . It also suggests that, as a component of a larger measurement and assimilation system, the total column measurement could provide probabilistic constraints on PBL water vapor.

Figure 12

Data from HALO flights suggest that the PBL dominates the total column variability. Correlation coefficients between the lowest atmospheric layer and the total column water vapor for different capping altitudes.

[Figure omitted. See PDF]

6 Conclusion

This paper describes an approach for mapping column-average atmospheric water vapor at sub-kilometer spatial scales with remote VSWIR imaging spectroscopy. We validate the method by comparison with in situ AERONET observations. We then map column water vapor using reflected solar sunglint over ocean surfaces in several flight lines from an airborne campaign where favorable solar angles permitted a uniquely high-resolution measurement. We find scaling behavior broadly consistent with prior studies, with high repeatability across different observations of the same air mass. The experiment is limited in some respects; for example, it provides only a total-column measurement, making it less compatible with height-resolved structure functions . Moreover, stringent requirements on observation geometry limited the dataset to a few representative flight lines. Nevertheless, it provides a proof of concept that VSWIR imaging spectroscopy, primarily used for measurements of surface phenomena, can also provide accurate water vapor maps for atmospheric studies. The key innovation afforded by VSWIR water vapor retrievals will be the spatial scale, accuracy, and retrieval capability over land surfaces – all of which improve on various aspects of MODIS bi-spectral or passive microwave imaging techniques.

Future measurement campaigns aiming to observe water vapor with this approach should consider both the solar geometry and signal-to-noise level, either of which could limit the resolution. In this study, we limited the influence of solar geometry by careful selection of observing conditions. When this is not possible, computer techniques might be used to remove the spatial blurring effect of low solar angles. The effective spatial response function could be calculated and deconvolved from the measured vapor field. It may also be possible to reduce the impact of solar zenith angles by measuring the structure function orthogonally to the solar direction, a premise that will be explored in future work. Regardless of whether such compensation is possible, the next generation of orbital VSWIR imaging spectrometers should significantly increase the data available for fine-scale mapping of atmospheric water vapor.

Code and data availability

Code used for the calculations described in this paper is available in the ISOFIT code repository under the Apache 2.0 open-source license (https://github.com/isofit/isofit, last access: 21 March 2021; DOI: 10.5281/zenodo.4614338, ). The airborne datasets described in this paper are available from the Jet Propulsion Laboratory (https://avirisng.jpl.nasa.gov, last access: 21 March 2021, , for AVIRIS-C). AERONET data are available from the AERONET web site (https://aeronet.gsfc.nasa.gov/data_push/V3/AOD/AOD_Level20_All_Points_V3.tar.gz, last access: 23 March 2021, , 10.1016/S0034-4257(98)00031-5).

Author contributions

DRT contributed to the writing, algorithm design, experimental concept, and execution. PGB contributed to code development, algorithm design and implementation, interpretation of the experiments, and writing. ROG planned, coordinated, and managed the AVIRIS-C and AVIRIS-NG campaigns. BHK provided scientific motivation, analysis of MERRA-2 data, and writing. MDL contributed to the measurement concept, data interpretation, and writing. MR contributed to the measurement concept, code evaluation and validation, analysis of HALO data, and writing.

Competing interests

The authors have no competing interests with the investigation described in this document.

Special issue statement

This article is part of the special issue “Analysis of atmospheric water vapour observations and their uncertainties for climate applications (ACP/AMT/ESSD/HESS inter-journal SI)”. It is not associated with a conference.

Acknowledgements

We thank the AERONET station principal investigator(s) and co-investigator(s) and their staff for establishing and maintaining the sites used in this investigation. Some of the Cimel sun photometer data were collected by the U.S. Department of Energy as part of the Atmospheric Radiation Measurement Program User Facility (ARM) and processed by the National Aeronautics and Space Administration's Aerosol Robotic Network (AERONET). We thank the members of the AVIRIS-NG team who participated in data acquisition and analysis, including Michael Eastwood, Sven Geier, Mark Helmlinger, Winston Olson-Duvall, and Sarah Lundeen. AVIRIS-NG is sponsored by the National Aeronautics and Space Administration (NASA) Earth Science Division. We also thank Amin R. Nehrir and Brian Carroll, who assisted with HALO data analysis. This research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. We also acknowledge the NASA Earth Science Division's AVIRIS-NG instrument and the data analysis program “Utilization of Airborne Visible/Infrared Imaging Spectrometer Next Generation Data from an Airborne Campaign in India”. Any use of trade, firm, or product names in this publication is for descriptive purposes only and does not imply endorsement by the US government. Copyright 2021 California Institute of Technology. Government sponsorship acknowledged.

Financial support

This research has been supported by the National Aeronautics and Space Administration (NASA) (grant no. AO NNH16ZDA001N-AVRSNG).

Review statement

This paper was edited by Martina Krämer and reviewed by Christoph Kiemle and David Adams.