the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 31 Jul 2026
A continental benchmark dataset for evaluating ecosystem gradient-flux approaches across 47 NEON flux towers
Abstract. Benchmark datasets for evaluating profile-based ecosystem flux methods across environmental gradients are currently lacking. This limits method evaluation, cross-site intercomparison, and the expansion of tower-based monitoring for gases not routinely measured by the eddy covariance (EC) method. Here, we present a benchmark dataset derived from 47 terrestrial towers in the National Ecological Observatory Network (NEON), integrating co-located EC, concentration profiles, tower geometry, and canopy structural metrics. We developed a dataset to evaluate the performance of three widely used gradient flux approaches – the modified Bowen ratio (MBR), aerodynamic (AE), and wind-profile (WP) methods – against co-located EC measurements of CO₂ and H₂O. We evaluate how canopy structure, sensor height configuration, and data filtering influence agreement with EC across ecosystems, with canopy heights ranging from 0.15 to 53 m. Performance varied strongly among approaches and measurement configurations. Across all ecosystems, 11% of height-pair combinations for CO₂ and 19% for H₂O achieved moderate-to-strong concordance with EC, with the MBR approach providing the most consistent performance. Reliable estimates most often occurred when both sampling heights were above the canopy (AA) or when one level was above the canopy and the other within sparse canopies (AW). Concordance declined as canopy complexity increased, highlighting persistent challenges in tall, complex forests. Since filtering substantially reduced data availability, we developed an ensemble framework that combined reliable estimates across the three GF approaches and height pairs. Ensemble GF fluxes reproduced seasonal diel dynamics measured by EC across ecosystems (R2 = 0.58–0.99), capturing both the seasonal magnitude and direction of net ecosystem exchange. Benchmark analyses show that the GF method can provide robust ecosystem-scale flux information when deployed under favorable structural and micrometeorological conditions. The released benchmark dataset provides a resource for testing new flux algorithms, optimizing sensor placement, benchmarking profile methods, and extending multi-gas tower observations to be more inclusive of gases that are difficult to measure with fast-response EC instrumentation, including CH₄, N₂O, volatile organic compounds, and reactive pollutants.
- Preprint
(6322 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
- RC1: 'Comment on essd-2026-413', Anonymous Referee #1, 09 Sep 2026
-
RC2: 'Comment on essd-2026-413', Anonymous Referee #2, 03 Oct 2026
The study presented a profile-based CO2 and H2O flux dataset derived from NEON datasets, along with collocated eddy-covariance fluxes, concentration profiles, tower geometry, and canopy structural metrics. It tested and summarized the performance of three gradient flux approaches across 47 NEON sites and proposed an ensemble data product for further research.
I think the research initiative, generated datasets, and the paper itself are important in many ways, both for research and application. In particular, I appreciate the authors’ effort and comprehension in designing, implementing, and discussing the study and its outcome. I would recommend the paper for consideration for publication in Earth System Science Data after careful revision. Here, I provide some general and specific comments to help improve the presentation and, most importantly, the usability of the dataset in future research. Since this is a data-focused journal, I don’t insist that all comments be addressed directly with additional analyses. However, I suggest the authors consider adding these points to the discussion and/or providing the necessary data and metadata to enable further research.
[1] As the authors pointed out, the goal of this paper is to test when and where these gradient-based flux estimates work, and why. I think it is equally important to present the estimated flux dataset and also metadata that explain the validity or robustness of the methods and dataset. In particular, I suggest including information or data on the canopy’s vertical structure. The extracted canopy metrics are valuable for exploratory analyses across sites, but they are limited if one intends to dive further into the details of each case and generalize best practices for future deployment. Instead, the vertical profiles of structures (e.g., LAI, biomass) would be much more informative. Figure 9 is great for its abstraction but still contains a lot of information. But looking at the results between short- and tall-canopy systems, I wonder whether it makes more sense to use the vertical structure profile to separate the canopy into crown and trunk space (esp. forests), potentially revising the groups in Table 1.
[2] Along these lines, I also suggest including information or data on the vertical distribution of sources/sinks for both CO2 and H2O. As a best guess, I assume mid-canopy (peak foliage) and soil to be two major source/sink levels. Diving deeper, this assumption depends on the target (CO2 vs H2O), time (day vs night, summer vs winter), and ecosystem (dense forest vs savanna/shrub). I think it would be valuable to generate and present monthly or seasonal average CO2/H2O concentration profiles and use them to inform the source/sink distribution. This could help explain performance differences among sites and between CO2 and H2O.
[2] Section 2.1: Consider adding a few more details to avoid confusion or for readers who are less familiar with NEON or gradient flux methodology: 1) I assume at least one of the sampling levels needs to be above the canopy for all pairs, since all categories begin with “A”. But I got confused when reading Figure 6C, which suggests cases where the top level is below canopy height; 2) Please define (quantitatively) the canopy heights and the relative location above/within/below the canopy heights. The categories in Table 1 are good for summary and grouping, but please provide the quantitative metrics. If possible, consider revising Figure 1B to show all pair categories. 3) How does the study address ecosystems that experience changing canopy, either in height or foliage? When were the canopy metrics measured, and how often were they resampled?
[3] Section 2.2: Consider adding or clarifying a few more details: 1) Was each CO2/H2O sampling level sampled once (for 6 min) over each 30-min period? 2) The description in Lines 169-170 is a bit unclear, particularly for data users who have no experience with this type of instrument; 3) Lines 174-175: It’s unclear where and which variable this alignment occurred; 4) A side note: The NEON profile was designed to have the highest level at the eddy covariance heights, but not above, right? Would this limit or somehow bias this research’s conclusion or implication?
[4] Section 2.4.1: Please specify the known flux and known gas gradient.
[5] Section 2.4.2: Please specify the data sources used in equation 6.
[6] Section 2.4.3: Please specify the derivation of d. Also, I assume it should be u measured at height z. Please be specific.
[7] Section 2.7: Most of these metrics focus on the structures, and are more useful in explaining the site-to-site difference. As noted above, I think some information on vertical structure and source/sink profiles may help this and future analyses. For example, what is the relative height of the 2nd sampling level compared with the peak LAI height? Or, what is the relative height of the 2nd sampling level compared to the maximal source/sink in the profile? Following this line of discussion, the potential seasonal changes in structure (e.g., canopy height, foliage) and in source/sink between day and night (e.g., photosynthesis vs respiration) may need to be considered.
[8] Figure 8 seems to be missing in the review version.
[9] Discussion paragraph (Line 488-497). This is a great discussion that provides context and background for this research and dataset. I suggest moving some of these points up front to help set expectations.
[10] Section 4.2: 1) Please double-check and consolidate the definition of Ha and Hb here and in Figure 6. 2) This is one of the places where I think additional information on the vertical profiles (structure or source/sink) would help the discussion about the distinct patterns between short and tall canopies. 3) I suggest the author discuss the implications of the findings at short vegetation sites, particularly for those planning to apply the method somewhere else. The 5-10 m measurement height above short vegetation may not always be logistically feasible and may suffer from limited fetch.
[11] Section 4.3: Reading the comparison between EC and gradient flux, I wonder whether the authors could also comment on whether it is feasible or practically useful to use flux gradient (e.g., slow-response sensors) as a complementary measurement at EC flux towers, e.g., for benchmarking, QA/QC, or gap-filling.
Citation: https://doi.org/10.5194/essd-2026-413-RC2 -
AC1: 'Comment on essd-2026-413 AC1', Sparkle Malone, 06 Nov 2026
Authors Response to RC1
We thank Referee 1 for the thoughtful assessment and detailed comments, which improved the paper's clarity.
Overview of the planned revision
Below, each referee comment is quoted in full and italicized, followed by our response and the changes we will make to the manuscript, shown in shading and quoted from our revised draft. Section, figure, table, and equation numbers refer to the revised manuscript. Line numbers will be provided in the tracked-changes manuscript submitted after the discussion.
While addressing both referees' comments, we extended the record, corrected and documented the processing, and reran every calculation. The revised manuscript will include five changes that affect most of the numbers in the paper:
Record extended from August 2021–July 2024 to August 2021–August 2026
Two more years at every site
Reliable sampling height pairs (RSHPs) at 42 sites for CO₂ and 43 for H₂O; 42 sites have them for both gases, compared with 30 in the discussion paper
Spike removal: values outside the median ± 6 scaled MAD per site, approach, gas, pair, day or night, and season
Isolated extreme values from near-zero tracer gradients or large eddy diffusivities (Referee 1, Fig. A13)
2.4% of about 5 million filtered values removed; median site maximum of |MBR CO₂| from 473 to 104 µmol m−2 s−1 (Sect. 2.5; new Appendix B)
Processing corrections (listed below)
Text and code disagreed in places, found while documenting Eqs. 2–7 for Referee 2 and in an independent code review by co-authors
All fluxes and evaluations recalculated; main conclusions hold, and the secondary conclusions that changed are listed below
Ensemble described as implemented (inverse-RMSE weighted) and released with the spread between pairs
The text said the mean ensemble was used; Referee 1 noted the value of the ensemble for gappy data
New 30-min ensemble table in the data release with the spread, the number of pairs and the approaches for every half-hour (Sects. 2.8 and 5)
Mean canopy height hc taken from the NEON airborne lidar canopy height model throughout
The discussion paper mixed lidar and NEON site-metadata heights
Site canopy heights span 0.03–38 m (previously 0.15–53 m), consistently across the AA/AW categories, the random forest and the placement analysis (Sect. 2.1; Table 1)
Processing corrections. We will describe the methods as implemented after the following corrections:
- MBR (Eq. 2): the target flux is the median of 1000 bootstrap draws of F1ΔC2/ΔC1, because the mean of a ratio is unstable when ΔC1 is close to zero; the 2.5–97.5 percentile range of the draws is the uncertainty (Sect. 2.4.1).
- Aerodynamic canopy height h: obtained by inverting the logarithmic wind profile at the top inlet, using the height of the wind measurement, with d = 0.66h, z0 = 0.1h and a centered 20-day window (Sect. 2.4.2).
- Effective height: z is the geometric mean of each pair's two inlet heights; ζ = (z − d)/L is computed for each pair; k = 0.4 throughout (Eq. 6 used 0.41).
- Wind-profile method (Eq. 7): Kwp is set to missing when z − d ≤ z0, where the equation is singular (Sect. 2.4.3).
- Ensemble: weights use the direct root-mean-square difference between each pair's GF and EC, and each GF value is assigned to the half-hour containing the midpoint between its two inlet samples, the time to which EC is interpolated (Sect. 2.8).
- Diel curves: absolute diel values and daily totals in the discussion paper were half their correct values because of a unit conversion error; they are now in g m−2 h−1 (Sect. 2.9; Table A2). Correlations, slopes, relative biases and signs were not affected.
Key results, discussion paper to revised manuscript.
Sites with an RSHP: CO₂ / both gases
38 / 30
42 / 42 (43 for H₂O)
Reliable height-pair combinations: CO₂ / H₂O
11% / 19%
16% / 23%
RSHPs per site: CO₂ / H₂O
3.8 / 5.7
5.1 / 7.1
Mean CCC for CO₂, MBR / AE / WP
0.27 / 0.09 / 0.13
0.34 / 0.14 / 0.15
Share of CO₂ RSHPs that are AW
28%
38%
Random forest accuracy
95% (no-information rate 85%)
89.6% on held-out pairs (no-information rate 80.8%)
Ensemble used; mean site CCC across the six methods
mean; 0.65–0.66
inverse-RMSE weighted; 0.66–0.68
Sites with seasonal diel curves: CO₂ / H₂O
31 / 43
40 / 42
R² of GF vs EC diel curves
0.58–0.99
median 0.91 (CO₂), 0.99 (H₂O)
The main conclusions are unchanged: the MBR is the most reliable approach, concordance declines with canopy height and complexity, and the ensemble reproduces the seasonal diel cycles measured by EC. Three secondary conclusions changed and will be corrected. (1) The best AE and WP lower inlets in tall forests lie in the upper crown, not in the trunk space. (2) The statement that AW configurations raise AE and WP concordance from about 0.18 to 0.64 in tall forests did not hold and will be removed. (3) The random forest is evaluated on held-out pairs rather than in-sample.
The revision will also add a glossary (new Table 1, so later tables are renumbered by one), three appendix figures (A20 inlet heights relative to the canopy, A21 seasonal concentration profiles, and A22 seasonal aerodynamic canopy height), and three appendices: B on gradient-flux spikes, C on selecting pairs when no EC benchmark exists (relevant for CH₄), and D on seasonal day- and night-time CO₂ and H₂O concentration profiles at every site. Panel labels will be lowercase, following ESSD style. Our responses to Referee 2 are posted separately.
1.1 General comments
This manuscript uses 47 NEON sites to compare the eddy covariance and gradient flux methods for estimating CO₂ and H₂O fluxes. This analysis is valuable for assessing the conditions under which the two methods agree or disagree and is unique in that NEON infrastructure and data make a systematic comparison possible. A better understanding of how well gradient fluxes capture ecosystem gas exchange is particularly valuable for trace gas fluxes, which are more difficult or impossible to measure with eddy covariance. Overall, the analysis is extremely thorough and contains a wealth of information about how seasonality and ecosystem structure affect agreement between the two methods, as well as technically useful details such as measurement height for gradient fluxes. One conclusion that I found particularly interesting is the ensemble gradient flux approach, in which multiple gradient flux methods, each with different assumptions or limitations, can be combined to produce a much more temporally complete and informative dataset. Rather than selecting a single method, the ensemble approach is intriguing and, in my opinion, a good way to approach gappy data.
Author response
We agree that the ensemble is one of the most useful outcomes for users of gappy records. In the discussion paper, we initially used the mean ensemble. Now we use the inverse-RMSE weights. In the revision, the ensemble will be described as implemented. We will also release the spread between pairs for every 30-min value, so users can see where pairs disagree. The spread ranks uncertainty: the absolute difference between the ensemble and EC increases with it (Spearman ρ = 0.43 for CO₂ and 0.58 for H₂O). However, it is not a calibrated interval, and we will say so. With the extended record, the ensemble produced values in a median of 1.6 (CO₂) and 1.5 (H₂O) times as many 30-min windows as the best single pair at each site.
Changes in the manuscript
Sect. 2.8: “For each ensemble value, we report the number of contributing pairs and values, the contributing approaches, and the spread between pairs, defined as the standard deviation of the 30-min mean GFs across contributing approach and sampling-height-pair combinations (undefined when one combination contributes). The spread ranks uncertainty: the absolute difference between the ensemble and EC increased with the spread (Spearman ρ = 0.43 for CO₂ and 0.58 for H₂O). However, it is not a calibrated interval because EC fell within ±1 spread of the ensemble in only 48–51 % of windows and within ±2 spreads in 72 %.”
Sect. 5: “The revised release (August 2021 to August 2026) adds, for each site, a 30-min ensemble table ({SITE}_GF_ENSEMBLE_30min.csv) with the inverse-RMSE-weighted ensemble flux, the spread between pairs, the number of contributing pairs and values, the contributing approaches, and the matching EC flux, with times in local time and UTC. The table lists every half-hour of the record with a flag (ok, single_pair or no_reliable_GF); no gap filling is applied.”
1.2 Abbreviations and acronyms
The analyses and figures in this manuscript are excellent, and I am impressed with the figures that synthesize huge amounts of information and add conceptual or diagrammatic elements that help with a more intuitive understanding of the information. I did struggle with the many abbreviations and acronyms for the flux methods, statistical methods, sensor placement, particularly when a sentence had multiple acronyms. I am not sure the best way around this, but two ideas are to either drop some of the acronyms and write out the term, or include a glossary table. My recommendation is to publish this manuscript with minor revisions.
Author response
Thank you. We will do both. A glossary of abbreviations and symbols will become the new Table 1 and will be referred to at the end of the Introduction. We will also write terms out where several acronyms meet in one sentence, for example in the description of the AE (see 1.3) and in Sect. 3.2. Tables after Table 1 will be renumbered by one.
Changes in the manuscript
New Table 1 (Abbreviations and symbols): Defines AA and AW, AE, CCC, d, EC, GF, h, hc, K, L, MAD, MBR, RSHP, SNR and TSNR, u*, ū, WP, zA and zB, z, z0, ζ and φh, with units.
Sect. 1, last sentence: “Abbreviations and symbols are defined in Table 1.”
1.3 Line 100
‘With the AE, K is derived from wind and turbulence statistics’ – this was one early sentence I struggled with
Author response
We will rewrite the sentence to say what K depends on. We will also rewrite the following sentence on the WP so that the difference between the two methods is clear.
Changes in the manuscript
Sect. 1: “With the AE, K is calculated from the friction velocity (u*) measured by the sonic anemometer, the height of the measurement above the zero-plane displacement height, and a correction for atmospheric stability (Sect. 2.4.2). […] The WP approach simplifies the AE by deriving K from mean wind speed and surface roughness length rather than from u*.”
1.4 Lines 129–130
I don’t fully understand this. I assume that for gradient flux the gas concentrations are measured at multiple heights (between 4-8 sampling levels) while eddy covariance is above the canopy. I am therefore confused which sensors are both above the canopy (AA) and which sensor would be within the canopy for (AW). Since gradient has multiple sensors, while eddy covariance is one sensor I assume that the eddy covariance measure is always above (A), but I don’t understand what the other paired measurement is. In figure 1B which category do JORN, KONZ, GUAN, and HARV fall into? Are JORN and KONZ both AA because there is not enough canopy to be ‘within’ while GUAN and HARV are AW because some instruments are above and some within?
Author response
The referee's later reading (comment 1.12) is correct. AA and AW describe the two profile inlets used for one GF estimate, not the EC system, which is always at the tower top. Every GF pair has at least one inlet above the mean canopy height (hc), so no pair is entirely within the canopy. A pair is AA when its lower inlet zA is above hc, and AW when zA is at or below hc. JORN and KONZ (hc ≤ 1 m) have only AA pairs because every inlet is above the short canopy. GUAN and HARV have both AA and AW pairs. We will state this directly in Sect. 2.1. We will also add Figure A20, which shows every inlet relative to the canopy at each tower.
Changes in the manuscript
Sect. 2.1: “The categories AA and AW describe the two profile inlets used for a GF estimate, not the EC system, which is always at the tower top. Every pair includes at least one inlet above the mean canopy height (hc), so no pair has both inlets within the canopy. We write zA for the lower and zB for the upper inlet of a pair; a pair is AA when zA > hc and AW when zA ≤ hc. […] Of the sites in Fig. 1b, JORN and KONZ (hc ≤ 1 m) have only AA pairs because every inlet is above the short canopy, whereas GUAN and HARV have both AA and AW pairs.”
New Fig. A20 caption: “Inlet heights at each NEON tower (points; above-canopy inlets in blue, within-canopy inlets in red) and the mean canopy height (green bars), with sites ordered by canopy height (square-root scale). Every tower has at least two inlets above the canopy, short towers have only above-canopy inlets, and tall forests have one to three inlets within the canopy; the top inlet sits at the EC measurement height.”
1.5 Section 2.4.1
In this dataset are the target gas and tracer different combinations of CO₂ and H₂O? For example to calculate CO₂ flux with MBR, the H₂O is the tracer, and to calculate H₂O flux, CO₂ is the tracer?
Author response
Yes. CO₂ fluxes use H₂O as the tracer, and H₂O fluxes use CO₂. Both tracer fluxes are the 30-min EC values interpolated to the midpoint of the paired sampling window. We will state this in Sect. 2.4.1.
Changes in the manuscript
Sect. 2.4.1: “In this dataset, the tracer is the other gas measured on the same profile: CO₂ fluxes are calculated with H₂O as the tracer (F1 is the EC water vapor flux derived from the latent heat flux, and ΔC1 is the H₂O mixing-ratio difference between the same two inlets), and H₂O fluxes with CO₂ as the tracer (F1 is the EC turbulent CO₂ flux, and ΔC1 is the CO₂ difference). F1 is the 30-min EC flux interpolated to the midpoint of the paired sampling window (Sect. 2.2).”
1.6 Figure 2
This figure is helpful! A very minor aesthetic comment, but I suggest aligning the edge of the Q3 box with the Q2 box to show that they both emerge from Q1. I presume the idea is that Q2 and Q3 are based on the outcomes of Q1?
Author response
Yes, both Q2 and Q3 use the reliable pairs identified under Q1. We will extend the Q3 (Ecosystem Assessment) box so that its right edge aligns with the Q2 box. We will also list the filtering steps in the caption, including the new spike removal.
Changes in the manuscript
Fig. 2: Q3 box aligned with the Q2 box. Caption: “Filtering includes u*, SNR, and TSNR screening and despiking (Sect. 2.5).”
1.7 Line 357
I don’t understand how the ability of GF to capture different ecosystem components was evaluated. Does the footprint of EC and GF differ? Why would GF capture different ecosystem components from EC?
Author response
We did not evaluate this, and we will remove the sentence. The footprints do differ. Every GF represents the layer between two inlets below the EC sensor, so its footprint is smaller and closer to the tower than the EC footprint. We will describe this in the limitations, where we note that the mismatch can only lower agreement, so the reported CCC values are conservative.
Changes in the manuscript
Sect. 3, opening (replaces the sentence on ecosystem components): “The challenge was to determine which sampling height pairs represented the ecosystem-scale flux measured by EC, which we assessed by their concordance with the EC method.”
Sect. 4.4: “First, the profile does not extend above the EC system: the top inlet is at the tower top, close to the EC measurement height, so every GF represents a layer below the EC sensor, with a smaller and closer footprint. This mismatch can only reduce agreement, so the CCC values are conservative […].”
1.8 Figure 3
It would be helpful to remind readers in the caption that CCC varies from 1 to -1 which represent strong concordance, while values close to 0 represent weak concordance.
Author response
We will add this to the caption.
Changes in the manuscript
Fig. 3 caption: “CCC ranges from −1 to 1: values near 1 indicate strong agreement, values near 0 weak agreement, and negative values GF that vary opposite to EC.”
1.9 Line 408
very minor but the font size changed in this section
Author response
Thank you. The first paragraph of Sect. 3.2 carried a different paragraph style. We will rewrite it in the body style.
Changes in the manuscript
Sect. 3.2, first paragraph: rewritten in the body text style.
1.10 Figure 6F
it’s interesting that TOOL is the only site with low canopy and CCC<0.5. Do the authors have any thoughts on why this site has lower concordance, given the low canopy?
Author response
We examined TOOL in detail. With the two additional years, the corrected MBR estimator and spike removal, three MBR CO₂ pairs at TOOL now pass the threshold (CCC = 0.56–0.76). Its AE and WP CO₂ pairs remain poor (CCC ≤ 0.16). Only 54–531 CO₂ values per pair pass the filters, compared with up to 4,174 for H₂O. CO₂ fluxes at this Arctic tundra site are small for most of the year, so the CO₂ gradient rarely exceeds three times its noise. In addition, over a very smooth, often snow-covered surface, the K estimates of Eqs. 3 and 7 are sensitive to small errors in d and z0. NIWO, an alpine tundra site, shows the same pattern. We conclude that where fluxes are small, a low signal-to-noise ratio limits the GF method more than canopy structure does. We will add a paragraph on this to Sect. 4.1.
Changes in the manuscript
Sect. 4.1 (new paragraph): “A short canopy alone is not sufficient. At the Arctic tundra site TOOL (hc = 0.04 m), H₂O pairs were reliable, but for CO₂ three MBR pairs passed the threshold (CCC = 0.56–0.76), while the AE and WP remained poor (CCC ≤ 0.16), and only 54–531 CO₂ values per pair passed the filters, compared with up to 4,174 for H₂O. CO₂ fluxes at TOOL were small for most of the year, so CO₂ gradients rarely exceeded three times their noise, and over a very smooth, often snow-covered surface the K estimates of Eqs. (3) and (7) were likely sensitive to small errors in d and z0. NIWO, an alpine tundra site, similarly had nine reliable H₂O pairs but only one reliable CO₂ pair (AE, CCC = 0.70). Where fluxes were small, low SNR rather than canopy structure limited the GF method.”
1.11 Figure 8
the figure appears to be missing from the main text, I can only see the caption.
Author response
The figure appears to have been dropped when the discussion paper was produced. We will retain Figure 8 with its caption, regenerate it from the 2021–2026 data, and check that it appears in the typeset version. In regenerating it, we found that the GF and EC axis titles of panel (a) had been swapped in the plotting code; this will be corrected. GUAN will no longer be shown in panels (b) and (c) because it has no reliable pair.
Changes in the manuscript
Fig. 8: Figure regenerated from the 2021–2026 data; GF and EC axis titles of panel (a) corrected. Caption: “GUAN (Fig. 1b) is not shown because it has no reliable sampling height pair.”
1.12 Lines 515–525
I think I understand the AA and AW terminology better from this section. How I understand it now is that AA and AW is not refering to EC vs GF, it’s about the paired sample height comparisons for a GF estimate. AA and AW is referring to pairs of sampling heights for gradient flux calculations where the paired gradient measurements are either both above or above and within the canopy.
Author response
That is correct. This definition will appear where AA and AW are first introduced (Sect. 2.1) and in the glossary (Table 1); see our response to 1.4.
Changes in the manuscript
Table 1, AA and AW: “Canopy level of a sampling height pair: both inlets above the mean canopy height (AA), or the upper inlet above and the lower inlet at or below it (AW); + and − mark the inlet closest to the canopy top from above and below (Table 2).”
1.13 Lines 554–555
It’s unclear to me how GF in itself helps with decomposition into source contributions? Is this because below canopy gradient heights can represent primarily respiration/evaporation contributions, while above-canopy would include both respiration and photosynthesis/transpiration? I am not sure that the analysis in this paper particularly helps highlight that possibility with GF and found the statements about partitioning here, and above (line 357), more confusing than helpful.
Author response
We agree that our analysis does not show this. We will remove both statements (here and at line 357; see 1.7) and say explicitly that this was not evaluated in the work presented. Source partitioning will remain only as future work in the Conclusions.
Changes in the manuscript
Sect. 4.4 (replaces the statement on decomposing fluxes by source): “We did not evaluate whether within-canopy pairs can separate soil and canopy contributions; doing so would require independent soil or understory flux measurements.”
1.14 Figures A11 and 6
Figure A11 and Figure 6 look almost like the same figure to me, though there seems to be some slight difference in the underlying data as panels C, D, E, F are slightly different.
Author response
Figure A11 is the H₂O version of Figure 6 (CO₂). Panel (a) is the same schematic in both, which made them look alike. The caption of Figure A11 will begin “As Fig. 6, for H₂O” and will state how the H₂O results differ. While checking the two figures, we found that the x-axis title of panel (c) read upper sensor height. It should read lower inlet height minus canopy height, and it will be corrected in both figures.
Changes in the manuscript
Fig. A11 caption: “As Fig. 6, for H₂O: (a) schematic of AA and AW pairs; (b) CCC by canopy level and approach; (c) mean CCC against the height of the lower inlet above the canopy (zA − hc); (d) mean CCC against mean canopy height; (e) correlations between CCC and canopy and geometry metrics; (f) highest CCC at each site against hc, colored by VAI. […] The patterns match those for CO₂, but more H₂O pairs exceed 0.5, including AW pairs of the AE and WP in tall canopies (Sect. 3.2).”
1.15 Figure A13
some of the ranges of flux values seem very large eg: CPER with a values at -1000, but even other sites with axis values at +/- 200.
Author response
These values came from GFs computed when the tracer gradient in the MBR approached zero, or when K in the AE and WP was large; filtering by signal-to-noise ratio did not remove all of them. We made two changes. First, the MBR flux is now the median of bootstrap draws, which is much less sensitive to a near-zero tracer gradient (Sect. 2.4.1). Second, we added a spike-removal step (Sect. 2.5) and a new appendix that describes where and why spikes occur (Appendix B).
Spike removal took out 2.4 % of the filtered values. For the MBR CO₂ flux at CPER, the range fell from −23 to 828 µmol m−2 s−1 after filtering to −22 to 26 after spike removal. A few isolated ensemble values remain, for example up to 133 µmol m−2 s−1 at WREF, against a 99.5th percentile of 75. They still stretch some axes, and the revised caption will say so.
Changes in the manuscript
Sect. 2.5 (new paragraph): “Filtered GFs still contained isolated extreme values (up to −4,700 µmol m−2 s−1 for the AE and 830 µmol m−2 s−1 for the MBR CO₂ flux at CPER) produced when K was large or the tracer gradient in Eq. (2) approached zero. We therefore removed values outside the median ± 6 MAD (scaled MAD, 1.4826 × median absolute deviation) of the GF, computed separately for each site, approach, gas, sampling height pair, time of day (day, 07:00–18:59 local time; night otherwise) and season; strata with fewer than 20 values were left unchanged, and EC data were not used. Despiking removed 2.4 % of the approximately 5 million filtered values. […] Across sites, the median of the largest absolute MBR CO₂ flux fell from 473 µmol m−2 s−1 after filtering to 104 µmol m−2 s−1 after despiking (at CPER, from a range of −23 to 828 to a range of −22 to 26 µmol m−2 s−1).”
Fig. A13 caption: “After despiking (Sect. 2.5), the extreme values present before despiking (e.g., up to about 830 µmol m−2 s−1 for the MBR at CPER, from near-zero tracer gradients) are removed, and the bulk of the GF and EC distributions overlap; a few isolated GF values remain and stretch some axes (e.g., WREF, up to 133 µmol m−2 s−1, against a 99.5th percentile of 75).”
New Appendix B: “Gradient-flux spikes: frequency, conditions and causes” (spike rates by approach, gas and canopy level; conditions under which spikes occur; mechanisms; comparison with EC).
1.16 Figure A15
I was curious to see the relationship for TOOL but it’s missing entirely from this graph and A13. It does appear in the H₂O fluxes so perhaps all CO₂ flux measurements were lost to filtering.
Author response
That is right for the discussion paper. TOOL had no CO₂ pair with CCC ≥ 0.5, so it had no CO₂ ensemble. In the revised record, TOOL has three reliable MBR pairs and 887 half-hourly CO₂ ensemble values, so it will appear in Figures A13 and A15. See also our response to 1.10.
Changes in the manuscript
Fig. A15 caption: “Sites without a reliable CO₂ pair have empty panels; TOOL has 887 half-hours.”
Citation: https://doi.org/10.5194/essd-2026-413-AC1 -
AC2: 'Comment on essd-2026-413 AC2', Sparkle Malone, 06 Nov 2026
Authors Response to RC2
We thank Referee 2 for the constructive review and for focusing on what dataset users need. We addressed most comments with new analyses, figures, and data, not only with text. The new material comprises:
- Figure A20: inlet heights relative to the canopy at every tower.
- Figure A21 and Appendix D: seasonal day and night concentration profiles at all 47 sites.
- Figure A22: seasonal change in the aerodynamic canopy height.
- New paragraphs in Sects. 2.1, 2.2, 2.4, 4.2, 4.3 and 4.4.
- A 30-min ensemble table with the spread between pairs in the data release.
Overview of the planned revision
Below, each referee comment is quoted in full (in italics), followed by our response and by the changes we will make to the manuscript (shaded), quoted from our revised draft. Section, figure, table, and equation numbers refer to the revised manuscript. Line numbers will be given with the tracked-changes manuscript submitted after the discussion.
While addressing both referees' comments, we extended the record, corrected and documented the processing, and reran every calculation. The revised manuscript will include five changes that affect most of the numbers in the paper:
Record extended from August 2021–July 2024 to August 2021–August 2026
Two more years at every site
Reliable sampling height pairs (RSHPs) at 42 sites for CO₂ and 43 for H₂O; 42 sites have them for both gases, compared with 30 in the discussion paper
Spike removal: values outside the median ± 6 scaled MAD per site, approach, gas, pair, day or night, and season
Isolated extreme values from near-zero tracer gradients or large eddy diffusivities (Referee 1, Fig. A13)
2.4 % of about 5 million filtered values removed; median site maximum of |MBR CO₂| from 473 to 104 µmol m−2 s−1 (Sect. 2.5; new Appendix B)
Processing corrections (listed below)
Text and code disagreed in places, found while documenting Eqs. 2–7 for Referee 2 and in an independent code review by co-author J. D. Gewirtzman
All fluxes and evaluations recalculated; main conclusions hold, and the secondary conclusions that changed are listed below
Ensemble described as implemented (inverse-RMSE weighted) and released with the spread between pairs
The text said the mean ensemble was used; Referee 1 noted the value of the ensemble for gappy data
New 30-min ensemble table in the data release with the spread, the number of pairs and the approaches for every half-hour (Sects. 2.8 and 5)
Mean canopy height hc taken from the NEON airborne lidar canopy height model throughout
The discussion paper mixed lidar and NEON site-metadata heights
Site canopy heights span 0.03–38 m (previously 0.15–53 m), consistently across the AA/AW categories, the random forest and the placement analysis (Sect. 2.1; Table 1)
Processing corrections. We will describe the methods as implemented after the following corrections:
- MBR (Eq. 2): the target flux is the median of 1000 bootstrap draws of F1ΔC2/ΔC1, because the mean of a ratio is unstable when ΔC1 is close to zero; the 2.5–97.5 percentile range of the draws is the uncertainty (Sect. 2.4.1).
- Aerodynamic canopy height h: obtained by inverting the logarithmic wind profile at the top inlet, using the wind-measurement height, with d = 0.66h, z0 = 0.1h, and a centered 20-day window (Sect. 2.4.2).
- Effective height: z is the geometric mean of each pair's two inlet heights; ζ = (z − d)/L is computed for each pair; k = 0.4 throughout (Eq. 6 used 0.41).
- Wind-profile method (Eq. 7): Kwp is set to missing when z − d ≤ z0, where the equation is singular (Sect. 2.4.3).
- Ensemble: weights use the direct root-mean-square difference between each pair's GF and EC, and each GF value is assigned to the half-hour containing the midpoint between its two inlet samples, the time to which EC is interpolated (Sect. 2.8).
- Diel curves: absolute diel values and daily totals in the discussion paper were half their correct values because of a unit conversion error; they are now in g m−2 h−1 (Sect. 2.9; Table A2). Correlations, slopes, relative biases and signs were not affected.
Key results, discussion paper to revised manuscript.
Sites with an RSHP: CO₂ / both gases
38 / 30
42 / 42 (43 for H₂O)
Reliable height-pair combinations: CO₂ / H₂O
11% / 19%
16% / 23%
RSHPs per site: CO₂ / H₂O
3.8 / 5.7
5.1 / 7.1
Mean CCC for CO₂, MBR / AE / WP
0.27 / 0.09 / 0.13
0.34 / 0.14 / 0.15
Share of CO₂ RSHPs that are AW
28%
38%
Random forest accuracy
95% (no-information rate 85%)
89.6% on held-out pairs (no-information rate 80.8%)
Ensemble used; mean site CCC across the six methods
mean; 0.65–0.66
inverse-RMSE weighted; 0.66–0.68
Sites with seasonal diel curves: CO₂ / H₂O
31 / 43
40 / 42
R² of GF vs EC diel curves
0.58–0.99
median 0.91 (CO₂), 0.99 (H₂O)
The main conclusions are unchanged: the MBR is the most reliable approach, concordance declines with canopy height and complexity, and the ensemble reproduces the seasonal diel cycles measured by EC. Three secondary conclusions changed and will be corrected. (1) The best AE and WP lower inlets in tall forests lie in the upper crown, not in the trunk space. (2) The statement that AW configurations raise AE and WP concordance from about 0.18 to 0.64 in tall forests did not hold and will be removed. (3) The random forest is evaluated on held-out pairs rather than in-sample.
The revision will also add a glossary (new Table 1, so later tables are renumbered by one), three appendix figures (A20 inlet heights relative to the canopy, A21 seasonal concentration profiles, A22 seasonal aerodynamic canopy height) and three appendices: B on gradient-flux spikes, C on selecting pairs where no EC benchmark exists (relevant for CH₄), and D on seasonal day and night CO₂ and H₂O concentration profiles at every site. Panel labels will be lowercase, following ESSD style. Our responses to Referee 1 are posted separately.
Responses to RC 2
2.0 General comments
The study presented a profile-based CO₂ and H₂O flux dataset derived from NEON datasets, along with collocated eddy-covariance fluxes, concentration profiles, tower geometry, and canopy structural metrics. It tested and summarized the performance of three gradient flux approaches across 47 NEON sites and proposed an ensemble data product for further research.
I think the research initiative, generated datasets, and the paper itself are important in many ways, both for research and application. In particular, I appreciate the authors’ effort and comprehension in designing, implementing, and discussing the study and its outcome. I would recommend the paper for consideration for publication in Earth System Science Data after careful revision. Here, I provide some general and specific comments to help improve the presentation and, most importantly, the usability of the dataset in future research. Since this is a data-focused journal, I don’t insist that all comments be addressed directly with additional analyses. However, I suggest the authors consider adding these points to the discussion and/or providing the necessary data and metadata to enable further research.
Author response
Thank you. Where an analysis was possible with the NEON data, we carried it out; where it was not, we added the point to the discussion and released the data needed for others to do it. The data release will include the 30-min ensemble table with the spread between pairs, the site canopy attributes and the inlet heights of every tower (Sect. 5).
2.1 [1] Vertical canopy structure
As the authors pointed out, the goal of this paper is to test when and where these gradient-based flux estimates work, and why. I think it is equally important to present the estimated flux dataset and also metadata that explain the validity or robustness of the methods and dataset. In particular, I suggest including information or data on the canopy’s vertical structure. The extracted canopy metrics are valuable for exploratory analyses across sites, but they are limited if one intends to dive further into the details of each case and generalize best practices for future deployment. Instead, the vertical profiles of structures (e.g., LAI, biomass) would be much more informative. Figure 9 is great for its abstraction but still contains a lot of information. But looking at the results between short- and tall-canopy systems, I wonder whether it makes more sense to use the vertical structure profile to separate the canopy into crown and trunk space (esp. forests), potentially revising the groups in Table 1.
Author response
We agree that vertical structure is the missing piece for generalizing placement rules, and we did three things.
First, the new Figure A20 shows the height of every inlet relative to the mean canopy height at all 47 towers, so users can see which inlets sit in the crown and which sit below it.
Second, we re-examined the placement analysis behind Figure 9 with the corrected data and lidar canopy heights. In tall forests (hc≥ 15 m), the best MBR pairs had the lower inlet in the upper crown, a few meters below the canopy top (median zA − hc = −2.8 m for CO₂ and −2.6 m for H₂O). The best AE and WP lower inlets lay deeper (median −6.3 to −8.7 m), but still in the upper crown, not in the trunk space as the discussion paper stated. We will correct this in Sect. 4.2 and Figure 9.
Third, we will state as a limitation that our lidar metrics describe structure, not the vertical distribution of leaf area. We will also point users to the NEON discrete-return point clouds, from which leaf-area profiles can be derived.
We did not split the AW category (Table 2 in the revision, formerly Table 1) into crown and trunk space. Without leaf-area profiles at all sites, we cannot consistently place the crown base. The random forest already uses the continuous height of the lower inlet relative to the canopy (Relative Dist A, zA − hc), which a crown/trunk split would only turn into categories. This was one of the two most important predictors. For the MBR, the predicted probability of a reliable pair rose from about 0.13 with the lower inlet 30 m below the canopy top to about 0.35 at or above it (Fig. 5).
Changes in the manuscript
Sect. 4.2: “For tall canopies (hc ≥ 15 m), the MBR performed best with the lower inlet in the upper crown, a few meters below the canopy top (median zA − hc = −2.8 m for CO₂ and −2.6 m for H₂O; interquartile range about −7 to +4 m), and the upper inlet 14–16 m above the canopy (top-quartile CCC = 0.67 and 0.70). The best AE and WP pairs in tall canopies placed the lower inlet deeper in the crown (median −6.3 to −8.7 m), and the upper inlet 12–17 m above the canopy, but their top-quartile CCC remained low (0.22–0.43) […]. These positions are in the upper crown rather than the trunk space.”
Sect. 4.4: “Third, the canopy metrics describe structure, not the vertical distribution of leaf area or of sources and sinks. The seasonal concentration profiles (Fig. A21; Appendix D) are a first proxy for source and sink layers; vertical leaf-area profiles from the NEON discrete-return lidar point clouds (DP1.30003.001) would allow future analyses to relate each inlet to the height of peak foliage.”
New Fig. A20: Inlet heights at every NEON tower relative to the mean canopy height (see caption under 2.4).
2.2 [2] Vertical distribution of sources and sinks
Along these lines, I also suggest including information or data on the vertical distribution of sources/sinks for both CO₂ and H₂O. As a best guess, I assume mid-canopy (peak foliage) and soil to be two major source/sink levels. Diving deeper, this assumption depends on the target (CO₂ vs H₂O), time (day vs night, summer vs winter), and ecosystem (dense forest vs savanna/shrub). I think it would be valuable to generate and present monthly or seasonal average CO₂/H₂O concentration profiles and use them to inform the source/sink distribution. This could help explain performance differences among sites and between CO₂ and H₂O.
Author response
We computed seasonal mean day and night concentration profiles for both gases at all 47 sites. Four representative sites will appear in the new Figure A21, and all sites in every season in the new Appendix D (Figs. D1 and D2, Table D1).
The profiles partly support the referee's expectation. Across the 25 sites with hc ≥ 5 m, night-time CO₂ and daytime H₂O concentrations were highest at the lowest inlet at 96–100 % of sites in every season, so the soil and lower canopy are the dominant source. In contrast, leaf uptake rarely produced a distinct daytime CO₂ minimum inside the canopy (4–20 % of these sites, highest in summer). In the growing season, the second-highest inlet lay above the concentration extreme at 96–100 % of these sites. The uppermost pair therefore samples air above the main source or sink layer, while pairs with a near-ground inlet span it. We will use these results in Sect. 4.2 to explain the differences between gases and canopy heights.
Changes in the manuscript
Sect. 4.2: “The vertical distribution of sources and sinks helps explain the patterns. Seasonal mean concentration profiles (Fig. A21; all sites and seasons in Appendix D) show that nighttime CO₂ and daytime H₂O concentrations were highest at the lowest inlet at 96–100 % of the 25 sites with hc ≥ 5 m in every season, consistent with respiration and evaporation near the ground. In contrast, even in summer the daytime CO₂ minimum lay within the canopy at only 20 % of these sites, so leaf uptake did not create a distinct concentration minimum in most forests. […] For H₂O, with sources at the ground and throughout the crown, the gradient is steepest across the canopy, which favors cross-canopy pairs for the K-based methods. For CO₂, uptake by leaves and release from soil and stems partly cancel within the canopy, so the MBR, in which K cancels between tracer and target, benefits most from clearance above the crown.”
New Appendix D: “Seasonal CO₂ and H₂O concentration profiles”: method (D.1), where the sources and sinks are (D.2), Table D1 (percentage of sites where the concentration extreme is at the lowest inlet, and where the second-highest inlet lies above it, by gas, time of day, canopy class and season), and Figs. D1–D2 (profiles at all 47 sites).
2.3 [2] Section 2.1, item 1
Section 2.1: Consider adding a few more details to avoid confusion or for readers who are less familiar with NEON or gradient flux methodology: 1) I assume at least one of the sampling levels needs to be above the canopy for all pairs, since all categories begin with “A”. But I got confused when reading Figure 6C, which suggests cases where the top level is below canopy height;
Author response
The assumption is correct: every pair includes at least one inlet above the mean canopy height. The confusion came from an error in Figure 6c. Its x-axis was labeled “upper sensor height”, but it shows the height of the lower inlet relative to the canopy (zA − hc), which is negative for AW pairs. We will correct the label in Figures 6 and A11 and in their captions, and state the rule in Sect. 2.1.
Changes in the manuscript
Sect. 2.1: “Every pair includes at least one inlet above the mean canopy height (hc), so no pair has both inlets within the canopy.”
Fig. 6 caption: “(c) mean CCC against the height of the lower inlet above the canopy (zA − hc) by approach”.
2.4 [2] Section 2.1, item 2
2) Please define (quantitatively) the canopy heights and the relative location above/within/below the canopy heights. The categories in Table 1 are good for summary and grouping, but please provide the quantitative metrics. If possible, consider revising Figure 1B to show all pair categories.
Author response
Sect. 2.1 will define the canopy height and the categories quantitatively, with the range of inlet positions in each category. Rather than crowd Figure 1b, we added Figure A20, which shows every inlet at every tower relative to hc. Table 2 (formerly Table 1) will be cited where the categories are defined.
Changes in the manuscript
Sect. 2.1: “We write zA for the lower and zB for the upper inlet of a pair; a pair is AA when zA > hc and AW when zA ≤ hc. Relative to hc, the lower inlet was 0.1 to 29 m above the canopy in AA pairs and 0.3 to 37 m below it in AW pairs (Table 2; Fig. A20). […] The hcis the mean of the NEON airborne canopy height model within the tower footprint, averaged over the 2–7 flight years available at each site, and is treated as a fixed site attribute; at the four sites without lidar coverage (DCFS, KONA, LAJA and TREE) we used the NEON-reported site canopy height.”
New Fig. A20 caption: “Inlet heights at each NEON tower (points; above-canopy inlets in blue, within-canopy inlets in red) and the mean canopy height (green bars), with sites ordered by canopy height (square-root scale). Every tower has at least two inlets above the canopy, short towers have only above-canopy inlets, and tall forests have one to three inlets within the canopy; the top inlet sits at the EC measurement height.”
2.5 [2] Section 2.1, item 3
3) How does the study address ecosystems that experience changing canopy, either in height or foliage? When were the canopy metrics measured, and how often were they resampled?
Author response
The lidar metrics come from NEON airborne flights, two to seven per site, averaged into one value per site; they are fixed site attributes. Seasonal change enters the AE and WP through the aerodynamic canopy height. It is recomputed from turbulence statistics in a centered 20-day window and sets d and z0. The new Figure A22 shows that this height varied seasonally by a median of 28% of its site median, and by more than 25% at 30 sites, mostly short-canopy sites. We will note as a limitation that the AA/AW categories and the random forest predictors do not change with season.
Changes in the manuscript
Sect. 2.1: “Seasonal changes in canopy height and foliage enter the GF calculations through the aerodynamic canopy height (Sect. 2.4.2), which varied seasonally by a median of 28 % of its site median and by more than 25 % at 30 sites (Fig. A22).”
Sect. 4.4: “Second, the canopy categories and the lidar metrics are fixed site attributes. Seasonal changes in canopy height and foliage enter the AE and WP through the aerodynamic canopy height (Fig. A22), but not through the AA/AW categories or the random forest predictors.”
New Fig. A22 caption: “Monthly aerodynamic canopy height relative to its site median, for sites with (a) a mean canopy height < 5 m (orange) and (b) a mean canopy height ≥ 5 m (blue). Short canopies change height seasonally by up to an order of magnitude (median amplitude 28 % across all sites; > 25 % at 30 sites), whereas tall canopies are nearly constant; d and z0 in the AE and WP follow these changes.”
2.6 [3] Section 2.2, item 1
Section 2.2: Consider adding or clarifying a few more details: 1) Was each CO₂/H₂O sampling level sampled once (for 6 min) over each 30-min period?
Author response
No. Each inlet is sampled once per profile cycle in a 10-min slot, so each level is revisited every 40–80 min, depending on the number of levels. A level is therefore not sampled in every 30-min period. We will state this in Sect. 2.2.
Because the two inlets of a pair are not sampled at the same time, we also quantified how the time between them affects agreement with EC. Within a pair, agreement declined as the separation grew, mostly for the WP and MBR. Separations of more than 30 min were rare (5.5–7.0 % of values), so they did not determine which pairs were reliable. We will add this to Sect. 4.4.
Changes in the manuscript
Sect. 2.2: “Each inlet is sampled once per profile cycle in a 10-min slot (about 4 min of line flushing followed by a 6-min averaging period), so each level is revisited every 10 min × the number of levels: every 40 min at four-level towers, 50–60 min at five- and six-level towers, and 70–80 min at the two tallest towers. Therefore, a level is not sampled in every 30-min period.”
Sect. 4.4: “Within a sampling height pair, agreement with EC declined as the time between the two samples grew: compared with samples 10 min apart, the CCC of values 20–30 min apart was lower by a median of 0.07 for the MBR (both gases), 0.12 (CO₂) and 0.03 (H₂O) for the WP and 0.02 or less for the AE, and the CCC of WP values more than 30 min apart was lower by 0.25 (CO₂) and 0.11 (H₂O). Such long separations were rare (7.0 % of CO₂ and 5.5 % of H₂O values), so they did not determine which pairs were reliable.”
2.7 [3] Section 2.2, item 2
2) The description in Lines 169-170 is a bit unclear, particularly for data users who have no experience with this type of instrument;
Author response
We will rewrite the paragraph to explain how the two samples of a pair are chosen, how far apart they are, and how quality flags are applied. We also found a typo: “4-min flux” should have read “4-min flush”.
Changes in the manuscript
Sect. 2.2: “A GF value combines the 6-min means of the two inlets sampled closest in time, in either order, and pairs were dropped when either concentration carried a NEON final quality flag. The centers of the two samples were 10–60 min apart: 23–29 % of filtered values were 10 min apart, 66–70 % 20–30 min apart and 5.5–7.0 % more than 30 min apart (Sect. 4.4).”
2.8 [3] Section 2.2, item 3
3) Lines 174-175: It’s unclear where and which variable this alignment occurred;
Author response
Two alignments are applied. Meteorological variables are averaged over the combined window of the pair. EC fluxes and u* are interpolated to the midpoint between the two samples. Both will be named in Sect. 2.2.
Changes in the manuscript
Sect. 2.2: “Meteorological data were available as 1-min averages (2-min for horizontal wind speed), and each variable was averaged over the combined window of the pair, from the start of the first 6-min sample to the end of the second. […] The 30-min EC turbulent fluxes of CO₂, latent heat (converted to H₂O flux) and sensible heat, and the friction velocity were filtered and linearly interpolated to the midpoint between the two samples (Hall et al., 1999). A GF could therefore be calculated only when the EC half-hours before and after that midpoint were available.”
2.9 [3] Section 2.2, item 4
4) A side note: The NEON profile was designed to have the highest level at the eddy covariance heights, but not above, right? Would this limit or somehow bias this research’s conclusion or implication?
Author response
Yes, the top inlet is at the tower top, close to the EC measurement height. Every GF therefore represents a layer below the EC sensor, with a smaller and closer footprint than EC. This mismatch can only lower agreement, so the CCC values are conservative. The design also means that we could not test pairs with the upper inlet above the EC height. Such pairs may benefit tall forests, where the best upper inlets were 12–17 m above the canopy. We will add this to the limitations.
Changes in the manuscript
Sect. 4.4: “Three features of the NEON design also limit what the benchmark can show. First, the profile does not extend above the EC system: the top inlet is at the tower top, close to the EC measurement height, so every GF represents a layer below the EC sensor, with a smaller and closer footprint. This mismatch can only reduce agreement, so the CCC values are conservative, and we could not test configurations with the upper inlet above the EC height, which may benefit tall forests.”
2.10 [4] Section 2.4.1
Section 2.4.1: Please specify the known flux and known gas gradient.
Author response
We will specify both. In documenting Eq. 2, we also found that a single ratio of means is unstable when the tracer gradient is close to zero. The revised MBR flux is therefore the median of bootstrap draws, and the spread of the draws gives its uncertainty.
Changes in the manuscript
Sect. 2.4.1: “Because the mean of a ratio is unstable when ΔC1 is close to zero, we computed F2 as the median of 1000 bootstrap draws of the right-hand side of Eq. (2), using a normal distribution with the measured mean and standard deviation for each profile concentration. The 2.5–97.5 percentile range represents the uncertainty. […] In this dataset, the tracer is the other gas measured on the same profile: CO₂ fluxes are calculated with H₂O as the tracer (F1 is the EC water vapor flux derived from the latent heat flux, and ΔC1 is the H₂O mixing-ratio difference between the same two inlets), and H₂O fluxes with CO₂ as the tracer (F1 is the EC turbulent CO₂ flux, and ΔC1 is the CO₂ difference). F1 is the 30-min EC flux interpolated to the midpoint of the paired sampling window (Sect. 2.2).”
2.11 [5] Section 2.4.2
Section 2.4.2: Please specify the data sources used in equation 6.
Author response
We will specify every input. In doing so, we found that the discussion paper wrongly gave the temperature as measured at the lowest level; it is measured at the tower top. We will also state the stability range for which Eqs. 4–5 were derived, and how few values fall outside it.
Changes in the manuscript
Sect. 2.4.2: “where T is the air temperature (K) measured at the tower top, u* is the EC friction velocity, k = 0.4, g is gravitational acceleration (9.81 m s−2), and the virtual potential temperature flux is computed from the EC sensible (H) and latent (LE) heat fluxes as H/(ρcp) + 0.61 T LE/(ρλ), where ρ is the air density from the mid-tower pressure and the tower-top temperature, cp the specific heat of air, and λ the latent heat of vaporization. All EC quantities are 30-min values interpolated to the midpoint of the paired sampling window (Sect. 2.2).”
Sect. 2.4.2, after Eqs. 4–5: “We derived these forms for approximately −5 < ζ < 0 and 0 ≤ ζ < 1; we applied them to all values, and fewer than 0.2 % of the filtered AE and 0.4 % of the filtered WP values fell outside these ranges.”
2.12 [6] Section 2.4.3
Section 2.4.3: Please specify the derivation of d. Also, I assume it should be u measured at height z. Please be specific.
Author response
We will describe the derivation of h, from which d = 0.66h and z0 = 0.1h follow (Sect. 2.4.2). The wind speed in Eq. 7 is the mean horizontal wind speed at the upper inlet of the pair, not at z, because there is no anemometer at z. We will state this and the bias it implies.
Documenting Eqs. 3–7 also led us to corrections that are applied throughout the revision. First, the effective height now uses the geometric mean of each pair's two inlet heights, as the text stated. Second, the canopy-height inversion uses the height at which the wind was measured. Third, ζ is computed for each pair. Fourth, k = 0.4 everywhere. Finally, Kwp is set to missing when z − d ≤ z0, where Eq. 7 is singular. Near-singular values that remain (1 < (z − d)/z0 ≤ 1.5) account for 18–19 % of the WP spikes described in Appendix B.
Changes in the manuscript
Sect. 2.4.2: “Canopy height (h) was calculated aerodynamically by inverting the logarithmic wind profile at the top inlet (Pennypacker and Baldocchi, 2016), using u*, the mean wind speed and the height of the top inlet, with d = 0.66h and z0 = 0.1h, in near-neutral conditions (|z/L| < 0.1) with 0.35 < u* < 0.5 m s−1 and wind speed > 1 m s−1; daily medians were averaged over a centered 20-day window and interpolated to each pair window, so that d and z0 follow seasonal changes in canopy height and snow cover (Fig. A22).”
Sect. 2.4.3: “where k is the von Kármán constant, z − d is the effective height (Sect. 2.4.2) with z the geometric mean of the two inlet heights, ū is the mean horizontal wind speed measured by the 2-D sonic anemometer at the upper inlet of the pair, averaged over the pair window, and z0 is the roughness length estimated as 0.1h. Because ū is measured at the upper inlet rather than at z, Kwp is biased high when the wind profile is logarithmic, by the ratio ln((zB − d)/z0)/ln((z − d)/z0). Equation (7) is singular when z − d approaches z0, because the logarithm tends to zero; we therefore set Kwp to be missing when z − d ≤ z0 (using z − d without the 0.1 m floor), which removes pairs whose geometric-mean height lies close to the displacement height.”
2.13 [7] Section 2.7
Section 2.7: Most of these metrics focus on the structures, and are more useful in explaining the site-to-site difference. As noted above, I think some information on vertical structure and source/sink profiles may help this and future analyses. For example, what is the relative height of the 2nd sampling level compared with the peak LAI height? Or, what is the relative height of the 2nd sampling level compared to the maximal source/sink in the profile? Following this line of discussion, the potential seasonal changes in structure (e.g., canopy height, foliage) and in source/sink between day and night (e.g., photosynthesis vs respiration) may need to be considered.
Author response
We computed the second of these quantities for every site, season, and time of day, defining the 2nd sampling level as the second-highest inlet (Appendix D, Table D1). In summer, it was above the height of the largest concentration departure from the top inlet at 96–100% of sites with hc ≥ 5 m and at 73–100% of short-canopy sites (lowest for daytime CO₂). The uppermost pair therefore samples above the main source or sink layer, while pairs with a near-ground inlet span it. In winter, when leaves take up little CO₂, the daytime CO₂ minimum lies at or above the canopy top at 96 % of the taller sites. We could not compute the height of peak LAI with the metrics in the dataset (see 2.1).
Seasonal change in structure enters through the aerodynamic canopy height (see 2.5). Day–night differences now enter the processing through the spike removal, which is stratified by day and night. Spikes were more frequent at night: a night value was 1.4 (MBR), 2.0 (AE), and 3.5 (WP) times as likely to be a spike as a daytime one (Appendix B). The remaining static predictors will be discussed as a limitation (Sect. 4.4).
Changes in the manuscript
Sect. 4.2: “For nighttime CO₂ and daytime H₂O in every season, and for daytime CO₂ in spring and summer, the second-highest inlet was above the height of the concentration extreme at 96–100 % of these sites, so the uppermost pair samples air above the dominant source or sink layer, whereas pairs whose lower inlet is near the ground span it. In winter, when leaves take up little CO₂, the daytime CO₂ minimum lies at or above the canopy top at 96 % of these sites (Appendix D).”
Appendix B, B.3: “Holding other conditions fixed, a night GF was 1.4 (MBR), 2.0 (AE), and 3.5 (WP) times as likely to be a spike as a daytime one.”
2.14 [8] Figure 8
Figure 8 seems to be missing in the review version.
Author response
The figure appears to have been dropped when the discussion paper was produced. We will retain Figure 8 with its caption, regenerate it from the 2021–2026 data, and check that it appears in the typeset version. In regenerating it, we found that the GF and EC axis titles of panel (a) had been swapped in the plotting code; this will be corrected.
Changes in the manuscript
Fig. 8: regenerated from the 2021–2026 data; axis titles of panel (a) corrected.
2.15 [9] Discussion, Lines 488–497
Discussion paragraph (Line 488-497). This is a great discussion that provides context and background for this research and dataset. I suggest moving some of these points up front to help set expectations.
Author response
We will move the main points to the end of the Introduction and state the expectations that follow. The opening paragraph of the Discussion will refer to them instead of repeating them.
Changes in the manuscript
Sect. 1: “We aim not to test whether GF methods can replace EC everywhere, but to identify where and when they provide operationally useful ecosystem-scale fluxes, so profile systems on existing towers can be used for gases for which fast-response EC instrumentation remains expensive, technically demanding, or unavailable. […] We expect agreement with EC to be strongest in short, structurally simple canopies and to decline in tall, complex forests, where the roughness sublayer and vertically distributed sources and sinks violate the assumptions of flux-gradient theory, and we expect that combining approaches and height pairs would improve the temporal coverage of otherwise sparse GF records.”
Sect. 4, opening: “These results support the expectations set out in Sect. 1: profile systems on existing towers are a practical route to ecosystem fluxes of gases that are difficult to measure by EC, provided they are deployed where the structural conditions identified here are met.”
2.16 [10] Section 4.2, item 1
Section 4.2: 1) Please double-check and consolidate the definition of Ha and Hb here and in Figure 6.
Author response
We will use zA for the lower and zB for the upper inlet of a pair throughout the text, Figures 6, 9 and A11, Table 4 and the glossary. In checking this, we found two inconsistencies. First, the x-axis of Figure 6c called zA the upper sensor. Second, Table 4 labeled both lidar distance predictors “Relative Dist A” and described A as the upper inlet. Relative Dist A is in fact the lower inlet and Relative Dist B the upper inlet. The data files follow a different convention (in the pair label dLevelsAminusB, A is the upper inlet), and we will state this in Sect. 5.
Changes in the manuscript
Table 1: “zA, zB: Heights of the lower (A) and upper (B) inlet of a sampling height pair.”
Table 4: “Lower inlet relative to the canopy (Relative Dist A): Height of the lower inlet minus the mean canopy height (zA − hc). Upper inlet relative to the canopy (Relative Dist B): Height of the upper inlet minus the mean canopy height (zB − hc).”
Sect. 5: “In the EDI files, level A of the pair label dLevelsAminusB is the upper inlet; in this paper, zA denotes the lower inlet.”
2.17 [10] Section 4.2, item 2
2) This is one of the places where I think additional information on the vertical profiles (structure or source/sink) would help the discussion about the distinct patterns between short and tall canopies.
Author response
Sect. 4.2 will use the concentration profiles (Fig. A21; Appendix D) to explain the differences between gases and canopy heights; see our response and the proposed text under 2.2.
Changes in the manuscript
Sect. 4.2: see the text quoted under 2.2 and 2.13.
2.18 [10] Section 4.2, item 3
3) I suggest the author discuss the implications of the findings at short vegetation sites, particularly for those planning to apply the method somewhere else. The 5-10 m measurement height above short vegetation may not always be logistically feasible and may suffer from limited fetch.
Author response
We will add a paragraph on this to Sect. 4.2.
Changes in the manuscript
Sect. 4.2 (new paragraph): “For short vegetation, these results suggest that inlets do not need to be high. The best pairs had the lower inlet within 0.5 m of the canopy top and the upper inlet 6.5–8 m above it, but within-site concordance increased with inlet separation at most sites (Sect. 3.2), and the AA pairs at short NEON towers span separations from about 1 to 9 m (Fig. A20). Where a mast of 5–10 m is impractical, or where the upwind fetch of homogeneous vegetation is limited, a lower upper inlet reduces the gradient, and therefore the SNR, so fewer values will pass the filters. In that case, we recommend favoring the MBR, which uses the ratio of two gradients rather than an absolute K, and placing the upper inlet as high as the fetch allows; a common rule of thumb is a fetch of about 100 times the measurement height above d.”
2.19 [11] Section 4.3
Section 4.3: Reading the comparison between EC and gradient flux, I wonder whether the authors could also comment on whether it is feasible or practically useful to use flux gradient (e.g., slow-response sensors) as a complementary measurement at EC flux towers, e.g., for benchmarking, QA/QC, or gap-filling.
Author response
We will add a paragraph to Sect. 4.3. In short, GF is useful for benchmarking and quality control. For gap-filling, it helps only partly, because every GF needs EC turbulence data. It is best suited to adding gases measured by slow-response analyzers on the same profile.
For such gases, pairs cannot be selected against EC. The new Appendix C therefore tests two selection rules on CO₂ and H₂O, using EC only to score them. Where another gas on the same inlets has EC, as CO₂ does at NEON, the pairs that are reliable for that gas gave H₂O ensembles with a median site CCC of 0.66, close to the 0.68 obtained when pairs are selected on H₂O itself. Elsewhere, an ensemble of above-canopy MBR pairs can serve as the reference for screening and rescaling AE and WP pairs.
Changes in the manuscript
Sect. 4.3 (new paragraph): “The comparison also shows how GF can complement EC at existing flux towers. First, it can serve as a benchmark and quality-control tool. Because the MBR uses the EC flux of the tracer gas, persistent disagreement between the MBR flux and EC for the other gas, or between GF pairs, flags problems in either system, and the between-pair spread published with the ensemble supports such checks at 30-min resolution. Second, it can help with gap-filling, but only partly. The ensemble produced values in 1.5–1.6 times as many 30-min windows as the best single pair, yet every GF requires EC turbulence measurements (u* for the AE and WP, and the tracer flux for the MBR) and so is unavailable during most EC outages. GF is therefore better suited to filling periods when only the target-gas analyzer of the EC system fails, and to extending EC towers to additional gases measured by slow-response analyzers on the same profile, than to replacing EC during complete outages. Third, profile analyzers already deployed for storage-flux calculations, as at NEON, provide this capability at little extra cost.”
New Appendix C: “Choosing sampling-height pairs without an eddy-covariance reference”, including C.5 “When another gas has EC: use the pairs reliable for CO₂” and a recommended protocol (C.6).
Citation: https://doi.org/10.5194/essd-2026-413-AC2
Data sets
Benchmark dataset for comparing eddy covariance measurements and gradient flux calculations across 47 National Ecological Observatory Network (NEON) flux towers, 2021–2024 Sparkle L. Malone et al. https://portal.edirepository.org/nis/mapbrowse?scope=edi&identifier=2369&revision=1&key=_LHRQgDvJXjj7acnumjrl-j3mEQ
Model code and software
lterwg-flux-gradient Sparkle L. Malone et al. https://github.com/lter/lterwg-flux-gradient
Interactive computing environment
lterwg-flux-gradient-eval Sparkle L. Malone et al. https://github.com/lter/lterwg-flux-gradient-eval
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 434 | 232 | 173 | 839 | 219 | 178 |
- HTML: 434
- PDF: 232
- XML: 173
- Total: 839
- BibTeX: 219
- EndNote: 178
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
This manuscript uses 47 NEON sites to compare the eddy covariance and gradient flux methods for estimating CO2 and H2O fluxes. This analysis is valuable for considering the conditions under which the two methods agree and disagree and is unique in the sense that NEON infrastructure and data makes a systematic comparison possible. A better understanding of how well gradient fluxes capture ecosystem gas exchange is particularly valuable for trace gas fluxes which are more difficult or impossible to measure with eddy covariance. Overall the analysis is extremely thorough and contains multitudes of information about how seasonality, and ecosystem structure affect agreement between the two methods, as well as technically useful information such as measurement height for gradient fluxes. One conclusion that I found particularly interesting is the ensemble gradient flux approach in which multiple gradient flux methods, each with different assumptions or limitations, can be combined to produce a much more temporally complete and informative dataset. Rather than selecting a single method, the ensemble approach is intriguing and, in my opinion, a good way to approach gappy data.
The analyses and figures in this manuscript are excellent, and I am impressed with the figures that synthesise huge amounts of information and add conceptual or diagrammatic elements that help with a more intuitive understanding of the information. I did struggle with the many abbreviations and acronyms for the flux methods, statistical methods, sensor placement, particularly when a sentence had multiple acronyms. I am not sure the best way around this, but two ideas are to either drop some of the acronyms and write out the term, or include a glossary table. My recommendation is to publish this manuscript with minor revisions.
Minor and additional comments and suggestions are listed below.
Line 100: ‘With the AE, K is derived from wind and turbulence statistics’ – this was one early sentence I struggled with
Line 129-130: I don’t fully understand this. I assume that for gradient flux the gas concentrations are measured at multiple heights (between 4-8 sampling levels) while eddy covariance is above the canopy. I am therefore confused which sensors are both above the canopy (AA) and which sensor would be within the canopy for (AW). Since gradient has multiple sensors, while eddy covariance is one sensor I assume that the eddy covariance measure is always above (A), but I don’t understand what the other paired measurement is. In figure 1B which category do JORN, KONZ, GUAN, and HARV fall into? Are JORN and KONZ both AA because there is not enough canopy to be ‘within’ while GUAN and HARV are AW because some instruments are above and some within?
Section 2.4.1
In this dataset are the target gas and tracer different combinations of CO and H2O? For example to calculate CO2 flux with MBR, the H2O is the tracer, and to calculate H2O flux, CO2 is the tracer?
Figure 2: This figure is helpful! A very minor aesthetic comment, but I suggest aligning the edge of the Q3 box with the Q2 box to show that they both emerge from Q1. I presume the idea is that Q2 and Q3 are based on the outcomes of Q1?
Line 357: I don’t understand how the ability of GF to capture different ecosystem components was evaluated. Does the footprint of EC and GF differ? Why would GF capture different ecosystem components from EC?
Figure 3: It would be helpful to remind readers in the caption that CCC varies from 1 to -1 which represent strong concordance, while values close to 0 represent weak concordance.
Line 408: very minor but the font size changed in this section
Figure 6F: it’s interesting that TOOL is the only site with low canopy and CCC<0.5. Do the authors have any thoughts on why this site has lower concordance, given the low canopy?
Figure 8: the figure appears to be missing from the main text, I can only see the caption.
Lines 515-525: I think I understand the AA and AW terminology better from this section. How I understand it now is that AA and AW is not refering to EC vs GF, it’s about the paired sample height comparisons for a GF estimate. AA and AW is referring to pairs of sampling heights for gradient flux calculations where the paired gradient measurements are either both above or above and within the canopy.
Line 554-555: It’s unclear to me how GF in itself helps with decomposition into source contributions? Is this because below canopy gradient heights can represent primarily respiration/evaporation contributions, while above-canopy would include both respiration and photosynthesis/transpiration? I am not sure that the analysis in this paper particularly helps highlight that possibility with GF and found the statements about partitioning here, and above (line 357), more confusing than helpful.
Figure A11 and Figure 6 look almost like the same figure to me, though there seems to be some slight difference in the underlying data as panels C, D, E, F are slightly different.
Figure A13: some of the ranges of flux values seem very large eg: CPER with a values at -1000, but even other sites with axis values at +/- 200.
Figure A15: I was curious to see the relationship for TOOL but it’s missing entirely from this graph and A13. It does appear in the H2O fluxes so perhaps all CO2 flux measurements were lost to filtering.