Climate & environment scientist: MQL recipes

Recipes

Mauna Loa CO2 trend and 2030 projection
Temperature anomalies against a reference decade
Household carbon footprint
Warming trend and a hidden cycle
Zero-dimensional energy balance model after a CO2 doubling
Acceleration of the sea level rise by QR least squares
Planck feedback and climate sensitivity
Hot summers: shift of the mean and return period
Amplitude and phase of the annual temperature cycle
Dates of the remaining 1.5 C carbon budget

executed on MentDB

Mauna Loa CO2 trend and 2030 projection

Fit the annual Mauna Loa CO2 means with a line and an exponential, project 2030, and convert the result to radiative forcing and atmospheric carbon stock.

Pipeline

01 · trendLinear and exponential fitsstat fit→02 · projectCO2 in 2030math eval→03 · forcing5.35 ln C/C0math calc

Script

admin

#1. Mauna Loa annual mean CO2 in ppm, years counted from 2000;
-> "[yr]" "[14, 15, 16, 17, 18, 19, 20, 21, 22, 23]";
-> "[co2]" "[397.12, 399.41, 403.08, 405.22, 407.61, 410.07, 412.44, 414.70, 418.53, 421.08]";
#2. Linear trend in ppm per year and exponential fit for the projection;
-> "[lin]" (stat fit [yr] [co2] "linear");
json load "rx_co2lin" [lin];
-> "[exp]" (stat fit [yr] [co2] "exp");
json load "rx_co2exp" [exp];
-> "[c2030]" (math eval (json select "rx_co2exp" /expression) "{\"x\": 30}");
#3. Radiative forcing 5.35 ln C/C0 against 280 ppm and atmospheric carbon stock 2.124 GtC per ppm;
-> "[dF]" (math calc "5.35*ln([c2030]/280)");
concat "trend = " (math round (json select "rx_co2lin" /b) 3) " ppm/yr, growth = " (math round (math calc (concat "100*" (json select "rx_co2exp" /b))) 3) " %/yr, 2030 = " (math round [c2030] 1) " ppm, forcing = " (math round [dF] 2) " W/m2, stock = " (math round (math calc "[c2030]*2.124") 0) " GtC"

Result

mentdb

trend = 2.616 ppm/yr, growth = 0.64 %/yr, 2030 = 440.1 ppm, forcing = 2.42 W/m2, stock = 935 GtC

Go further

Remove the seasonal cycle of monthly data with dsp detrend and dsp spectrum
Turn the forcing into warming with a climate sensitivity

executed on MentDB

Temperature anomalies against a reference decade

Compute the baseline of a reference decade, the anomalies of the recent decade, and test with a Welch t-test whether the warming is significant.

Pipeline

01 · baselineReference meanstat mean→02 · anomalyRecent minus baselinevector sub→03 · testWelch t-teststat ttest_2

Script

admin

#1. Annual mean temperature of a station in degC, reference decade 1994-2003 and recent decade 2014-2023;
-> "[ref]" "[11.8, 12.1, 11.5, 12.0, 12.3, 12.2, 12.4, 11.9, 12.5, 12.8]";
-> "[now]" "[12.9, 13.2, 13.1, 12.8, 13.4, 13.3, 13.0, 13.6, 13.9, 13.7]";
#2. Baseline and anomalies of the recent years;
-> "[base]" (stat mean [ref]);
-> "[anom]" (vector sub [now] (vector scale "[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]" [base]));
#3. Is the warming significant? Welch t-test between the two decades;
-> "[t]" (stat ttest_2 [now] [ref]);
json load "rx_ttest" [t];
concat "baseline = " (math round [base] 2) " degC, mean anomaly = +" (math round (stat mean [anom]) 2) " degC, max anomaly = +" (math round (stat max [anom]) 2) " degC, t = " (math round (json select "rx_ttest" /t) 2) ", p = " (math round_sig (json select "rx_ttest" /p_value) 2)

Result

mentdb

baseline = 12.15 degC, mean anomaly = +1.14 degC, max anomaly = +1.75 degC, t = 6.93, p = 0.0000018

Go further

Get the trend in degC per decade with stat fit
Smooth the series with dsp moving_average

executed on MentDB

Household carbon footprint

Turn a year of car fuel, flights, electricity, gas heating billed in BTU and methane into CO2 equivalent, then into tonnes of carbon.

Pipeline

01 · convertBTU to kWhunit convert→02 · factorsEmission factors and GWP100math calc→03 · totalSum and largest itemstat sum

Script

admin

#1. Yearly activity data of a household: car, two long-haul flights, electricity, gas heating billed in BTU;
-> "[car]" (math calc "12000*6.5/100*2.31");
-> "[fly]" (math calc "2*6000*0.152");
-> "[elec]" (math calc "4200*0.056");
-> "[gas]" (math calc (concat (unit convert 30.7e6 "BTU" "kWh") "*0.227"));
#2. Methane from waste, 15 kg CH4 converted with the GWP100 of 27.9;
-> "[ch4]" (math calc "15*27.9");
#3. Total in tonnes of CO2 equivalent and in tonnes of carbon, share of the largest item;
-> "[items]" (concat "[" [car] ", " [fly] ", " [elec] ", " [gas] ", " [ch4] "]");
-> "[tot]" (stat sum [items]);
concat "total = " (math round (math calc "[tot]/1000") 2) " t CO2e = " (math round (math calc "[tot]/1000*12/44") 2) " t C, largest item = " (math round (math calc (concat "100*" (stat max [items]) "/[tot]")) 1) " %"

Result

mentdb

total = 6.32 t CO2e = 1.72 t C, largest item = 32.3 %

Go further

Compare with the 2 t per person target
Put the items in a JSON object and sort them with json sort

executed on MentDB

Warming trend and a hidden cycle

An annual temperature anomaly series is fitted by a linear trend with its 95 % confidence interval and the share of variance it explains; the trend is then removed and the periodogram of the residuals reveals a 4-year cycle (ENSO-like).

Pipeline

01 · trendLinear regressionpa rl load_from_json→02 · ciTrend per decade and CIpa rl slope_confidence_interval→03 · variance1 - SSE/SSTpa rl total_sum_squares→04 · cycleDetrend and periodogramdsp periodogram

Script

admin

#1. Annual temperature anomaly 1980-2023 in degC, t = years since 1980;
-> "[data]" (num table "0.018*t+0.12*sin(2*pi*t/4)+0.04*sin(13.7*t)" "t" 0 43 1);
pa rl load_from_json "rx_trend" [data];
#2. Warming trend per decade with its 95 % confidence interval;
-> "[slope]" (math calc (concat "10*" (pa rl slope "rx_trend")));
-> "[ci]" (math calc (concat "10*" (pa rl slope_confidence_interval "rx_trend")));
#3. Share of the variance explained by the trend: 1 - SSE/SST;
-> "[share]" (math calc (concat "1-" (pa rl sum_squared_errors "rx_trend") "/" (pa rl total_sum_squares "rx_trend")));
#4. Remove the trend, then look for a cycle in the residuals;
json load "rx_per" (dsp periodogram (dsp detrend (matrix col [data] 2)) 1);
json load "rx_cyc" (dsp peaks (json select "rx_per" /psd) "{\"distance\": 100}");
-> "[f]" (json select "rx_per" (concat "/f[" (json select "rx_cyc" /indices[0]) "]"));
-> "[r]" (concat "trend = " (math round [slope] 3) " +- " (math round [ci] 3) " degC/decade, r = " (math round (pa rl r "rx_trend") 3) ", residual sd = " (math round (math calc (concat "sqrt(" (pa rl mean_square_error "rx_trend") ")")) 3) " degC, trend explains " (math round (math calc "100*[share]") 1) " % of the variance, residual cycle of " (math round (math calc "1/[f]") 2) " years");
pa rl close "rx_trend";
[r]

Result

mentdb

trend = 0.175 +- 0.022 degC/decade, r = 0.927, residual sd = 0.092 degC, trend explains 85.9 % of the variance, residual cycle of 4 years

Go further

Use dsp welch for a smoother spectrum on a longer series.
Draw the residuals with plot signal (dsp detrend ...).

executed on MentDB

Zero-dimensional energy balance model after a CO2 doubling

Compute the emission temperature of the Earth with the Stefan-Boltzmann law, the radiative forcing of a CO2 doubling (Myhre formula), then integrate the energy balance C dT/dt = F - lambda T of a 100 m ocean mixed layer to see the warming after 10, 30 and 60 years.

Pipeline

01 · constantStefan-Boltzmann sigmaconst get→02 · forcing5.35 ln C/C0 and heat capacitymath calc→03 · integrateC dT/dt = F - lambda Tnum ode→04 · readT at 10, 30, 60 yearsjson select

Script

admin

#1. Emission temperature of the Earth: solar constant 1361 W/m2, albedo 0.30;
-> "[sigma]" (const get "sigma");
-> "[Te]" (math calc "(1361*(1-0.30)/(4*[sigma]))^0.25");
#2. Radiative forcing of an abrupt CO2 doubling, Myhre formula 5.35 ln C/C0;
-> "[F]" (math calc "5.35*ln(560/280)");
#3. Feedback parameter and mixed layer of 100 m: heat capacity per year;
-> "[lambda]" 1.2;
-> "[C]" (math calc "1025*3990*100/3.156e7");
#4. Energy balance model C dT/dt = F - lambda T over 60 years;
json load "rx_ebm" (num ode "{\"T\": \"([F]-[lambda]*T)/[C]\"}" "{\"T\": 0}" 0 60 "{\"points\": 7}");
#5. Warming after 10, 30 and 60 years against the equilibrium;
-> "[ecs]" (/ [F] [lambda]);
concat "Te = " (math round [Te] 1) " K | forcing 2xCO2 = " (math round [F] 2) " W/m2 | ECS = " (math round [ecs] 2) " K, tau = " (math round (/ [C] [lambda]) 1) " years | T10 = " (math round (json select "rx_ebm" "/T/[1]") 2) " K, T30 = " (math round (json select "rx_ebm" "/T/[3]") 2) " K, T60 = " (math round (json select "rx_ebm" "/T/[6]") 2) " K"

Result

mentdb

Te = 254.6 K | forcing 2xCO2 = 3.71 W/m2 | ECS = 3.09 K, tau = 10.8 years | T10 = 1.87 K, T30 = 2.9 K, T60 = 3.08 K

Go further

Draw the warming curve with plot ode on the result of num ode.
Add a deep ocean layer as a second variable of the system to get the slow response.

executed on MentDB

Acceleration of the sea level rise by QR least squares

The satellite record of the global mean sea level is fitted with a quadratic a + b t + c t^2. The script builds the design matrix, solves the least squares problem with a Householder QR decomposition (R b = Q^T y, stable even when t^2 is large), and gets the rate, the acceleration 2c, the residual and the level reached in 2050.

Pipeline

01 · designDesign matrix [1, t, t^2]atom get→02 · qrHouseholder decompositionmatrix qr→03 · solveR b = Q^T ymatrix solve→04 · rateRate, acceleration, 2050vector dot→05 · residualNorm of the residualvector norm

Script

admin

#1. Global mean sea level by satellite altimetry, mm, every 3 years since 1993;
-> "[t]" "0,3,6,9,12,15,18,21,24,27,30";
-> "[y]" "[0, 8.1, 15.0, 24.2, 33.0, 40.9, 51.2, 61.3, 73.0, 86.1, 99.4]";
#2. Design matrix of the quadratic model a + b t + c t^2;
-> "[X]" "";
for (-> "[i]" 1) (<= [i] (atom size [t] ",")) (++ "[i]") {
	-> "[ti]" (atom get [t] [i] ",");
	-> "[X]" (concat [X] (if (> [i] 1) {","} {""}) "[1," [ti] "," (math calc "[ti]^2") "]");
};
-> "[X]" (concat "[" [X] "]");
#3. Householder QR, then the triangular system R b = Q^T y;
json load "rx_p4e_qr" (matrix qr [X]);
-> "[b]" (matrix solve (json select "rx_p4e_qr" "/R") (matrix mul (matrix transpose (json select "rx_p4e_qr" "/Q")) [y]));
json unload "rx_p4e_qr";
#4. Rate b + 2 c t in 2008, acceleration 2c, residual and level in 2050, t = 57;
-> "[rate]" (vector dot [b] "[0, 1, 30]");
-> "[acc]" (vector dot [b] "[0, 0, 2]");
-> "[res]" (vector norm (vector sub (matrix mul [X] [b]) [y]));
concat "coefficients a, b, c = " [b] [_n_] "rate in 2008 = " (math round [rate] 2) " mm/yr, acceleration = " (math round [acc] 4) " mm/yr2, residual norm = " (math round [res] 2) " mm, level in 2050 = " (math round (vector dot [b] "[1, 57, 3249]") 0) " mm above 1993"

Result

mentdb

coefficients a, b, c = [0.9055944055944333,2.138352758352755,0.037347837347837434]
rate in 2008 = 3.26 mm/yr, acceleration = 0.0747 mm/yr2, residual norm = 2.37 mm, level in 2050 = 244 mm above 1993

Go further

Compare with the normal equations X^T X b = X^T y: matrix cond shows why QR is preferred.
Add a sine of period 1 year to the design matrix to remove the seasonal cycle of monthly data.

executed on MentDB

Planck feedback and climate sensitivity

The script extracts the Stefan-Boltzmann constant from the physics library, gets the effective emission temperature of the Earth with a fourth root, the Planck feedback 4 sigma Te^3, the net feedback after water vapour, lapse rate, albedo and cloud feedbacks, then the forcing of a CO2 doubling, the equilibrium climate sensitivity and the warming committed at 420 ppm.

Pipeline

01 · constantStefan-Boltzmannconstant physics_stefan_boltzmann→02 · extractValue out of the textstring regexp_extract→03 · temperatureFourth rootmath nroot→04 · forcing5.35 ln C/C0math log

Script

admin

#1. Stefan-Boltzmann constant of the physics library, value extracted from the text;
-> "[sigma]" (string regexp_extract (constant physics_stefan_boltzmann) "[0-9.]+(E-?[0-9]+)?" 0);
#2. Effective emission temperature Te = fourth root of S 1 - albedo / 4 sigma;
-> "[Te]" (math nroot (math calc "1361*(1-0.30)/(4*[sigma])") 4);
#3. Planck feedback 4 sigma Te^3 and the other feedbacks: water vapour, lapse rate, albedo, clouds, W/m2/K;
-> "[l0]" (math calc "4*[sigma]*[Te]^3");
-> "[lambda]" (math calc "[l0]-(1.8-0.5+0.35+0.42)");
#4. Forcing of a CO2 doubling and of 420 ppm against 280 ppm, 5.35 ln C/C0;
-> "[F2x]" (math calc "5.35*ln(2)");
-> "[F]" (math calc (concat "5.35*" (math log (math calc "420/280"))));
concat "sigma = " [sigma] ", Te = " (math round [Te] 2) " K, Planck feedback = " (math round [l0] 3) " W/m2/K, net feedback = " (math round [lambda] 3) " W/m2/K, F2x = " (math round [F2x] 2) " W/m2, ECS = " (math round (math calc "[F2x]/[lambda]") 2) " K, equilibrium warming at 420 ppm = " (math round (math calc "[F]/[lambda]") 2) " K"

Result

mentdb

sigma = 5.6703744191844314E-8, Te = 254.58 K, Planck feedback = 3.742 W/m2/K, net feedback = 1.672 W/m2/K, F2x = 3.71 W/m2, ECS = 2.22 K, equilibrium warming at 420 ppm = 1.3 K

Go further

Draw the warming against the cloud feedback from 0 to 1 W/m2/K with plot function: the tail of high sensitivities.
Add the ocean heat uptake to get the transient response with num ode.

executed on MentDB

Hot summers: shift of the mean and return period

Ten recent summer anomalies of a station are ranked, tested against the 1961-1990 climate with a z-test (the climatological standard deviation is known), and the return period of a +2.1 K summer is computed in the old and in the current climate with the complementary error function: the risk ratio of an attribution study.

Pipeline

01 · rankOrdered anomaliesmath order_double→02 · testz-test, known sigmastat ztest→03 · tailNormal tail 0.5 erfcmath erfc→04 · ratioReturn periodsmath calc

Script

admin

#1. Summer JJA temperature anomalies of a station against 1961-1990, last 10 years, known sigma 0.6 K;
-> "[jja]" "[0.9, 1.4, 0.7, 2.1, 1.1, 1.6, 0.8, 2.4, 1.3, 1.9]";
-> "[sigma]" 0.6;
#2. Ranking of the summers, hottest last;
-> "[rank]" (math order_double [jja]);
#3. z-test: has the mean shifted from 0?;
-> "[z]" (stat ztest [jja] 0 [sigma]);
#4. Probability of a summer above +2.1 K in the old climate and now, 0.5 erfc x / sigma sqrt 2;
-> "[p_old]" (math calc (concat "0.5*" (math erfc (math calc "2.1/([sigma]*sqrt(2))"))));
-> "[p_new]" (math calc (concat "0.5*" (math erfc (math calc (concat "(2.1-" (stat mean [jja]) ")/([sigma]*sqrt(2))")))));
concat "ordered = " [rank] [_n_] "z-test = " [z] [_n_] "summer above +2.1 K: return period " (math round (math calc "1/[p_old]") 0) " years in 1961-1990, " (math round (math calc "1/[p_new]") 1) " years now, risk ratio = " (math round (math calc "[p_new]/[p_old]") 0)

Result

mentdb

ordered = [0.7,0.8,0.9,1.1,1.3,1.4,1.6,1.9,2.1,2.4]
z-test = {"z": 7.484057129065165, "p_value": 7.20626055562942E-14}
summer above +2.1 K: return period 4299 years in 1961-1990, 7.8 years now, risk ratio = 553

Go further

Use stat ttest if the standard deviation must be estimated from the sample.
Fit a GEV law to the yearly maxima instead of a normal law on the summer means.

executed on MentDB

Amplitude and phase of the annual temperature cycle

The twelve monthly normals of a station are projected on the first harmonic. The script gets the annual mean, the amplitude with hypot, the phase with atan2, the date of the warmest day of the fitted cycle and its lag after the summer solstice, a signature of the thermal inertia of the site.

Pipeline

01 · harmonicProjection on cos and sinmath calc→02 · amplitudeAmplitudemath hypot→03 · phasePhasemath atan2→04 · dateWarmest day and lagdate add→05 · lagDays after the solsticedate nb_day

Script

admin

#1. Monthly mean temperatures of a mid-latitude station, degC, January to December;
-> "[T]" "2.1,3.4,6.9,10.2,14.1,17.6,19.8,19.3,15.9,11.4,6.3,2.9";
-> "[a]" 0; -> "[b]" 0;
#2. First harmonic: a1 and b1 by projection on cos and sin at mid-month;
for (-> "[m]" 1) (<= [m] 12) (++ "[m]") {
	-> "[x]" (math calc "2*pi*([m]-0.5)/12");
	-> "[Tm]" (atom get [T] [m] ",");
	-> "[a]" (math calc "[a]+[Tm]*cos([x])/6");
	-> "[b]" (math calc "[b]+[Tm]*sin([x])/6");
};
#3. Amplitude and phase of the annual cycle, date of the maximum;
-> "[amp]" (math hypot [a] [b]);
-> "[phase]" (math calc (concat "mod(" (math atan2 [b] [a]) "+2*pi, 2*pi)"));
-> "[day]" (math round (math calc "[phase]/(2*pi)*365") 0);
#4. Annual mean and date of the warmest day of the fitted cycle;
concat "annual mean = " (math round (math avg (concat "[" [T] "]")) 2) " degC, amplitude = " (math round [amp] 2) " K, phase = " (math round (math rad_to_deg [phase]) 1) " deg, warmest day of the fit = " (date add "2025-01-01" "DAY" [day]) ", lag after the solstice = " (date nb_day "2025-06-21" (date add "2025-01-01" "DAY" [day])) " days"

Result

mentdb

annual mean = 10.83 degC, amplitude = 8.78 K, phase = 198.1 deg, warmest day of the fit = 2025-07-21, lag after the solstice = 30 days

Go further

Compare a coastal and a continental station: the coast has a smaller amplitude and a longer lag.
Get the same harmonic with dsp spectrum on several years of daily data.

executed on MentDB

Dates of the remaining 1.5 C carbon budget

From the remaining carbon budget for 1.5 C and the current CO2 emissions, the script gets the date when the budget is used at constant emissions, the net zero date of a straight line path that stays within the budget, the yearly cut it needs, and the warming that the budget represents with the TCRE.

Pipeline

01 · yearsYears of budgetmath calc→02 · dateExhaustion datedate add→03 · splitYears and monthsmath frac_part→04 · yearNet zero yeardate year→05 · daysWhole daysmath trunc

Script

admin

#1. Remaining carbon budget for 1.5 C, 50 percent chance, 130 GtCO2 from 1 January 2025, emissions 41.6 GtCO2/yr;
-> "[B]" 130; -> "[E]" 41.6;
#2. Constant emissions: budget gone after B/E years, converted to a date;
-> "[years]" (math calc "[B]/[E]");
-> "[gone]" (date add "2025-01-01" "DAY" (math trunc (math calc "[years]*365.25")));
#3. Straight line to net zero within the budget: T = 2B/E;
-> "[T]" (math calc "2*[B]/[E]");
-> "[nz]" (date add "2025-01-01" "DAY" (math trunc (math calc "[T]*365.25")));
#4. Warming from the TCRE, 0.45 C per 1000 GtCO2, and the yearly cut needed;
concat "budget gone on " [gone] " at constant emissions, " (math trunc [years]) " years and " (math round (math calc (concat "12*" (math frac_part [years]))) 1) " months; linear path: net zero on " [nz] " in year " (date year [nz]) ", cut of " (math round (math calc "[E]/[T]") 2) " GtCO2 per year; warming of the budget " (math round (math calc "0.45*[B]/1000") 3) " C"

Result

mentdb

budget gone on 2028-02-16 at constant emissions, 3 years and 1.5 months; linear path: net zero on 2031-04-02 in year 2031, cut of 6.66 GtCO2 per year; warming of the budget 0.059 C

Go further