the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 29 Jun 2026
Beyond the simple mean: an alternative way to improve multi-model bottom-up wetland CH4 estimates
Abstract. Wetlands are the largest natural source of atmospheric methane, yet bottom-up estimates of their emissions remain highly uncertain due to structural differences, parametric uncertainties among process-based models and strong environmental heterogeneity. In addition, most wetland emission models are not systematically calibrated against methane flux measurements, limiting their ability to capture realistic spatial and temporal dynamics. The increasing availability of eddy-covariance towers measuring CH4 fluxes across diverse wetland types now offers new opportunities to evaluate and constrain model ensembles using observational data. Multi-model ensembles are commonly used to quantify uncertainty, but the widespread use of simple model averaging implicitly assumes that all models contribute equally and optimally across sites, an assumption that is rarely justified. Here, we present a data-driven framework that moves beyond the simple mean to derive site-adaptive ensemble estimates and to characterize spatial patterns of model disagreement, and how they are reshaped by ensemble weighting, in wetland CH4 emissions.
Using flux observations from a global network of 44 wetland sites, ranging from boreal arctic to tropical wetland ecosystems, and simulations from sixteen global wetland biogeochemistry models, we estimate site-specific optimal ensemble weights via a Bayesian model averaging framework fitted by Expectation–Maximization (EM-BMA). To improve robustness, weights are stabilized through resampling, and sites are clustered based on their stabilized weight signatures in compositional space, yielding groups of locations with similar model skill structures. Within each cluster, we perform cross-validated predictions and compare EM-BMA against simple model averaging (SMA) using standard performance metrics (R2, normalized RMSE, and mean bias).
To interpret these clusters and to delineate the conditions under which different model combinations are best fitting the site measurements, we relate cluster membership to environmental predictors. We characterize cluster-specific validity domains in predictor space using low-dimensional projections and geometric and probabilistic envelopes, and we identify the most influential predictors using a machine-learning classifier with back-projected feature importance.
We show that optimal model combinations vary across sites and that Bayesian model averaging outperforms simple model averaging in cross-validation. When propagated beyond the measurement sites using environmental predictors and applied to model diagnostic outputs, the resulting global wetland CH4 emission estimates differ only slightly from those obtained with simple averaging (less than 5 %). However, substantial differences emerge at local and regional scales, highlighting the importance of accounting for spatial heterogeneity in model skill. This framework therefore provides a transparent and reproducible alternative to equal-weight ensemble means for improving bottom-up wetland CH4 estimates across heterogeneous environments.
- Preprint
(7581 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
- RC1: 'Comment on egusphere-2026-2604', Anonymous Referee #1, 02 Sep 2026
-
RC2: 'Comment on egusphere-2026-2604', Anonymous Referee #2, 21 Oct 2026
This is a timely, well-motivated contribution that I expect will be publishable after revision. Besic et al. replace the equal-weight ensemble mean of 16 global wetland CH4 models with site-specific weights fitted by Expectation–Maximization Bayesian model averaging (EM-BMA) against monthly fluxes at 44 sites. They stabilize the weights by resampling, cluster sites in CLR-transformed weight space (four clusters), relate clusters to environmental predictors with a PCA + random forest classifier (accuracy ~0.70), and map cluster-specific mixtures globally. The major results are clear and useful: BMA outperforms simple model averaging (SMA) in cross-validation (R² 0.62 vs 0.06), global totals change by less than 5%, but regional emissions shift by tens of percent, and between-model spread is redistributed rather than reduced. The candid framing as a proof of concept is appropriate.
My main concerns are about whether the evaluation design supports the strength of the BMA-vs-SMA claim, the scale and units of the model–tower comparison, and how far the site-derived clusters can be trusted when extrapolated to under-sampled wetland types. None of these undermine the idea; addressing them would make the framework considerably more convincing and reusable.
Major Strengths
- This work addresses an important, under-addressed problem. Equal weighting of wetland CH4 models is the default in global budget work despite well-documented structural differences.
- Leveraging FLUXNET-CH4 and BAWLD to constrain an ensemble, rather than single models, is the right direction and makes a concrete case for the value of tower networks.
- Treating weights as compositional data, stabilizing weights by resampling, and reporting weight uncertainty are all sound.
- Separating between-model from within-model variance moves beyond min–max envelopes, and the finding that BMA redistributes rather than reduces structural spread is honest and informative.
- That global totals barely change while regional patterns shift substantially is an important message for budget users who might otherwise assume weighting only matters globally.
- Code and data are on Zenodo, and the workflow is easy to follow.
Major comments
- The BMA-vs-SMA comparison is not like-for-like, and the cross-validation design needs clarification.
BMA fits up to 32 parameters per site (16 weights, 16 variances) to the observations, while SMA fits none. A large gain over SMA (R² 0.62 vs 0.06; MBE 1.38 vs 23.11) is therefore expected and partly reflects calibration rather than a better combination rule. I suggest:
- Adding stronger baselines: best single model per cluster, skill-weighted means (e.g., inverse-RMSE), a single global BMA weight vector, and a bias-corrected SMA. This would show how much of the gain comes from site-adaptivity specifically.
- Stating how the K folds are formed (Sect. 3.1.4). If folds are random site-months, temporal autocorrelation and same-site leakage will inflate skill. Leave-one-site-out (within cluster) is the test that matters for upscaling.
- Evaluating the full upscaling chain end to end: hold out a site, assign its cluster via the classifier from predictors only, apply that cluster's weights, and compare against SMA. With classifier accuracy ~0.70, this is the realistic skill of the global product.
- Revisiting Sect. 5.5. Clusters are defined from weights fitted to all data, including later validation months, so the dependence is not obviously "limited". A nested scheme (clustering re-derived within each training fold) would resolve this.
- Lines 125–127 state that model outputs already include wetland fractional area, which implies grid-cell-mean fluxes at 0.5°. Towers measure flux per unit wetland area over footprints of hundreds of metres. If grid-cell fluxes were compared directly with tower fluxes, the weights partly reward models with larger simulated wetland fraction in the site pixel rather than better flux physics, and prognostic and prescribed-extent models are affected differently. Please state exactly which model variable was extracted (flux density per wetland area vs per grid area), how sites were matched to cells, and how prognostic vs prescribed extent affects the weights. A short sensitivity test using per-wetland-area flux would be very helpful.
Relatedly, units should be made consistent: the text gives mg m⁻² s⁻¹, figures use µg m⁻² s⁻¹, the rescaling factor is not reported, and the MBE values in Fig. 3c (1.38, 23.11) appear to be in rescaled units that do not match the plotted axes (0–1.6).
- Monthly CH4 fluxes are non-negative, right-skewed and heteroscedastic, while Eq. (1) uses Gaussian kernels with constant σk and no bias-correction terms (unlike Raftery et al., 2005). With 16 strongly co-varying seasonal model outputs, EM mixture weights can be poorly identified, and the sparsity seen in Fig. 2 may partly reflect this rather than true skill differences. Please (a) clarify whether any additive or multiplicative bias correction is applied, (b) consider a log or other variance-stabilizing transform as a sensitivity test, (c) report the model-to-model correlation structure at sites, and (d) say how many resamples were used and whether the bootstrap preserve’s temporal structure (block or seasonal resampling). Note also that EM yields maximum-likelihood weights; "posterior" weights should be described accordingly.
The near-zero within-model variance for Clusters 1 and 2 in Fig. 4a is surprising and may indicate σk collapsing to the regularization floor. Please report the floor values and how often they bind, and consider normalizing variances by observed variance so clusters with different flux magnitudes can be compared.
- Removing zero-flux site-months (Lines 112–116) focuses the evaluation on emitting periods, but the weights are then applied to annual totals that include cold-season and shoulder-season emissions, which matter for boreal budgets. Please clarify whether negative fluxes were also removed, how many site-months were dropped per site and season, and discuss how this choice could bias weights toward models that perform well in summer.
- With 26 boreal, 14 temperate and 4 tropical sites, the clusters assigned across tropical wetlands (Fig. 7) rest on very few observations, yet the tropics contribute most of the global flux and most of the BMA–SMA difference (Fig. 8c–d; e.g., the Sudd). I recommend:
- Using the validity domains from Sect. 3.2.2 operationally: mask or hatch pixels that fall outside all cluster envelopes in Figs. 7, 8 and C1, and report the fraction of global wetland area and emissions affected.
- Given that Fig. B1 shows assignment probabilities of 0.3–0.6 over much of the domain, replacing hard cluster assignment with a probability-weighted mixture of cluster weights, or at least showing how the global and regional results change.
- Discussing wetland types absent from the training set. Coastal and saline wetlands (mangroves, tidal marshes), where sulfate suppresses methanogenesis, are not represented by the predictors (no salinity, tidal or hydrogeomorphic variable) yet are included in the WAD2M mask in places. The same applies to rice-adjacent and floodplain systems.
- The classifier and importance analysis deserve more scrutiny given n = 44. Gravel content is the top predictor (Fig. 6b), which is hard to interpret physically for wetlands and may act as a site-identity proxy. Please report variability of importances across CV folds or repeated fits, justify the PCA step before the random forest with only nine predictors, explain how importances are back-projected through PCA loadings, and clarify "k=5" (Line 315). The text in Sect. 4.7 cites bulk density, which does not appear in Fig. 6b. Predictors relevant to CH4 that are missing include soil organic carbon or peat depth, inundation dynamics and wetland type; HWSD is weak in organic soils. The predictor period (2003–2022) and climate source (ERA5) also differ from the model period (2000–2020) and forcing (CRU, GSWP3); a sentence on why this does not matter would help.
- Near-uniform weights can arise because all models perform comparably well, or because all perform comparably poorly and EM cannot discriminate. Cluster 3 has the largest between-model variance relative to within-model variance (Fig. 4a), which better supports the latter. Please report per-model skill within Cluster 3 and temper the conclusion at Lines 450–451 accordingly. Since Cluster 3 dominates boreal and tropical wetlands in Fig. 7, this interpretation matters for the paper's main message.
- Fig. 8c shows BMA consistently ~5–8 Tg CH4 yr⁻¹ above SMA (roughly 4–5%), driven by the temperate and tropical bands. Please report these numbers explicitly by latitude band and forcing, and discuss whether BMA moves the bottom-up estimate toward or away from top-down inversions (Saunois et al., 2025). This comparison would substantially raise the paper's relevance. Please also explain how the between-model variance envelopes in Fig. 8c–d are aggregated across pixels with different cluster weights, and why the time series end in 2016 when simulations cover 2000–2020.
Minor comments
- Please explain how monthly series were obtained for the 8 BAWLD sites, and clarify the role of ABCFlux v2 (Virkkala et al., 2025), which is cited at Line 86 but not described. A supplementary table of sites (ID, coordinates, wetland type, years, number of months used, cluster) would be very useful.
- Cluster number. Show the elbow curve (supplement) and report cluster sizes (from Fig. 6a they appear to be 9, 12, 6 and 17).
- Fig. 4b: Specify whether the PCA is on raw or CLR-transformed weights; the overlapping labels are hard to read.
- 2. The observed axes in the right-hand scatter panels run in decreasing order, which is unusual and confusing. Please also state the time aggregation shown in the middle panels.
- Fig. 3c: This is pooled cross-validated output; please provide units for MBE.
- Fig. 3b:. Absolute standard deviations are hard to judge when centroid weights differ by an order of magnitude; a coefficient of variation or a rank-stability metric would be more informative.
- Fig. 8: Add a small table of regional BMA vs SMA means and percentage differences for both forcings, so readers do not have to read values off the panels.
Citation: https://doi.org/10.5194/egusphere-2026-2604-RC2
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 302 | 152 | 51 | 505 | 56 | 51 |
- HTML: 302
- PDF: 152
- XML: 51
- Total: 505
- BibTeX: 56
- EndNote: 51
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
This manuscript presents a site-adaptive ensemble framework for estimating wetland methane emissions from 16 process-based models and observations at 44 wetland sites. The study addresses an important problem. Equal weighting remains common in multi-model wetland CH₄ assessments despite substantial spatial variation in model performance. The attempt to combine site observations, multi-model simulations, compositional clustering, environmental classification, and spatial extrapolation is potentially valuable. The manuscript is generally well organized, the workflow is transparent, and the recognition that model skill may vary among environmental regimes is scientifically reasonable.
However, several methodological issues currently limit confidence in the reported predictive improvement and especially in the global extrapolation. The most important concerns are the potential dependence between cluster construction and cross-validation, the very limited observational basis for global classification, the treatment of zero and non-emitting observations, the mismatch between site measurements and coarse-grid model outputs, and insufficient quantification of uncertainty in spatial cluster assignment. These issues require substantial additional analysis and clarification. I therefore recommend major revision.
Major comments
Line 99: Add a comma after “Table 1”: “These models, presented in Table 1, differ…”
Lines 151–153: GPP does not itself provide “proxies of methane emissions.” Suggested wording: “GPP was included as a proxy for vegetation productivity and potential carbon-substrate supply…”
ine 153: “vegetation activities” should be “vegetation activity.”
Line 230: Replace “As evoked in Sec. 2.3” with “As described in Sect. 2.3.”
Line 296: “This confirmes” should be “This confirms.”
Lines 318–320: Replace “This confirms” with “This suggests.”
Lines 455–457: “allow the propagation” could be revised to “enable the propagation.”
Lines 457–458: “Bayesian model averaged emissions” should be hyphenated as “Bayesian-model-averaged emissions,” or simply “BMA emissions.”