Retrieve raw CMIP data for surface air temperature (tas), incoming shortwave radiation (rsdt), outgoing shortwave radiation (rsut),
outgoing longwave radiation (rlut), and grid cell area (areacella) from the Earth System Grid Federation.
Download data associated with the following experiments: piControl, abrupt-2xCO2, abrupt-4xCO2, 1pctCO2, historical, ssp119, ssp245, ssp370, ssp460, and ssp585.
For a simple implementation of cmip data retrieval, see this repository.
Save the downloaded files to the ./data_raw/CMIP6 directory.
ls ./data_raw/CMIP6 | grep MIROC6areacella_fx_MIROC6_1pctCO2_r1i1p1f1_gn.nc areacella_fx_MIROC6_abrupt-2xCO2_r1i1p1f1_gn.nc areacella_fx_MIROC6_abrupt-4xCO2_r1i1p1f1_gn.nc areacella_fx_MIROC6_hist-nat_r1i1p1f1_gn.nc areacella_fx_MIROC6_historical_r1i1p1f1_gn.nc areacella_fx_MIROC6_historical_r2i1p1f1_gn.nc areacella_fx_MIROC6_historical_r3i1p1f1_gn.nc areacella_fx_MIROC6_historical_r4i1p1f1_gn.nc areacella_fx_MIROC6_piControl_r1i1p1f1_gn.nc areacella_fx_MIROC6_ssp119_r1i1p1f1_gn.nc areacella_fx_MIROC6_ssp245_r1i1p1f1_gn.nc areacella_fx_MIROC6_ssp370_r1i1p1f1_gn.nc areacella_fx_MIROC6_ssp460_r1i1p1f1_gn.nc areacella_fx_MIROC6_ssp585_r1i1p1f1_gn.nc rlut_Amon_MIROC6_1pctCO2_r1i1p1f1_gn_320001-329912.nc rlut_Amon_MIROC6_1pctCO2_r1i1p1f1_gn_330001-334912.nc rlut_Amon_MIROC6_abrupt-2xCO2_r1i1p1f1_gn_320001-329912.nc rlut_Amon_MIROC6_abrupt-2xCO2_r1i1p1f1_gn_330001-339912.nc rlut_Amon_MIROC6_abrupt-2xCO2_r1i1p1f1_gn_340001-344912.nc rlut_Amon_MIROC6_abrupt-4xCO2_r1i1p1f1_gn_320001-329912.nc rlut_Amon_MIROC6_abrupt-4xCO2_r1i1p1f1_gn_330001-334912.nc rlut_Amon_MIROC6_abrupt-4xCO2_r1i1p1f1_gn_335001-344912.nc rlut_Amon_MIROC6_piControl_r1i1p1f1_gn_320001-329912.nc rlut_Amon_MIROC6_piControl_r1i1p1f1_gn_330001-339912.nc rlut_Amon_MIROC6_piControl_r1i1p1f1_gn_340001-349912.nc rlut_Amon_MIROC6_piControl_r1i1p1f1_gn_350001-359912.nc rlut_Amon_MIROC6_piControl_r1i1p1f1_gn_360001-369912.nc rlut_Amon_MIROC6_piControl_r1i1p1f1_gn_370001-379912.nc rlut_Amon_MIROC6_piControl_r1i1p1f1_gn_380001-389912.nc rlut_Amon_MIROC6_piControl_r1i1p1f1_gn_390001-399912.nc rsdt_Amon_MIROC6_1pctCO2_r1i1p1f1_gn_320001-329912.nc rsdt_Amon_MIROC6_1pctCO2_r1i1p1f1_gn_330001-334912.nc rsdt_Amon_MIROC6_abrupt-2xCO2_r1i1p1f1_gn_320001-329912.nc rsdt_Amon_MIROC6_abrupt-2xCO2_r1i1p1f1_gn_330001-339912.nc rsdt_Amon_MIROC6_abrupt-2xCO2_r1i1p1f1_gn_340001-344912.nc rsdt_Amon_MIROC6_abrupt-4xCO2_r1i1p1f1_gn_320001-329912.nc rsdt_Amon_MIROC6_abrupt-4xCO2_r1i1p1f1_gn_330001-334912.nc rsdt_Amon_MIROC6_abrupt-4xCO2_r1i1p1f1_gn_335001-344912.nc rsdt_Amon_MIROC6_historical_r1i1p1f1_gn_185001-194912.nc rsdt_Amon_MIROC6_historical_r1i1p1f1_gn_195001-201412.nc rsdt_Amon_MIROC6_piControl_r1i1p1f1_gn_320001-329912.nc rsdt_Amon_MIROC6_piControl_r1i1p1f1_gn_330001-339912.nc rsdt_Amon_MIROC6_piControl_r1i1p1f1_gn_340001-349912.nc rsdt_Amon_MIROC6_piControl_r1i1p1f1_gn_350001-359912.nc rsdt_Amon_MIROC6_piControl_r1i1p1f1_gn_360001-369912.nc rsdt_Amon_MIROC6_piControl_r1i1p1f1_gn_370001-379912.nc rsdt_Amon_MIROC6_piControl_r1i1p1f1_gn_380001-389912.nc rsdt_Amon_MIROC6_piControl_r1i1p1f1_gn_390001-399912.nc rsut_Amon_MIROC6_1pctCO2_r1i1p1f1_gn_320001-329912.nc rsut_Amon_MIROC6_1pctCO2_r1i1p1f1_gn_330001-334912.nc rsut_Amon_MIROC6_abrupt-2xCO2_r1i1p1f1_gn_320001-329912.nc rsut_Amon_MIROC6_abrupt-2xCO2_r1i1p1f1_gn_330001-339912.nc rsut_Amon_MIROC6_abrupt-2xCO2_r1i1p1f1_gn_340001-344912.nc rsut_Amon_MIROC6_abrupt-4xCO2_r1i1p1f1_gn_320001-329912.nc rsut_Amon_MIROC6_abrupt-4xCO2_r1i1p1f1_gn_330001-334912.nc rsut_Amon_MIROC6_abrupt-4xCO2_r1i1p1f1_gn_335001-344912.nc rsut_Amon_MIROC6_historical_r1i1p1f1_gn_185001-194912.nc rsut_Amon_MIROC6_historical_r1i1p1f1_gn_195001-201412.nc rsut_Amon_MIROC6_piControl_r1i1p1f1_gn_320001-329912.nc rsut_Amon_MIROC6_piControl_r1i1p1f1_gn_330001-339912.nc rsut_Amon_MIROC6_piControl_r1i1p1f1_gn_340001-349912.nc rsut_Amon_MIROC6_piControl_r1i1p1f1_gn_350001-359912.nc rsut_Amon_MIROC6_piControl_r1i1p1f1_gn_360001-369912.nc rsut_Amon_MIROC6_piControl_r1i1p1f1_gn_370001-379912.nc rsut_Amon_MIROC6_piControl_r1i1p1f1_gn_380001-389912.nc rsut_Amon_MIROC6_piControl_r1i1p1f1_gn_390001-399912.nc tas_Amon_MIROC6_1pctCO2_r1i1p1f1_gn_320001-329912.nc tas_Amon_MIROC6_1pctCO2_r1i1p1f1_gn_330001-334912.nc tas_Amon_MIROC6_abrupt-2xCO2_r1i1p1f1_gn_320001-329912.nc tas_Amon_MIROC6_abrupt-2xCO2_r1i1p1f1_gn_330001-339912.nc tas_Amon_MIROC6_abrupt-2xCO2_r1i1p1f1_gn_340001-344912.nc tas_Amon_MIROC6_abrupt-4xCO2_r1i1p1f1_gn_320001-329912.nc tas_Amon_MIROC6_abrupt-4xCO2_r1i1p1f1_gn_330001-334912.nc tas_Amon_MIROC6_abrupt-4xCO2_r1i1p1f1_gn_335001-344912.nc tas_Amon_MIROC6_hist-nat_r1i1p1f1_gn_185001-194912.nc tas_Amon_MIROC6_hist-nat_r1i1p1f1_gn_195001-202012.nc tas_Amon_MIROC6_historical_r1i1p1f1_gn_185001-194912.nc tas_Amon_MIROC6_historical_r1i1p1f1_gn_195001-201412.nc tas_Amon_MIROC6_historical_r2i1p1f1_gn_185001-194912.nc tas_Amon_MIROC6_historical_r2i1p1f1_gn_195001-201412.nc tas_Amon_MIROC6_historical_r3i1p1f1_gn_185001-194912.nc tas_Amon_MIROC6_historical_r3i1p1f1_gn_195001-201412.nc tas_Amon_MIROC6_historical_r4i1p1f1_gn_185001-194912.nc tas_Amon_MIROC6_historical_r4i1p1f1_gn_195001-201412.nc tas_Amon_MIROC6_piControl_r1i1p1f1_gn_320001-329912.nc tas_Amon_MIROC6_piControl_r1i1p1f1_gn_330001-339912.nc tas_Amon_MIROC6_piControl_r1i1p1f1_gn_340001-349912.nc tas_Amon_MIROC6_piControl_r1i1p1f1_gn_350001-359912.nc tas_Amon_MIROC6_piControl_r1i1p1f1_gn_360001-369912.nc tas_Amon_MIROC6_piControl_r1i1p1f1_gn_370001-379912.nc tas_Amon_MIROC6_piControl_r1i1p1f1_gn_380001-389912.nc tas_Amon_MIROC6_piControl_r1i1p1f1_gn_390001-399912.nc tas_Amon_MIROC6_ssp119_r1i1p1f1_gn_201501-210012.nc tas_Amon_MIROC6_ssp245_r1i1p1f1_gn_201501-210012.nc tas_Amon_MIROC6_ssp370_r1i1p1f1_gn_201501-210012.nc tas_Amon_MIROC6_ssp460_r1i1p1f1_gn_201501-210012.nc tas_Amon_MIROC6_ssp585_r1i1p1f1_gn_201501-210012.nc
Next, process the temperature and radiation data to generate global-mean time series, outputting the results as csv files. This can be done using the provided Python script:
python process_cmip_data.py [--model_id MIROC6]The --model_id argument allows for specifying a particular climate model, with MIROC6 as the default.
The processed data will be stored in ./data_processed directory.
ls ./data_processed | grep MIROC6rlut_MIROC6_1pctCO2_r1i1p1f1.csv rlut_MIROC6_abrupt-2xCO2_r1i1p1f1.csv rlut_MIROC6_abrupt-4xCO2_r1i1p1f1.csv rlut_MIROC6_piControl_r1i1p1f1.csv rsdt_MIROC6_1pctCO2_r1i1p1f1.csv rsdt_MIROC6_abrupt-2xCO2_r1i1p1f1.csv rsdt_MIROC6_abrupt-4xCO2_r1i1p1f1.csv rsdt_MIROC6_piControl_r1i1p1f1.csv rsut_MIROC6_1pctCO2_r1i1p1f1.csv rsut_MIROC6_abrupt-2xCO2_r1i1p1f1.csv rsut_MIROC6_abrupt-4xCO2_r1i1p1f1.csv rsut_MIROC6_piControl_r1i1p1f1.csv tas_MIROC6_1pctCO2_r1i1p1f1.csv tas_MIROC6_abrupt-2xCO2_r1i1p1f1.csv tas_MIROC6_abrupt-4xCO2_r1i1p1f1.csv tas_MIROC6_historical_r1i1p1f1.csv tas_MIROC6_piControl_r1i1p1f1.csv tas_MIROC6_ssp119_r1i1p1f1.csv tas_MIROC6_ssp245_r1i1p1f1.csv tas_MIROC6_ssp370_r1i1p1f1.csv tas_MIROC6_ssp460_r1i1p1f1.csv tas_MIROC6_ssp585_r1i1p1f1.csv
After preprocessing the CMIP data, visualizations can be generated to examine the results. For example:
python plot_experiment.py [--model_id MIROC6]
python plot_historical_tas.py [--model_id MIROC6]
python plot_scenario_tas.py [--model_id MIROC6]
The generated figures will be stored in ./output directory.
The model is a two-layer energy balance model a la Cummins et al. (2020).
I calibrate the model based on the abrupt-4xCO2 experiment:
python calibrate_emulator.py [MIROC6]--- BFGS fvalue: 5.815386799558981 (attempt 1) status: Success in 3.691662073135376 seconds message: Optimization terminated successfully. estimated parameters (vs initial guess): 1.8550 +-0.4758 (1.9937) 4.5689 +-0.8068 (5.1617) 354.4277 +-62.7558 (356.3721) 1.5986 +-0.1792 (1.4604) 1.0627 +-0.0882 (1.0577) 0.4566 +-0.1749 (0.3517) 0.7309 +-0.0923 (0.7535) 0.9718 +-0.1588 (1.0787) 0.0013 +-1.8005 (0.5765) 10.4458 +-0.8570 (9.6197) --- SLSQP fvalue: 5.815386799558981 (attempt 1) status: Success in 0.04993581771850586 seconds message: Optimization terminated successfully estimated parameters (vs initial guess): 1.8550 +-0.4758 (1.9937) 4.5689 +-0.8068 (5.1617) 354.4277 +-62.7558 (356.3721) 1.5986 +-0.1792 (1.4604) 1.0627 +-0.0882 (1.0577) 0.4566 +-0.1749 (0.3517) 0.7309 +-0.0923 (0.7535) 0.9718 +-0.1588 (1.0787) 0.0013 +-1.8005 (0.5765) 10.4458 +-0.8570 (9.6197) --- Nelder-Mead (Best method) fvalue: 5.815386254608917 (attempt 1) status: Success in 12.236296892166138 seconds message: Optimization terminated successfully. estimated parameters (vs initial guess): 1.8550 +-0.4758 (1.9937) 4.5689 +-0.8068 (5.1617) 354.4272 +-62.7558 (356.3721) 1.5986 +-0.1792 (1.4604) 1.0627 +-0.0882 (1.0577) 0.4566 +-0.1749 (0.3517) 0.7309 +-0.0923 (0.7535) 0.9718 +-0.1588 (1.0787) 0.0000 +-1.8005 (0.5765) 10.4459 +-0.8570 (9.6197)
The estimated parameter values of the model will be stored in ./output directory.
cat ./output/parameter_MIROC6_abrupt-4xCO2_r1i1p1f1.csvparameter,value gamma,1.854986624583049 chi1,4.5689113459737865 chi2,354.4271598484283 kappa1,1.5986036054916877 kappa2,1.0626953280905416 epsilon,0.45658815881797343 sigma1,0.7309323634764261 sigma2,0.9717890994755948 sigma3,5.887261283027005e-06 Fbar,10.445851486241933
Evaluate the internal validity (against abrupt-4xCO2)
and the external validity (against historical, ssp119, ssp245, ssp370, ssp460, ssp585)
of the calibrated model:
python evaluate_emulator.py [--model_id MIROC6]
For the carbon cycle, linear models adequately represent impulse responses given a background concentration level. However, these responses are known to be influenced by background concentration levels (Joos et al., 2013), meaning that nonlinear feedback mechanisms must also be considered.
I address this through a two-step approach. First, I employ a linear model and calibrate the model parameters based on the impulse response experiment of Joos et al. (2013). Then, taking the calibrated carbon cycle parameters given, I calibrate the scaling factor as a function of cumulative carbon emission to capture the non-linear feedback effects.
Let us begin with a four-layer linear carbon cycle model
and calibrate its parameters using experimental data (PI100, PD100, PI500) from Joos et al. (2013).
python calibrate_co2_cycle_linear.py
This generates three different sets of parameter values, each optimized for a particular background concentration level.
ls ./output/parameter_co2_cycle_linear*.csv./output/parameter_co2_cycle_linear_PD100.csv ./output/parameter_co2_cycle_linear_PI100.csv ./output/parameter_co2_cycle_linear_PI5000.csv
For example:
cat ./output/parameter_co2_cycle_linear_PD100.csvvar_id,delta21,delta31,delta12,delta32,delta13,delta43,delta34 co2,0.01854978884905409,0.01206756382904478,0.030861326590195495,1.4799916209866694e-14,0.011842834050843291,0.0022372207337764587,0.0011589024764729959
With the carbon cycle matrix, say A, calibrated as described, I incorporate a non-linear scaling factor, exp(-gamma0 - gamma1 * M) * A, where M represents cumulative carbon emissions. The parameter gamma1 is intended to be positive, reflecting the potential weakening of the carbon cycle as the total carbon in the system increases.
When feedback adjustments are disabled (gamma0=gamma1=0),
none of the calibrated linear models (PD100, PI100, PI5000) accurately reproduce the output of historical or scenario experiments (as compiled in RCMIP).
However, by calibrating gamma0 and gamma1, the non-linear model is able to effectively emulate the results from these experiments.
python calibrate_co2_cycle_nonlinear.py
This process yields three distinct sets of non-linear parameter values, derived from three different sets of linear model parameters (PD100, PI100, PI5000).
ls ./output/parameter_co2_cycle_nonlinear*.csv./output/parameter_co2_cycle_nonlinear_PD100.csv ./output/parameter_co2_cycle_nonlinear_PI100.csv ./output/parameter_co2_cycle_nonlinear_PI5000.csv
For example:
cat ./output/parameter_co2_cycle_nonlinear_PD100.csvvar_id,gamma0,gamma1 co2,-0.6691877083124073,0.0002591773687699655
The PI5000-based calibration appears to yield the most satisfying results among the three variants evaluated.
For CH4 and N2O, I use linear gas cycle models without considering interactions between the gases, which is not perfect, but appear sufficient for the moment.
python calibrate_gas_cycle.py
I employ a standard power function, forcing = phi * (u^zeta - 1)/zeta, where u represents the concentration of a well-mixed greenhouse gas (CO2, CH4, or N2O). The parameters phi and zeta are calibrated for each gas, individually, to align with RCMIP data.
For CO2, zeta is manually set to 0 (resulting in a logarithmic relationship) because this provides a good fit and sometimes simplifies analytical calculations.
python calibrate_forcing.py
Obtain the RFF socio-economic projections
and extract the pop_income and emissions datasets to the ./data_raw/RFF directory.
ls ./data_raw/RFF/emissions/*.csv./data_raw/RFF/emissions/rffsp_ch4_emissions.csv ./data_raw/RFF/emissions/rffsp_co2_emissions.csv ./data_raw/RFF/emissions/rffsp_n2o_emissions.csv
ls ./data_raw/RFF/pop_income/*.feather./data_raw/RFF/pop_income/rffsp_pop_income_run_1.feather ./data_raw/RFF/pop_income/rffsp_pop_income_run_2.feather ./data_raw/RFF/pop_income/rffsp_pop_income_run_3.feather ... ./data_raw/RFF/pop_income/rffsp_pop_income_run_10000.feather
Then integrate historical data with the RFF projections and produce a csv file for each of the 10,000 simulated samples.
python process_gdp_pop_data.py
python process_emission_data.pyThe output files will be stored in ./data_processed/gdp_pop and ./data_processed/emissions directories.
The following script allows us to visually inspect the generated files:
python plot_rff_projections.py
python calibrate_savings_rate.py
Download the replication package for Bauer and Rudebusch (2023) from Michael Bauer’s website.
Estimate the parameters of their interest rate model using the script estimate_uc.r, and save the results to ./data_processed/uc_estimates_y10.RData.
Then, use the estimated model to generate future interest rate projections:
Rscript simulate_interest_rate.r
I align the distributions of projected interest rate (Bauer and Rudebusch) and per-capita income growth (RFF) by selecting preference parameters (rho and eta) that minimize the discrepancy between these projections.
python calibrate_rho_eta.py
Store the climate damage point estimates from Barrage and Nordhaus (2024) in /data_raw/Barrage2024/dice2023.csv.
The damage function parameters are selected to globally replicate these estimated damages:
python calibrate_damage_function.py
I estimate the social cost of CO2, CH4, and N2O in several steps.
Create future projections of atmospheric concentrations and radiative forcing by combining RFF emission projections with gas cycle models:
python calculate_scc_1_emis_conc_forc.py
This step is independent of the climate emulator selected.
10,000 concentration trajectories are simulated for each gas, both with and without a unit emission pulse (1 GtCO2, 1 MtCH4, or 1 MtN2O) introduced in 2020. In rare cases, simulated gas concentrations reach negative values; these instances are replaced with a near-zero value.
Next, using an emulator calibrated against CMIP climate models, 10,000 temperature pathways are simulated, both with and without a unit emission pulse of either 1 GtCO2, 1 MtCH4, or 1 MtN2O.
python calculate_scc_2_forc_temp.pypython calculate_scc_3_delta_tas.pypython calculate_scc_4_damage.pypython calculate_scc_5_present_value.pysocial cost of co2 (MIROC6): 75.08007257153537 (USD/tCO2) social cost of ch4 (MIROC6): 377.2429176555695 (USD/tCH4) social cost of n2o (MIROC6): 12971.253584056005 (USD/tN2O)
- Barrage, L., & Nordhaus, W. (2024). Policies, projections, and the social cost of carbon: Results from the DICE-2023 model. Proceedings of the National Academy of Sciences of the United States of America, 121(13), e2312030121.
- Bauer, M. D., & Rudebusch, G. D. (2023). The rising cost of climate change: Evidence from the bond market. The Review of Economics and Statistics, 105(5), 1255–1270.
- Cummins, D. P., Stephenson, D. B., & Stott, P. A. (2020). Optimal Estimation of Stochastic Energy Balance Model Parameters. Journal of Climate, 33(18), 7909–7926.
- Joos, F., Roth, R., Fuglestvedt, J. S., Peters, G. P., Enting, I. G., Bloh, W. von, Brovkin, V., Burke, E. J., Eby, M., Edwards, N. R., & Others. (2013). Carbon dioxide and climate impulse response functions for the computation of greenhouse gas metrics: a multi-model analysis. Atmospheric Chemistry and Physics, 13(5), 2793–2825.
- Meinshausen, M., Nicholls, Z. R. J., Lewis, J., Gidden, M. J., Vogel, E., Freund, M., Beyerle, U., Gessner, C., Nauels, A., Bauer, N., Canadell, J. G., Daniel, J. S., John, A., Krummel, P. B., Luderer, G., Meinshausen, N., Montzka, S. A., Rayner, P. J., Reimann, S., Smith, S. J., van den Berg, M., Velders, G. J. M., Vollmer, M. K., and Wang, R. H. J. (2020). The shared socio-economic pathway (SSP) greenhouse gas concentrations and their extensions to 2500. Geoscientific Model Development, 13, 3571–3605.
- MPD version 2023: Bolt, Jutta and Jan Luiten van Zanden (2024). Maddison style estimates of the evolution of the world economy: A new 2023 update. Journal of Economic Surveys, 1–41.
- Osborn, T.J., Jones, P.D., Lister, D.H., Morice, C.P., Simpson, I.R., Winn, J.P., Hogan, E., and Harris, I.C., (2021). Land surface air temperature variations across the globe updated to 2019: the CRUTEM5 dataset. Journal of Geophysical Research: Atmospheres. 126, e2019JD032352,
- Kevin Rennert, Brian C. Prest, William A. Pizer, Richard G. Newell, David Anthoff, Cora Kingdon, Lisa Rennels, Roger Cooke, Adrian E. Raftery, Hana Ševčíková, & Frank Errickson. (2022). Resources for the Future Socioeconomic Projections (RFF-SPs) (Version 5). Zenodo.
- Smith, Christopher J. (2019, October 21). Effective Radiative Forcing from Shared Socioeconomic Pathways (Version v0.3.1). Zenodo.
- Tatebe, Hiroaki; Watanabe, Masahiro (2018). MIROC6 model output prepared for CMIP6. Earth System Grid Federation.
- World Bank. Gross savings (% of GDP) World Development Indicators. The World Bank Group.
- Zebedee Nicholls, & Jared Lewis. (2021). Reduced Complexity Model Intercomparison Project (RCMIP) protocol (v5.1.0). Zenodo.