Modified

24 July 2026

Level-1.1: Self-heating correction

TODO IN PROGRESS

We refined self-heating correction for open-path IRGAs described in earlier studies (Burba et al., 2006; Järvi et al., 2009, Kittler et al., 2017) and implemented the code in the Python library diive.

The refined correction is a specialized procedure in eddy covariance flux processing designed to correct a systematic measurement error in open-path \(\text{CO}_2\) analyzers, primarily due to solar heating of the instrument body. The entire workflow is modularized into three main stages: Physics (calculation of the unscaled flux correction term \(\text{FCT}_{\text{unsc}}\) , optimization (deriving scaling factors (\(\text{SF}\) for \(\text{FCT}_{\text{unsc}}\) ), and application (applying the final correction).

Physical basis

The core problem addressed is the dilution effect caused by the OP-IRGA sensor head heating up in the sun. This heating creates a thermal plume of warmer, less dense air around the sensor sampling path. Since the OP-IRGA measures molar density (\(q_c\), \(\text{mol m}^{-3}\)), this low-density plume causes the measured \(\text{CO}_2\) concentration to appear artificially low, often leading to a persistent negative flux bias (apparent \(\text{CO}_2\) uptake).

1. Calculation of the Unscaled Flux Correction Term (\(\text{FCT}_{\text{unsc}}\))

The ScopPhysics class calculates the initial correction term based on boundary-layer dynamics. It supports three methods, notably those from Järvi et al. (2009) and Burba et al. (2006, 2008).

The unscaled term (used by JAR09/BUR06 methods) is calculated as:

\[\text{FCT}_{\text{unsc}} = \frac{(T_s - T_a) q_c}{r_a T_a K} \left(1 + 1.6077 \frac{\rho_v}{\rho_d}\right)\]

Where:

  • \(T_s\) and \(T_a\): Instrument surface and ambient air temperature (\(\text{K}\) or \(\Delta T\)).
  • \(q_c\): \(\text{CO}_2\) molar density (\(\mu \text{mol m}^{-3}\)).
  • \(r_a\): Aerodynamic resistance to heat transfer (\(\text{s m}^{-1}\)), calculated from wind speed (\(u\)) and friction velocity (\(u*\)).
  • \(T_a K\): Ambient air temperature in Kelvin.
  • \(\left(1 + 1.6077 \frac{\rho_v}{\rho_d}\right)\): The water vapor correction factor, analogous to the WPL correction’s effect on temperature/density.

The _estimate_surface_temp_bur06 and _estimate_surface_temp_jar09 methods provide parameterized relationships to calculate \(T_s\) from \(T_a\) for this equation.

2. Optimization: Deriving the Scaling Factor (\(\text{SF}\))

The calculated \(\text{FCT}_{\text{unsc}}\) is derived from theoretical/empirical models, but it requires tuning for a specific instrument and site. The ScopOptimizer class performs this step, using parallel measurements from an Enclosed-Path (CP) IRGA, which is not subject to the self-heating effect, as the flux reference.

Cost Function and Bootstrapping

The goal is to find the Scaling Factor (\(\text{SF}\) or \(\xi\)) that minimizes the cumulative difference between the corrected OP flux and the CP reference flux. The optimization uses a custom vectorized Cost Function (_cost_function_numpy) with the objective:

\[\text{Minimize} \sum \left| \sum \left( \text{Flux}_{\text{OP, raw}} + \text{FCT}_{\text{unsc}} \cdot \xi \right) - \sum \left( \text{Flux}_{\text{CP, ref}} \right) \right|\]

This minimizes the \(\text{L}1\) norm of the difference in cumulative sums, which ensures the long-term, integrated \(\text{CO}_2\) budget matches between the two instruments. The optimization is performed using scipy.optimize.minimize_scalar within defined bounds. To ensure robustness, the process uses Circular Block Bootstrapping (_block_bootstrap_indices) to resample time-series data while preserving its autocorrelation, yielding a median \(\text{SF}\) and associated confidence intervals for each bin.

Binned Optimization

The optimization is performed independently for bins of a class variable (e.g., \(u*\)) and separated by daytime/nighttime conditions, acknowledging that instrument self-heating mechanisms vary significantly with solar radiation and turbulence.

3. Application: Final Correction

The ScopApplicator class takes the optimized scaling factors (scaling_factors_df) and applies the final correction to the entire time series.

The _assign_scaling_factors method maps the derived median \(\text{SF}\) from the optimization table back to the original time series based on the value of the class variable and the daytime flag, using an efficient pd.merge_asof lookup. Finally, the total correction term is calculated, and the corrected OP flux is computed:

\[\text{FCT} = \text{FCT}_{\text{unsc\_gf}} \cdot \text{SF}\]

\[\text{Flux}_{\text{OP, corr}} = \text{Flux}_{\text{OP, raw}} + \text{FCT}\]

The use of gap-filled (\(\text{gf}\)) \(\text{FCT}_{\text{unsc}}\) (provided by RandomForestTS in the ScopPhysics._gapfill method) ensures the correction can be applied even when auxiliary variables required for the physical model were missing.

Back to top