Multi-source land-use emissions reveal rising airborne fraction

Background

Whether the airborne fraction (AF) is increasing or approximately constant remains contested^1,2,8–12. In its classical form, AF is a yearly ratio of atmospheric growth to total anthropogenic emissions, computed as the sum of fossil fuel emissions and land-use and land-cover change emissions:

AF_t = \frac{G_t}{FF_t + LULC_t}, \tag{1}

where G_t is the annual atmospheric CO_2 growth, FF_t is fossil fuel emissions excluding carbonation, and LULC_t is land-use and land-cover change emissions. AF is a key carbon-cycle diagnostic, with implications for carbon-cycle feedbacks and near-term mitigation planning^1–3.

A persistent concern is that AF inference depends on the treatment of land-use and land-cover change (LULC) emissions, which are uncertain and model-dependent in annual carbon-budget accounting³. The Global Carbon Budget (GCB) 2025 dataset provides one column of LULC emissions as the average of three bookkeeping models (BLUE, OSCAR, LUCE), but a broader set of model-based LULC alternatives can be constructed from the same source data. If the information from this broader set of LULC measurements is not incorporated into AF trend inference, tests can be underpowered and inference on trend direction becomes less reliable.

Here we address that issue with a design that incorporates measurement uncertainty to obtain more reliable trend estimates^4–6. We present a two-stage estimator that propagates denominator uncertainty from repeated LULC measurements into annual AF variance and then estimates the AF trend by weighted least squares (WLS) with heteroskedasticity- and autocorrelation-consistent (HAC) inference⁷. The approach uses cross-measurement dispersion to weight years by precision, rather than treating all years as equally precise as in ordinary least squares (OLS). We find that WLS delivers statistically significant and more stable evidence of positive trend, including under exclusion of the final observation (2024), which shows a large AF jump. For comparison, we estimate OLS on the same AF series, with conventional and HAC standard errors showing weaker evidence of a positive trend and greater sensitivity to endpoint exclusion. This clarifies why evidence based on OLS does not support a clear conclusion about AF trends.

Data

We use annual Global Carbon Budget 2025 data for 1959-2024, with atmospheric growth G_t from NOAA/ESRL global concentration trends¹³, fossil emissions excluding carbonation FF_t from the Global Carbon Project fossil dataset³, and a panel of 69 LULC measurements per year: BLUE¹⁴, OSCAR¹⁵, LUCE¹⁶, and peat-augmented process-based land-model combinations drawn from the GCB model ensemble^17–37, combined with peat components^38–40 (see Methods). The data is shown in Figure 1.

Two LULC means are constructed from the panel. First, the GCB LULC mean, defined as the mean of the three bookkeeping models BLUE, OSCAR, and LUCE, corresponds to the LULC column used in the Global Carbon Budget. Second, the all‑LULC mean, defined as the cross-series mean from all 69 LULC measurements. Uncertainty from the multiple LULC measurements are estimated by the cross-measurement dispersion from the full panel.

Figure 1: Global Carbon Budget 2025 data for fossil emissions, atmospheric growth, and LULC emissions (1959-2024). For the LULC panel, we show all extracted and derived LULC measurements together with the GCB LULC column and the all-LULC cross-series mean.

Identification Strategy

The empirical question is whether AF has a positive linear trend. The key design choice is how to handle annual denominator uncertainty from multiple LULC measurements. OLS assigns equal precision to all years and therefore discards information from cross-measurement dispersion. Our preferred specification uses delta-method variance estimates^41,42 to construct year-specific weights and then estimates the trend by WLS, with HAC standard errors for inference^7,43 (see Methods).

Results

Full-sample trends

Table 1 reports full-sample trend estimates for the four specifications: OLS using GCB LULC, OLS using the all-LULC mean, WLS using the GCB LULC with all-LULC weights, and WLS using the all-LULC mean with all-LULC weights.

Table 1: Full-sample trend comparison across OLS and WLS specifications.

Trend, full sample	OLS (GCB)	OLS (all)	WLS (GCB)	WLS (all)
Estimate	0.001581	0.001172	0.002649	0.002305
Standard error	0.000771	0.000792	0.000022	0.000023
HAC standard error	0.000621	0.000643	0.000657	0.000672
p-value	0.040212	0.138912	0.000000	0.000000
HAC p-value	0.010837	0.068425	0.000056	0.000609
R-squared	0.061706	0.033086	0.148987	0.111070

The main finding is that incorporating denominator-measurement information via WLS substantially changes inference relative to both OLS benchmarks. In the full sample, all approaches produce a positive slope, and both WLS specifications provide much stronger evidence than OLS (HAC p-values 0.000056 and 0.000609 versus 0.010837 and 0.068425). The estimated WLS trend is also steeper (0.002305 per year versus 0.001581 and 0.001172), consistent with greater weight on years with lower denominator uncertainty. The WLS slope implies an increase of about +0.15 in AF over 1959-2024, compared with about +0.10 and +0.08 from the two OLS specifications. The R-squared is also higher for WLS, indicating a better fit when accounting for measurement uncertainty. Nevertheless, the R-squared values are small in absolute terms, signaling that AF is a noisy variable and that the trend is only one component of its variation.

Figure 2 shows the AF series using the mean of all LULC measurements together with the associated variance, the OLS trend, and the delta-method WLS trend when using the full sample. The gray bars show the relative weights (standardised so that the maximum weight is 1) assigned to each year in the WLS estimation, which are inversely proportional to the estimated variance of the AF estimate for that year.

Figure 2: AF series with OLS and delta-method WLS trend (full sample, 1959-2024). Grey bars show the relative weights assigned to each year in the WLS estimation.

Note that the WLS trends are steeper than the OLS trends, which is consistent with the numerical results in Table 1. The WLS specifications assign more weight to the more precisely measured years (e.g., 2000-2020 as shown by the width of the shaded area) and less weight to noisier years (e.g., 1960s-1970s), which helps reveal the positive trend signal in the data.

Endpoint Robustness

The final observation (2024) shows a large increase in AF. To ensure this point is not mechanically driving the result, we re-estimate the same four specifications on the subsample ending in 2023. Results are shown in Table 2 and Figure 3.

Table 2: Subsample (ending in 2023) trend comparison across OLS and WLS specifications.

Trend (up to 2023)	OLS (GCB)	OLS (all)	WLS (GCB)	WLS (all)
Estimate	0.001296	0.000878	0.002367	0.002015
Standard error	0.000777	0.000798	0.000022	0.000023
HAC standard error	0.000578	0.000600	0.000600	0.000615
p-value	0.095163	0.271089	0.000000	0.000000
HAC p-value	0.024954	0.143217	0.000081	0.001059
R-squared	0.042332	0.018863	0.125521	0.089297

Results remain qualitatively unchanged for the WLS specifications: both estimated slopes stay positive and statistically significant (0.002367 with HAC p-value = 0.000081, and 0.002015 with HAC p-value = 0.001059). For the sample ending in 2023, OLS using GCB LULC is positive but only marginal under conventional inference (p-value = 0.095), while OLS using the all-LULC mean is not significant (p-value = 0.271). Overall, the WLS results are robust to endpoint exclusion, while OLS results are more sensitive to denominator construction and provide weaker evidence of a positive trend.

Figure 3: AF series with OLS and delta-method WLS trend (sample ending in 2023). Grey bars show the relative weights assigned to each year in the WLS estimation.

Supplementary Numerical Results

While the main focus is on the slope estimates, we also report intercept estimates for completeness. The intercept captures the baseline AF level in 1959, and its estimation can also be affected by denominator measurement error.

Results for the intercept estimates are shown in Table 3.

Table 3: Intercept comparison across OLS and WLS specifications for full sample and sample up to 2023.

Intercept	OLS (GCB) full	OLS (all) full	WLS (GCB) full	WLS (all) full	OLS (GCB) 2023	OLS (all) 2023	WLS (GCB) 2023	WLS (all) 2023
Estimate	-2.7266	-1.8946	-4.9141	-4.2130	-2.1611	-1.3121	-4.3539	-3.6377
Standard error	1.5352	1.5777	0.0443	0.0456	1.5461	1.5888	0.0440	0.0453
HAC standard error	1.2391	1.2845	1.3212	1.3518	1.1553	1.1994	1.2087	1.2390
p-value	0.0757	0.2298	0.0000	0.0000	0.1622	0.4089	0.0000	0.0000
HAC p-value	0.0278	0.1402	0.0002	0.0018	0.0614	0.2739	0.0003	0.0033
R-squared	0.0617	0.0331	0.1490	0.1111	0.0423	0.0189	0.1255	0.0893

Discussion and conclusion

Using a measurement-error-aware framework, we find robust evidence that AF increased from 1959 to 2024. The key methodological point is that incorporating multi-source denominator uncertainty through delta-method WLS materially changes inference relative to plain OLS. The endpoint test shows that this conclusion is not driven by the 2024 jump.

In practice, the implication is clear: AF trend assessments should report measurement-error-aware weighted estimates, benchmark them against OLS, and document endpoint sensitivity.

A rising AF implies that a larger share of emitted CO_2 remains in the atmosphere over policy-relevant horizons, tightening near-term mitigation requirements for a given temperature objective^1–3. The evidence of an increasing AF is consistent with broader findings that the climate system is out of energy balance and has recently exhibited elevated warming and heating rates; these results should therefore be interpreted within that wider risk context^44–47.

References

1.

Canadell, J. G. et al. Contributions to accelerating atmospheric CO_2 growth from economic activity, carbon intensity, and efficiency of natural sinks. Proceedings of the National Academy of Sciences 104, 18866–18870 (2007).

2.

Raupach, M. R. et al. Global and regional drivers of accelerating CO_2 emissions. Proceedings of the National Academy of Sciences 104, 10288–10293 (2007).

3.

Friedlingstein, P. et al. Global carbon budget 2025. Earth System Science Data 18, 3211–3288 (2026).

4.

Fuller, W. A. Measurement Error Models. (Wiley, New York, 1987).

5.

Carroll, R. J., Ruppert, D., Stefanski, L. A. & Crainiceanu, C. M. Measurement Error in Nonlinear Models: A Modern Perspective. (Chapman & Hall/CRC, Boca Raton, 2006).

6.

Vera-Valdés, J. E. & Grivas, C. Robust estimation of carbon dioxide airborne fraction under measurement errors. Environmental Research Communications 7, 031009 (2025).

7.

Newey, W. K. & West, K. D. A simple, positive semi-definite, heteroskedasticity and autocorrelation consistent covariance matrix. Econometrica 55, 703–708 (1987).

8.

Knorr, W. Is the airborne fraction of anthropogenic CO_2 emissions increasing? Geophysical Research Letters 36, L21710 (2009).

9.

Ballantyne, A. P., Alden, C. B., Miller, J. B., Tans, P. P. & White, J. W. C. Increase in observed net carbon dioxide uptake by land and oceans during the past 50 years. Nature 488, 70–72 (2012).

10.

Le Quéré, C. et al. Trends in the sources and sinks of carbon dioxide. Nature Geoscience 2, 831–836 (2009).

11.

Bennedsen, M., Hillebrand, E. & Koopman, S. J. On the evidence of a trend in the CO2 airborne fraction. Nature 616, E1–E3 (2023).

12.

Bennedsen, M., Hillebrand, E. & Koopman, S. J. A regression-based approach to the CO2 airborne fraction. Nature Communications 15, 8507 (2024).

13.

Lan, X., Tans, P. & Thoning, K. W. Trends in globally-averaged CO2 determined from NOAA global monitoring laboratory measurements. (2026) doi:10.15138/9N0H-ZH07.

14.

Hansis, E., Davis, S. J. & Pongratz, J. Relevance of methodological choices for accounting of land use change carbon fluxes. Global Biogeochemical Cycles 29, 1230–1246 (2015).

15.

Gasser, T. et al. Historical CO2 emissions from land use and land cover change and their uncertainty. Biogeosciences 17, 4075–4101 (2020).

16.

Qin, Z. et al. Global spatially explicit carbon emissions from land-use change over the past six decades (1961-2020). One Earth 7, 835–847 (2024).

17.

Haverd, V. et al. A new version of the CABLE land surface model (subversion revision r4601) incorporating land use and land cover change, woody vegetation demography, and a novel optimisation-based approach to plant coordination of photosynthesis. Geoscientific Model Development 11, 2995–3026 (2018).

18.

Melton, J. R. et al. CLASSIC v1.0: The open-source community successor to the canadian land surface scheme (CLASS) and the canadian terrestrial ecosystem model (CTEM) – part 1: Model framework and site-level performance. Geoscientific Model Development 13, 2825–2850 (2020).

19.

Lawrence, D. M. et al. The community land model version 5: Description of new features, benchmarking, and impact of forcing uncertainty. Journal of Advances in Modeling Earth Systems 11, 4245–4287 (2019).

20.

Fisher, R. A. et al. Taking off the training wheels: The properties of a dynamic vegetation model without climate envelopes, CLM4.5(ED). Geoscientific Model Development 8, 3593–3619 (2015).

21.

Tian, H. et al. North american terrestrial CO2 uptake largely offset by CH4 and N2O emissions: Toward a full accounting of the greenhouse gas budget. Climatic Change 129, 423–426 (2015).

22.

Ma, L. et al. Global evaluation of the ecosystem demography model (ED v3.0). Geoscientific Model Development 15, 1971–1994 (2022).

23.

Yang, X., Thornton, P., Ricciuto, D., Wang, Y. & Hoffman, F. Global evaluation of terrestrial biogeochemistry in the energy exascale earth system model (E3SM) and the role of the phosphorus cycle in the historical terrestrial carbon balance. Biogeosciences 20, 2813–2836 (2023).

24.

Needham, J. et al. Vertical canopy gradients of respiration drive plant carbon budgets and leaf area index. New Phytologist 246, 144–157 (2025).

25.

Felzer, B. S. & Jiang, M. Effect of land use and land cover change in context of growth enhancements in the united states since 1700: Net source or sink? Journal of Geophysical Research: Biogeosciences 123, 3439–3457 (2018).

26.

Xia, J. Z. et al. The carbon budget of china: 1980–2021. Science Bulletin 69, 114–124 (2024).

27.

Yue, X. et al. Development and evaluation of the interactive model for air pollution and land ecosystems (iMAPLE) version 1.0. Geoscientific Model Development 17, 4621–4642 (2024).

28.

Shu, S., Jain, A. K., Koven, C. D. & Mishra, U. Estimation of permafrost SOC stock and turnover time using a land surface model with vertical heterogeneity of permafrost soils. Global Biogeochemical Cycles 34, e2020GB006585 (2020).

29.

Reick, C. H. et al. JSBACH 3 - the land component of the MPI earth system model: Documentation of version 3.2. (2021) doi:10.17617/2.3279802.

30.

Poulter, B., Frank, D. C., Hodson, E. L. & Zimmermann, N. E. Impacts of land cover and climate data selection on understanding terrestrial carbon dynamics and the CO2 airborne fraction. Biogeosciences 8, 2027–2036 (2011).

31.

Smith, B. et al. Implications of incorporating n cycling and n limitations on primary production in an individual-based dynamic vegetation model. Biogeosciences 11, 2027–2054 (2014).

32.

Schaphoff, S. et al. LPJmL4 – a dynamic global vegetation model with managed land – part 1: Model description. Geoscientific Model Development 11, 1343–1375 (2018).

33.

Lienert, S. & Joos, F. A bayesian ensemble data assimilation to constrain model parameters and land-use carbon emissions. Biogeosciences 15, 2909–2930 (2018).

34.

Vuichard, N. et al. Accounting for carbon and nitrogen interactions in the global terrestrial ecosystem model ORCHIDEE (trunk version, rev 4999): Multi-scale evaluation of gross primary production. Geoscientific Model Development 12, 4751–4779 (2019).

35.

Walker, A. P. et al. The impact of alternative trait-scaling hypotheses for the maximum photosynthetic carboxylation rate (v-cmax) on global gross primary production. New Phytologist 215, 1370–1386 (2017).

36.

Kato, E., Kinoshita, T., Ito, A., Kawamiya, M. & Yamagata, Y. Evaluation of spatially explicit emission scenario of land-use change and biomass burning using a process-based biogeochemical model. Journal of Land Use Science 8, 104–122 (2013).

37.

Ito, A. Disequilibrium of terrestrial ecosystem CO2 budget caused by disturbance-induced emissions and non-CO2 carbon export flows: A global model assessment. Earth System Dynamics 10, 685–709 (2019).

38.

Conchedda, G. & Tubiello, F. N. Drainage of organic soils and GHG emissions: Validation with country data. Earth System Science Data 12, 3113–3137 (2020).

39.

Müller, J. & Joos, F. Committed and projected future changes in global peatlands – continued transient model simulations since the last glacial maximum. Biogeosciences 18, 3657–3687 (2021).

40.

Qiu, C. et al. Large historical carbon emissions from cultivated northern peatlands. Science Advances 7, eabf1332 (2021).

41.

Cochran, W. G. Sampling Techniques. (Wiley, New York, 1977).

42.

Oehlert, G. W. A note on the delta method. The American Statistician 46, 27–29 (1992).

43.

Aitken, A. C. On least squares and linear combination of observations. Proceedings of the Royal Society of Edinburgh 55, 42–48 (1935).

44.

Minière, A., Schuckmann, K. von, Sallée, J.-B. & Vogt, L. Robust acceleration of earth system heating observed over the past six decades. Scientific Reports 13, 22975 (2023).

45.

World Meteorological Organization. State of the global climate 2025. (2026) doi:10.59327/WMO/S/CRI/SOC/1.

46.

Storto, A. & Yang, C. Acceleration of the ocean warming from 1961 to 2022 unveiled by large-ensemble reanalyses. Nature Communications 15, 545 (2024).

47.

Foster, G. & Rahmstorf, S. Global warming has accelerated significantly. Geophysical Research Letters 53, e2025GL118804 (2026).

Methods

Our primary estimator is a two-stage measurement-error trend model using repeated yearly denominator measurements, followed by weighted trend estimation with HAC inference^4,5,7,42,43.

Construction of the LULC measurement panel

The repeated LULC measurements are built from the Global Carbon Budget 2025 dataset. We first extract the three bookkeeping series (BLUE, OSCAR, LUCE)^14–16. Then, for each process-based land-model LULC model in the GCB ensemble which does not already include peat emissions^17–37, we add the corresponding peat component to make it comparable to the bookkeeping models, which include peat emissions by construction. Specifically, for each of the 33 process-based land-model combinations in the GCB ensemble we create three peat-augmented variants by adding FAO_peat, LPX_Bern_peat, and ORCHIDEE_peat^38–40. This produces 66 derived series, and together with BLUE/OSCAR/LUCE gives a panel of 69 yearly LULC measurements.

In the estimation, the yearly denominator is computed as \hat C_t = FF_t + \bar{LULC}_t, \qquad \bar{LULC}_t = \frac{1}{n_t}\sum_{j=1}^{n_t} LULC_{tj},

where n_t=3 for the GCB LULC mean and n_t=69 for the all-LULC mean. The variance of \hat C_t is estimated from the cross-measurement dispersion across the 69 LULC values and is used to construct the year-specific AF variance estimate used for WLS weighting.

Assume for each time t=1,\dots,T we observe:

• a single numerator b_t (e.g., atmospheric CO_2 growth), and

• multiple denominator measurements c_{t1},\dots,c_{t n_t} with n_t \ge 2, which are noisy observations of a latent C_t.

Using the repeated denominator measurements, we can estimate the variance of the ratio estimator a_t = b_t/C_t via the delta method, which accounts for the variability in C_t due to measurement error. The procedure is described in detail next.

Two-step estimation approach

We use a two-step approach to estimate the trend in the ratio a_t = b_t/C_t over time, accounting for denominator measurement error.

Step 1: Estimate C_t and a_t

1. Estimate the “true” denominator at time t as the sample mean across available LULC measurements: \hat{C}_t = \frac{1}{n_t} \sum_{j=1}^{n_t} c_{tj}.

We consider two variants of this denominator estimate: one using the GCB LULC column (the mean of BLUE, OSCAR, LUCE) and one using the all-LULC mean (the mean across all 69 LULC measurements).

The variance of \hat{C}_t can be estimated from the multiple measurements.

With n_t \ge 2, an empirical estimator is \widehat{\operatorname{Var}}(\hat C_t)= s_{c,t}^2, \qquad s_{c,t}^2=\frac{1}{n_t-1}\sum_{j=1}^{n_t}(c_{tj}-\bar c_t)^2, \qquad \bar c_t=\hat C_t.

2. Form the ratio estimate \hat a_t = b_t/\hat C_t.

3. Use the multivariate delta method (shown below) to get an approximate variance of \hat a_t

\widehat{\operatorname{Var}}(\hat a_t)\approx\left(\frac{b_t}{\hat C_t^2}\right)^2\widehat{\operatorname{Var}}(\hat C_t).

Step 2: WLS regression of \hat{a}_t on time

For each time t, we obtain a point estimate \hat{a}_t and an estimated variance \widehat{\operatorname{Var}}(\hat a_t) that reflects the cross-measurement dispersion. We use these estimates to fit a linear trend via weighted least squares (WLS):

\hat a_t = \alpha + \beta t + \varepsilon_t, \quad \varepsilon_t \sim (0, \sigma_t^2), \quad \sigma_t^2 = \widehat{\operatorname{Var}}(\hat a_t).

In our application, the time index t is the year, and \hat a_t is the estimated AF for that year. The variance \sigma_t^2 captures the standard error variance and the uncertainty in the AF estimate due to denominator measurement error, which varies across years depending on the number and variability of the LULC measurements. We use this variance to weight the regression, giving more weight to years with more precise AF estimates.

Assuming negligible correlation in \varepsilon_t across time, this is WLS with weights w_t=\frac{1}{\sigma_t^2}. If we suspect serial correlation, we can use HAC standard errors for inference on \beta without changing the point estimate. In practice, we use a HAC covariance estimator with a Bartlett kernel and Andrews automatic bandwidth selection⁷. In the results, we report both conventional and HAC standard errors.

Testing for a time trend is then a test of \beta=0 versus \beta\neq 0 in this linear model, and using WLS we have used all the information in the repeated \hat{a}_t’s through the delta‑method approximation.

Delta method for ratio variance estimation

Step-by-step derivation (first-order delta method):

Define the random vector and the function of interest as X_t= \begin{bmatrix} b_t\\ \hat C_t \end{bmatrix}, \qquad g(X_t)=\frac{b_t}{\hat C_t}.

and let \sigma_{b,t}^2 = \operatorname{Var}(b_t),\qquad \sigma_{C,t}^2 = \operatorname{Var}(\hat C_t),\qquad \sigma_{bC,t} = \operatorname{Cov}(b_t,\hat C_t).

1. Linearize g around the mean vector (\mu_{b,t},\mu_{C,t})=\mathbb E[(b_t,\hat C_t)]: g(X_t)\approx g(\mu_{b,t},\mu_{C,t})+\nabla g(\mu_{b,t},\mu_{C,t})^\top (X_t-\mathbb E[X_t]).

2. Compute the gradient: \nabla g(b,C)= \begin{bmatrix} \partial g/\partial b\\ \partial g/\partial C \end{bmatrix} = \begin{bmatrix} 1/C\\ -b/C^2 \end{bmatrix}.

3. Write the covariance matrix of (b_t,\hat C_t): \Sigma_t= \begin{bmatrix} \sigma_{b,t}^2 & \sigma_{bC,t}\\ \sigma_{bC,t} & \sigma_{C,t}^2 \end{bmatrix}.

4. Apply the delta-method variance formula \operatorname{Var}(g(X_t))\approx \nabla g(\mu_{b,t},\mu_{C,t})^\top\Sigma_t\nabla g(\mu_{b,t},\mu_{C,t}), and with plug-in evaluation at (b_t,C_t), this becomes

\operatorname{Var}(\hat a_t)\approx\left(\frac{1}{C_t}\right)^2 \sigma_{b,t}^2+\left(\frac{b_t}{C_t^2}\right)^2 \sigma_{C,t}^2-2\frac{b_t}{C_t^3}\sigma_{bC,t}.

In practice, we use plug-in estimates (replace unknown moments by their empirical counterparts). This is the standard delta‑method approximation for the variance of a ratio estimator.

5. In the case analysed in this paper, we have a single numerator measurement from a robust source per time, so we treat \sigma_{b,t}^2 and \sigma_{bC,t} as negligible relative to the variance from the denominator measurement error, which is captured by \sigma_{C,t}^2.

The expression simplifies to \operatorname{Var}(\hat a_t)\approx\left(\frac{b_t}{C_t^2}\right)^2 \sigma_{C,t}^2.

Accordingly, a practical plug-in estimator is

\widehat{\operatorname{Var}}(\hat a_t)\approx\left(\frac{b_t}{\hat C_t^2}\right)^2\widehat{\operatorname{Var}}(\hat C_t).

Data availability

The study uses publicly available Global Carbon Budget 2025 data and derived yearly series generated from those source files. The processed analysis tables are provided in the repository under the results directory, and the extracted and derived LULC panel is available as a CSV file in the data directory.

Code availability

All code used to process data, estimate models, and generate figures is available in the repository under the scripts directory, including Quarto analysis files and Julia helper functions.

Competing interests

The author declares no competing interests.

Additional information

Supplementary information is not included.