Climate Prediction

Climate Prediction Presentation [Nathan Urban, LANL]


Hands-On

Simple Climate Model

# Simple climate model (Nathan Urban, nurban@lanl.gov)

import matplotlib.pyplot as plt; import numpy as np
infrared_loss = 1.7 # watts per square meter per degree Celsiusupper_ocean_depth = 100 # meters
water_heat_capacity = 4184000 # Joules per cubic meterheat_capacity = water_heat_capacity * upper_ocean_depthdt = 31557600 # seconds per year
# forcing for doubled CO2number_of_years = 51Forcing = np.full(number_of_years, 3.7) # watts per square meter for doubled CO2Time = np.arange(0, number_of_years, 1)
Temperature = np.zeros(number_of_years)
for year in range(0, number_of_years-1): heat_in = Forcing[year] * dt heat_out = infrared_loss * Temperature[year] * dt heat = heat_in - heat_out
 temperature_change = heat / heat_capacity Temperature[year+1] = Temperature[year] + temperature_change
plt.plot(Time, Temperature)plt.xlabel('Year')plt.ylabel('Temperature')plt.show()


Emissions Control with a Simple Climate Model


# Emissions control with a simple climate model (Nathan Urban, nurban@lanl.gov)
import matplotlib.pyplot as pltimport numpy as np
Year = np.arange(1850, 2101)number_of_years = len(Year)
# historical and future projected CO2 emissions (gigatons carbon per year)CO2_Emissions = np.array([0.508, 0.501, 0.556, 0.556, 0.565, 0.565, 0.575, 0.581, 0.587, 0.597, 0.610, 0.623, 0.571, 0.578, 0.586, 0.594, 0.597, 0.604, 0.607, 0.613, 0.617, 0.644, 0.740, 0.761, 0.758, 0.778, 0.788, 0.797, 0.805, 0.825, 0.857, 0.892, 0.861, 0.881, 0.887, 0.892, 0.898, 0.911, 0.942, 0.940, 0.968, 0.982, 0.997, 0.994, 1.023, 1.050, 1.065, 1.090, 1.116, 1.159, 1.187, 1.255, 1.269, 1.344, 1.372, 1.433, 1.504, 1.588, 1.561, 1.602, 1.641, 1.615, 1.627, 1.666, 1.565, 1.544, 1.609, 1.666, 1.650, 1.526, 1.653, 1.560, 1.593, 1.723, 1.720, 1.733, 1.744, 1.859, 1.865, 1.969, 1.949, 1.843, 1.659, 1.702, 1.770, 1.822, 1.931, 1.989, 1.924, 1.971, 2.066, 2.084, 2.111, 2.155, 2.151, 1.929, 2.080, 2.262, 2.344, 2.301, 2.522, 2.857, 2.906, 2.946, 3.017, 3.234, 3.423, 3.542, 3.650, 3.666, 3.769, 3.838, 3.936, 4.092, 4.260, 4.409, 4.599, 4.709, 4.823, 5.038, 5.273, 5.302, 5.446, 5.668, 5.683, 5.647, 5.974, 6.123, 6.195, 6.452, 6.357, 6.213, 6.317, 6.336, 6.550, 6.719, 6.897, 7.047, 7.281, 7.413, 7.463, 7.616, 7.422, 7.422, 7.532, 7.647, 7.761, 7.837, 7.803, 7.752, 7.884, 8.028, 8.181, 8.512, 8.915, 9.166, 9.326, 9.487, 9.648, 9.809, 9.969, 10.217, 10.464, 10.712, 10.959, 11.207, 11.454, 11.702, 11.949, 12.197, 12.444, 12.655, 12.866, 13.077, 13.288, 13.499, 13.710, 13.921, 14.132, 14.343, 14.554, 14.842, 15.130, 15.417, 15.705, 15.993, 16.281, 16.568, 16.856, 17.144, 17.432, 17.767, 18.101, 18.436, 18.771, 19.106, 19.441, 19.776, 20.111, 20.446, 20.781, 21.112, 21.444, 21.776, 22.107, 22.439, 22.770, 23.102, 23.434, 23.765, 24.097, 24.325, 24.553, 24.780, 25.008, 25.236, 25.463, 25.691, 25.919, 26.146, 26.374, 26.508, 26.642, 26.776, 26.911, 27.045, 27.179, 27.313, 27.447, 27.581, 27.715, 27.797, 27.878, 27.960, 28.041, 28.123, 28.204, 28.286, 28.368, 28.449, 28.531, 28.559, 28.588, 28.617, 28.645, 28.674, 28.702, 28.731, 28.760, 28.788, 28.817])
##### CHANGE THE CLIMATE ###### we will change CO2 emissions between 2020 ("year 170") to 2100 ("year 250") ... years start counting from 1850
# REMOVE '#' FROM THE LINE BELOW to hold CO2 emissions fixed at 2020 levels#CO2_Emissions[170:251] = CO2_Emissions[170]
# REMOVE '#' FROM THE LINE BELOW to decrease CO2 emissions to zero by 2100#CO2_Emissions[170:251] = np.linspace(CO2_Emissions[170], 0, 81)
# REMOVE '#' FROM THE LINE BELOW to stop emitting CO2 in 2020 (total decarbonization)#CO2_Emissions[170:251] = 0
#### END CHANGE THE CLIMATE #####
# convert CO2 emissions to concentrations (parts per million)# assume a fraction between 0.43 and 0.71 of emissions accumulate in atmosphereairborne_fraction = np.zeros(251) # fraction of emitted CO2 that remains in the atmosphere (rest is absorbed by plants and ocean)airborne_fraction[0:150] = 0.43 # airborne fraction is 0.43 until year 2000airborne_fraction[150:251] = np.linspace(0.43, 0.71, 101) # airborne fraction increases from 0.43 to 0.71 during 2000-2100gtc_per_ppm = 2.13 # gigatons carbon per ppm of CO2CO2_Concentration = 285 + np.cumsum(airborne_fraction*CO2_Emissions) / gtc_per_ppm # 285 = CO2 concentration in 1850 (ppm)
# radiative forcing for CO2F2xCO2 = 3.7 # watts per square meter of forcing when CO2 is doubled (no feedback effects)Forcing_CO2 = F2xCO2 * np.log2(CO2_Concentration/277) # 277 = pre-industrial CO2 concentration (ppm)
# radiative forcing for non-CO2 greenhouse gases, aerosols, solar cycle, and volcanoesForcing_OtherGHG = np.array([0.056, 0.057, 0.059, 0.060, 0.061, 0.062, 0.063, 0.064, 0.066, 0.067, 0.068, 0.069, 0.070, 0.071, 0.073, 0.074, 0.075, 0.076, 0.078, 0.079, 0.080, 0.081, 0.082, 0.083, 0.085, 0.086, 0.087, 0.088, 0.089, 0.090, 0.092, 0.093, 0.095, 0.096, 0.097, 0.099, 0.100, 0.101, 0.102, 0.104, 0.105, 0.107, 0.109, 0.110, 0.111, 0.113, 0.114, 0.116, 0.117, 0.119, 0.120, 0.123, 0.125, 0.128, 0.131, 0.133, 0.136, 0.139, 0.142, 0.145, 0.148, 0.151, 0.154, 0.157, 0.161, 0.164, 0.168, 0.171, 0.175, 0.179, 0.183, 0.187, 0.190, 0.194, 0.198, 0.202, 0.205, 0.209, 0.213, 0.216, 0.220, 0.223, 0.226, 0.230, 0.233, 0.236, 0.240, 0.243, 0.246, 0.249, 0.252, 0.256, 0.259, 0.262, 0.266, 0.269, 0.273, 0.278, 0.282, 0.287, 0.292, 0.297, 0.303, 0.309, 0.315, 0.321, 0.328, 0.335, 0.343, 0.351, 0.359, 0.367, 0.377, 0.386, 0.397, 0.408, 0.420, 0.433, 0.446, 0.461, 0.476, 0.492, 0.509, 0.527, 0.546, 0.566, 0.586, 0.607, 0.627, 0.648, 0.667, 0.688, 0.709, 0.729, 0.748, 0.767, 0.785, 0.803, 0.824, 0.845, 0.864, 0.881, 0.895, 0.905, 0.914, 0.922, 0.929, 0.936, 0.943, 0.950, 0.955, 0.959, 0.964, 0.969, 0.973, 0.976, 0.980, 0.984, 0.989, 0.995, 1.002, 1.009, 1.017, 1.025, 1.034, 1.043, 1.053, 1.063, 1.074, 1.085, 1.096, 1.107, 1.119, 1.130, 1.142, 1.154, 1.165, 1.177, 1.189, 1.200, 1.212, 1.224, 1.236, 1.248, 1.260, 1.272, 1.285, 1.297, 1.309, 1.321, 1.334, 1.345, 1.357, 1.369, 1.380, 1.392, 1.403, 1.415, 1.427, 1.439, 1.451, 1.463, 1.475, 1.486, 1.498, 1.509, 1.521, 1.532, 1.542, 1.553, 1.564, 1.575, 1.585, 1.595, 1.605, 1.614, 1.623, 1.632, 1.641, 1.649, 1.657, 1.665, 1.673, 1.681, 1.689, 1.696, 1.704, 1.712, 1.719, 1.726, 1.733, 1.741, 1.748, 1.755, 1.763, 1.770, 1.777, 1.785, 1.792, 1.800, 1.807, 1.814, 1.822, 1.829, 1.836, 1.843, 1.849, 1.855, 1.861, 1.867, 1.873])Forcing_Aerosol = np.array([-0.092, -0.094, -0.099, -0.105, -0.108, -0.111, -0.113, -0.116, -0.117, -0.119, -0.128, -0.141, -0.140, -0.136, -0.139, -0.143, -0.152, -0.161, -0.165, -0.163, -0.164, -0.173, -0.194, -0.211, -0.215, -0.215, -0.215, -0.211, -0.209, -0.213, -0.226, -0.235, -0.234, -0.237, -0.240, -0.238, -0.238, -0.241, -0.251, -0.257, -0.261, -0.269, -0.274, -0.279, -0.283, -0.289, -0.294, -0.296, -0.300, -0.305, -0.311, -0.317, -0.324, -0.331, -0.339, -0.346, -0.355, -0.364, -0.367, -0.368, -0.373, -0.376, -0.377, -0.380, -0.379, -0.375, -0.378, -0.383, -0.385, -0.382, -0.382, -0.382, -0.379, -0.391, -0.400, -0.402, -0.405, -0.410, -0.415, -0.420, -0.421, -0.416, -0.411, -0.411, -0.415, -0.419, -0.424, -0.428, -0.424, -0.422, -0.428, -0.435, -0.441, -0.446, -0.453, -0.458, -0.466, -0.479, -0.493, -0.503, -0.514, -0.531, -0.547, -0.561, -0.576, -0.590, -0.606, -0.622, -0.635, -0.646, -0.656, -0.672, -0.691, -0.705, -0.720, -0.733, -0.746, -0.761, -0.775, -0.786, -0.790, -0.795, -0.810, -0.828, -0.845, -0.864, -0.882, -0.900, -0.913, -0.924, -0.941, -0.961, -0.977, -0.993, -1.006, -1.016, -1.027, -1.038, -1.049, -1.056, -1.061, -1.064, -1.068, -1.068, -1.064, -1.057, -1.051, -1.058, -1.073, -1.085, -1.099, -1.108, -1.111, -1.114, -1.118, -1.117, -1.108, -1.095, -1.083, -1.069, -1.060, -1.057, -1.057, -1.057, -1.058, -1.058, -1.058, -1.059, -1.059, -1.059, -1.059, -1.056, -1.052, -1.048, -1.045, -1.041, -1.037, -1.034, -1.030, -1.026, -1.021, -1.012, -1.002, -0.991, -0.981, -0.970, -0.960, -0.949, -0.939, -0.928, -0.918, -0.908, -0.900, -0.891, -0.882, -0.873, -0.864, -0.855, -0.846, -0.836, -0.829, -0.823, -0.818, -0.814, -0.810, -0.805, -0.801, -0.796, -0.792, -0.787, -0.783, -0.779, -0.776, -0.772, -0.769, -0.766, -0.762, -0.759, -0.756, -0.752, -0.749, -0.746, -0.744, -0.741, -0.739, -0.736, -0.733, -0.731, -0.728, -0.726, -0.722, -0.717, -0.711, -0.704, -0.698, -0.692, -0.685, -0.679, -0.673, -0.666, -0.661, -0.657, -0.654, -0.651, -0.648, -0.646, -0.643, -0.640, -0.637, -0.634, -0.632])Forcing_Solar = np.array([0.047, 0.035, 0.023, 0.005, -0.014, -0.026, -0.025, -0.007, 0.025, 0.055, 0.066, 0.052, 0.031, 0.014, 0.001, -0.012, -0.024, -0.025, -0.003, 0.036, 0.064, 0.067, 0.052, 0.027, -0.000, -0.021, -0.032, -0.035, -0.037, -0.030, -0.009, 0.011, 0.021, 0.022, 0.011, -0.011, -0.032, -0.044, -0.051, -0.053, -0.040, -0.012, 0.013, 0.034, 0.044, 0.030, 0.002, -0.019, -0.030, -0.037, -0.045, -0.053, -0.046, -0.016, 0.009, 0.013, 0.016, 0.018, 0.014, 0.001, -0.019, -0.034, -0.041, -0.039, -0.019, 0.020, 0.058, 0.078, 0.072, 0.043, 0.014, -0.008, -0.021, -0.022, -0.009, 0.015, 0.040, 0.055, 0.051, 0.041, 0.033, 0.016, -0.004, -0.014, 0.000, 0.041, 0.082, 0.090, 0.081, 0.074, 0.064, 0.050, 0.031, 0.014, 0.021, 0.050, 0.081, 0.112, 0.129, 0.116, 0.081, 0.049, 0.032, 0.020, 0.022, 0.057, 0.125, 0.178, 0.180, 0.151, 0.116, 0.072, 0.039, 0.027, 0.027, 0.038, 0.064, 0.092, 0.108, 0.118, 0.108, 0.087, 0.072, 0.051, 0.025, 0.010, 0.019, 0.060, 0.123, 0.174, 0.193, 0.180, 0.142, 0.100, 0.055, 0.023, 0.022, 0.048, 0.105, 0.163, 0.174, 0.155, 0.125, 0.085, 0.049, 0.026, 0.020, 0.043, 0.093, 0.146, 0.181, 0.191, 0.174, 0.130, 0.086, 0.059, 0.043, 0.032, 0.026, 0.073, 0.146, 0.181, 0.191, 0.174, 0.130, 0.086, 0.059, 0.043, 0.032, 0.026, 0.073, 0.146, 0.181, 0.191, 0.174, 0.130, 0.086, 0.059, 0.043, 0.032, 0.026, 0.073, 0.146, 0.181, 0.191, 0.174, 0.130, 0.086, 0.059, 0.043, 0.032, 0.026, 0.073, 0.146, 0.181, 0.191, 0.174, 0.130, 0.086, 0.059, 0.043, 0.032, 0.026, 0.073, 0.146, 0.181, 0.191, 0.174, 0.130, 0.086, 0.059, 0.043, 0.032, 0.026, 0.073, 0.146, 0.181, 0.191, 0.174, 0.130, 0.086, 0.059, 0.043, 0.032, 0.026, 0.073, 0.146, 0.181, 0.191, 0.174, 0.130, 0.086, 0.059, 0.043, 0.032, 0.026, 0.073, 0.146, 0.181, 0.191, 0.174, 0.130, 0.086, 0.059, 0.043, 0.032, 0.026, 0.073, 0.146, 0.181, 0.191])Forcing_Volcanic = np.array([0.192, 0.193, 0.209, 0.224, 0.232, 0.101, -0.357, -0.610, -0.298, 0.028, 0.151, 0.140, 0.055, 0.068, 0.151, 0.199, 0.219, 0.229, 0.231, 0.223, 0.219, 0.219, 0.207, 0.193, 0.201, 0.194, 0.153, 0.148, 0.180, 0.202, 0.216, 0.223, 0.111, -0.803, -1.537, -0.944, -0.449, -0.321, -0.211, -0.278, -0.332, -0.262, -0.088, 0.060, 0.169, 0.162, 0.017, -0.065, 0.038, 0.151, 0.202, 0.192, -0.257, -0.651, -0.356, 0.002, 0.100, 0.076, 0.088, 0.154, 0.189, 0.195, -0.049, -0.157, 0.056, 0.159, 0.188, 0.196, 0.202, 0.195, 0.111, 0.096, 0.181, 0.202, 0.189, 0.189, 0.197, 0.195, 0.146, 0.096, 0.122, 0.146, 0.115, 0.124, 0.162, 0.169, 0.181, 0.179, 0.163, 0.164, 0.184, 0.193, 0.171, 0.166, 0.189, 0.200, 0.207, 0.199, 0.197, 0.189, 0.187, 0.191, 0.184, 0.176, 0.181, 0.203, 0.215, 0.226, 0.231, 0.227, 0.152, 0.061, -0.014, -0.495, -0.853, -0.532, -0.150, -0.022, -0.220, -0.296, -0.055, 0.114, 0.153, 0.106, -0.041, -0.141, -0.017, 0.134, 0.109, 0.094, 0.148, 0.080, -0.591, -0.876, -0.320, -0.023, 0.035, 0.066, 0.106, 0.129, 0.127, -0.815, -1.420, -0.682, -0.132, 0.059, 0.124, 0.153, 0.188, 0.216, 0.229, 0.232, 0.232, 0.232, 0.232, 0.184, 0.068, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000])
# total radiative forcing (in watts per square meter)Forcing = Forcing_CO2 + Forcing_OtherGHG + Forcing_Aerosol + Forcing_Solar + Forcing_Volcanic
# these are historical global temperature measurements (degrees Celsius above pre-industrial temperature)Observation_Year = np.array([1850, 1851, 1852, 1853, 1854, 1855, 1856, 1857, 1858, 1859, 1860, 1861, 1862, 1863, 1864, 1865, 1866, 1867, 1868, 1869, 1870, 1871, 1872, 1873, 1874, 1875, 1876, 1877, 1878, 1879, 1880, 1881, 1882, 1883, 1884, 1885, 1886, 1887, 1888, 1889, 1890, 1891, 1892, 1893, 1894, 1895, 1896, 1897, 1898, 1899, 1900, 1901, 1902, 1903, 1904, 1905, 1906, 1907, 1908, 1909, 1910, 1911, 1912, 1913, 1914, 1915, 1916, 1917, 1918, 1919, 1920, 1921, 1922, 1923, 1924, 1925, 1926, 1927, 1928, 1929, 1930, 1931, 1932, 1933, 1934, 1935, 1936, 1937, 1938, 1939, 1940, 1941, 1942, 1943, 1944, 1945, 1946, 1947, 1948, 1949, 1950, 1951, 1952, 1953, 1954, 1955, 1956, 1957, 1958, 1959, 1960, 1961, 1962, 1963, 1964, 1965, 1966, 1967, 1968, 1969, 1970, 1971, 1972, 1973, 1974, 1975, 1976, 1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018])Observed_Temperature = np.array([0.031, 0.186, 0.176, 0.135, 0.156, 0.132, 0.046, -0.057, -0.063, 0.120, 0.061, -0.003, -0.120, 0.126, -0.090, 0.125, 0.153, 0.083, 0.166, 0.142, 0.128, 0.069, 0.177, 0.100, 0.036, 0.009, 0.020, 0.329, 0.439, 0.174, 0.177, 0.204, 0.191, 0.108, -0.005, 0.015, 0.037, -0.014, 0.097, 0.233, -0.012, 0.074, -0.051, -0.069, -0.006, 0.014, 0.218, 0.198, -0.008, 0.115, 0.201, 0.145, 0.002, -0.075, -0.116, 0.027, 0.121, -0.061, -0.107, -0.118, -0.086, -0.140, -0.033, -0.020, 0.160, 0.263, 0.021, -0.064, 0.071, 0.129, 0.157, 0.217, 0.102, 0.128, 0.110, 0.189, 0.296, 0.194, 0.198, 0.054, 0.267, 0.317, 0.267, 0.131, 0.273, 0.226, 0.257, 0.378, 0.398, 0.352, 0.418, 0.424, 0.377, 0.400, 0.548, 0.429, 0.333, 0.366, 0.365, 0.330, 0.231, 0.352, 0.432, 0.501, 0.275, 0.214, 0.137, 0.397, 0.450, 0.421, 0.355, 0.442, 0.418, 0.452, 0.181, 0.264, 0.336, 0.330, 0.291, 0.436, 0.377, 0.218, 0.339, 0.466, 0.190, 0.255, 0.163, 0.451, 0.342, 0.461, 0.496, 0.544, 0.415, 0.598, 0.390, 0.374, 0.449, 0.596, 0.602, 0.522, 0.700, 0.658, 0.509, 0.552, 0.612, 0.729, 0.587, 0.794, 0.943, 0.710, 0.698, 0.845, 0.900, 0.909, 0.851, 0.949, 0.910, 0.895, 0.799, 0.910, 0.964, 0.829, 0.874, 0.918, 0.983, 1.167, 1.201, 1.080, 0.999])
# simple climate modelinfrared_loss = 1.7 # watts per square meter, per degree Celsiusupper_ocean_depth = 100 # meterswater_heat_capacity = 4184000 # Joules per cubic meterheat_capacity = water_heat_capacity * upper_ocean_depthdt = 31557600 # seconds per year
Temperature = np.zeros(number_of_years)
for year in range(0, number_of_years-1): heat_in = Forcing[year] * dt heat_out = infrared_loss * Temperature[year] * dt heat = heat_in - heat_out
 temperature_change = heat / heat_capacity Temperature[year+1] = Temperature[year] + temperature_change
# plot resultsplt.subplot(2,3,1)plt.plot(Year, CO2_Emissions, color='blue')plt.xlim([1850, 2100])plt.ylim([0, 30]) plt.title('CO2 emissions') # gigatons carbon per year
plt.subplot(2,3,2)plt.plot(Year, CO2_Concentration, color='blue')plt.xlim([1850, 2100])plt.ylim([0, 1000]) plt.title('CO2 level') # parts per million
plt.subplot(2,3,3)plt.plot(Year, Forcing_CO2, color='blue')plt.xlim([1850, 2100])plt.ylim([-2, 10]) plt.title('CO2 forcing') # watts per square meter
plt.subplot(2,3,4)plt.plot(Year, Forcing, color='blue')plt.xlim([1850, 2100])plt.ylim([-2, 10])plt.title('Total forcing') # watts per square meter
plt.subplot(2,3,5)plt.plot(Observation_Year, Observed_Temperature, color='red', linewidth=0.75)plt.plot(Year, Temperature, color='blue')plt.xlim([1850, 2020])plt.ylim([-0.2, 1.5])plt.title('Temperature') # degrees Celsius
plt.subplot(2,3,6)plt.plot(Observation_Year, Observed_Temperature, color='red', linewidth=0.75)plt.plot(Year, Temperature, color='blue')plt.xlim([1850, 2100])plt.ylim([-0.5, 5])plt.title('Temperature') # degrees Celsius
plt.tight_layout()plt.show()