A simple theoretical model for lags and asymmetries of surface temperature
Author
Gabriele Di Bona, Andrea Giacobbe
Title
A simple theoretical model for lags and asymmetries of surface temperature
Description
Creates datasets of average temperature and humidity and runs an ODE that simulates their evolution.
Category
Academic Articles & Supplements
Keywords
nonlinear dynamical systems, climate modelling, local climates, seasonal lag, diurnal lag, asymmetric evolution of daily temperatures, Earth and planetary climate
URL
http://www.notebookarchive.org/2021-05-1vzle31/
DOI
https://notebookarchive.org/2021-05-1vzle31
Date Added
2021-05-04
Date Last Modified
2021-05-04
File Size
106.36 kilobytes
Supplements
Rights
Redistribution rights reserved
Download
Open in Wolfram Cloud
A simple theoretical model for lags and asymmetries of surface temperature
A simple theoretical model for lags and asymmetries of surface temperature
The information in this notebook is a supplement for "A Simple Theoretical Model for Lags and Asymmetries of Surface Temperature" by Gabriele Di Bona and Andrea Giacobbe
Generation of the dataset
Generation of the dataset
It is well known that Meteorology and Climatology are different subjects, and manuals roughly define Climate as the average over a period of 30-50 years of meteorological phenomena. Let us put to work this definition.
In what follows, once chosen a particular WMO station (that typically is labelled with WMO followed by a five-digits number but that it can also be labelled as a ICAO code if it is situated in an airport) we produce timeseries in which time is time is expressed in days in UTC, starting from January 1st, 00:00 and ending in December 31st, 24:00 and the values are the temperature or the humidity measured by the given station.
The time step of the dataseries depends on the choice of a variable called bin is chosen automatically based on the number of available data by giving to a variable called bin the average value of data per day.
Les us set the working directory as the directory in which the notebook resides
In[]:=
SetDirectory[NotebookDirectory[]];
Once given a code (code) a beginning year (begin) a final year (end) and a number of measures per day (bin), the following function saves in the home directory two DAT files labelled tmp_code.dat and hmd_code.dat (code is variable) relative to the average temperature and humidity. The functions also prints, for every year, the number of data per day that are being retrieved (this is strictly related to the choice of bin).
In[]:=
datasettmphmd[code_,begin_,end_,bin_]:=Module{},year=begin;InitialYear=begin;FinalYear=end;listtmp={};listhmd={};Whileyear≤end,beginningTime=AbsoluteTime[{year,01,01,00,00,00}];endingTime=AbsoluteTime[{year,12,31,23,59,59}];tmp=WeatherData[code,"Temperature",{DateString[beginningTime],DateString[endingTime]}];Printyear,Floor;hmd=WeatherData[code,"Humidity",{DateString[beginningTime],DateString[endingTime]}];listtmp=listtmp~Join~SelectTransposeN,Map[QuantityMagnitude,tmp[[2,1,1]]],NumberQ[#[[2]]]&;listhmd=listhmd~Join~SelectTransposeN,hmd[[2,1,1]],NumberQ[#[[2]]]&;year=year+1;ExportStringJoin["tmp_",code,".dat"],NMapMean,SplitByMapRound#[[1]],,#[[2]]&,Sort[listtmp],#[[1]]&;ExportStringJoin["hmd_",code,".dat"],NMapMean,SplitByMapRound#[[1]],,#[[2]]&,Sort[listhmd],#[[1]]&;
Length[tmp[[2,2,1,1]]]
365
tmp[[2,2,1,1]]-beginningTime
60*60*24
hmd[[2,2,1,1]]-beginningTime
60*60*24
1
bin
1
bin
In our investigation we did analyze in details 6 places on Earth:
Haway, Hilo airport, PHTO or WMO91285
Lybia, Kufra airport, HLKF or WMO62271
Italy, Catania airport, LICC or WMO16460
USA, Lincoln airport, KLNK or WMO72551
Antartica, Vostok station, WMO89606 (non in the database I found)
Australia, Sydney airport, YSSY or WMO94767
We write the lines that will produce the dataset of our interest
Haway, Hilo airport, PHTO or WMO91285
Lybia, Kufra airport, HLKF or WMO62271
Italy, Catania airport, LICC or WMO16460
USA, Lincoln airport, KLNK or WMO72551
Antartica, Vostok station, WMO89606 (non in the database I found)
Australia, Sydney airport, YSSY or WMO94767
We write the lines that will produce the dataset of our interest
datasettmphmd["PHTO",1973,2020,24](*hilo*)
datasettmphmd["HLKF",1973,2020,8](*kufra*)
In[]:=
datasettmphmd["LICC",1973,2020,24](*catania*)
datasettmphmd["KLNK",1973,2020,24](*lincoln*)
In[]:=
datasettmphmd["WMO89606",1973,2020,4](*vostok*)
datasettmphmd["YSSY",1973,2020,24](*sydney*)
Simulation
Simulation
The system has a certain number of fixed parameters. We define them in groups. The first thing to do is to choose the code of the weather station (the same code we used in the previous chapter).
In[]:=
code="LICC"(*tobefilled*);
Then we fix the parameters. Invariable for the Earth are the astronomical parameters
In[]:=
σB=5.670373*;(*Stefan-BoltzmannConstant*)TSun=5778;(*SuperficialtemperatureoftheSun*)RS=695700000;(*Sunradius*)dESun=149600000000;(*Earth-Sundistance*)ϵ=0.0167;(*Earthorbiteccentricity.*)δ=*2π;(*Anglebetweenperihelionandwintersolstice*)γ=Re[23.27Degree];(*obliquity*)Day=3600*24;(*averagesolardayinseconds*)Year=Day*365.242199;Ω=2π+;(*2π/siderealdayexpressedinseconds*)RE=6378136;(*equatorialEarthradius*)RP=6356751.302;(*polarEarthradius*)
-8
10
13
365
1
Day
1
Year
Variable but dependent on the location are latitude and longitude
In[]:=
{lat,long}=WeatherData[code,"Coordinates"];φ=Re[latDegree];(*Latitudeofthesystem*)ψ=Re[longDegree];(*Longitudeofthesystem*)
Variable but also dependent on the location are the geothermal parameters. They have very little influence on simulations of Earth.
In[]:=
RW=RECos[φ]+RPSin[φ]-15000;(*OceanicMohodiscontinuitydepth*)RG=RECos[φ]+RPSin[φ]-35000;(*ContinentalMohodiscontinuitydepth*)T1G=700;(*SupericialtemperatureattheMohodiscontinuityunderthecontinentalcrust*)T1W=700;(*SupericialtemperatureattheMohodiscontinuityundertheoceaniccrust*)T2=300;(*averagesuperficialtemperatureonEarth*)kW=30;(*Thermalconductivityintheoceaniccrust*)kG=30;(*Thermalconductivityinthecontinentalcrust*)qW=N-;(*Thermalpowerofenergyfromthemantleforaunitofsurfaceontheoceaniccrust*)qG=N-;(*Thermalpowerofenergyfromthemantleforaunitofsurfaceonthecontinentalcrust*)
T1W-T2
1
4πkW
1
RW
1
RECos[φ]+RPSin[φ]
1
4π
2
RW
T1G-T2
1
4πkG
1
RG
1
RECos[φ]+RPSin[φ]
1
4π
2
RG
Parameters that depend very much on the climatic type of the region under investigation are the climatic parameters
In[]:=
p=0.6;(*fractionofground*)q=1-p;(*fractionofwater*)μ=2.8*;(*evaporationratefromsoilandwater*)ν=1*;(*rainfallrate*)
-5
10
-5
10
Also somehow variable are the thermal capacities. They depend on the depth of the layer and the type of material for the soil and the presence of iced water.
In[]:=
ℓ0=400;(*depthofatmosphere*)ρ0=1.225;(*densityofair*)Cl=1*;(*heatcapacityofsoil*)C0d=711.28*ρ0*ℓ0;(*heatcapacitydryair*)C0v=2050.16*ρ0*ℓ0;(*heatcapacitywatervapor*)ℓw=40(*depthofwater*);Cw=4214.4*1000*ℓw(*heatcapacityofwater*);
6
10
Invariable are the latent heats.
In[]:=
λ=2260000;(*Latentheatofvaporizationandofcondensation(donottouch)*)Λ=λ*ρ0*ℓ0;(*Latentheatofvaporizationandofcondensation(donottouch)*)
Also variable are the optical parameters. Reflectance and emissivity depend on the substance being irradiated
In[]:=
ϵW=0.94;(*emissivitywater*)αW=0.85;(*absorbancewater*)τW=0;(*transmittancewater(donottouch)*)rW=1-αW-τW;(*reflectancewater*)ϵG=0.94;(*emissivitysoil*)αG=0.8;(*absorbancesoil*)τG=0;(*trasmittancesoil(fixed)*)rG=1-αG-τG;(*reflectancesoil*)ϵAu=0.45;(*emissivityairupwards(donottouch)*)ϵAd=0.8;(*emissivityairdownwards(donottouch)*)τA=0.52;(*trasmittanceair(donottouch)*)rA=0.23;(*reflectanceair(donottouch)*)αA=1-τA-rA;(*absorbanceair(donottouch)*)αAT=0.86;(*absorbanceinfraredair*)
Finally, variable are the linear heat transfer coefficients. They model the way in which convection exchanges heat between the materials.
In[]:=
hAG=12;(*Heattransfercoefficientair-ground*)hAW=28.5;(*Heattransfercoefficientair-water*)
These two function are functions of time. Once the position of the Earth with respect to the Sun ϑ[t] and the position of the region on Earth are known, the function SolRad gives the intensity of solar radiation that hits a square meter of upper atmosphere.
ψ+tΩ+*2π
10Day
Year
In[]:=
scalprod[t_]:=-+-+-+;SolRad[t_]:=MaxσBscalprod[t],0;
RECos[γ]Cos[δ+ϑ[t]]Cos[φ]Cosψ+tΩ+*2π
10Day
Year
2
RE
2
Cos[φ]
2
RP
2
Sin[φ]
RPCos[δ+ϑ[t]]Sin[γ]Sin[φ]
2
RE
2
Cos[φ]
2
RP
2
Sin[φ]
RECos[φ]Sin[δ+ϑ[t]]Sinψ+tΩ+*2π
10Day
Year
2
RE
2
Cos[φ]
2
RP
2
Sin[φ]
4
TSun
2
RS
2
(1+ϵCos[ϑ[t]])
2
dESun
2
(1-)
2
ϵ
The dynamical system goes to a periodic equilibrium solution. We have to give the right initial data to make sure that the solution is at thermal equilibrium. The following procedure is needed to define the initial data. The last recurrence gives directly the solution of our model. The interpolating functions are stored in the variable solution
In[]:=
TW0=273.15+5;(*guessforinitialtemperatureofwater*)TA0=273.15+4;(*guessforinitialtemperatureofair*)TG0=273.15+5;(*guessforinitialtemperatureofsoil*)U0=0.2;(*guessforinitialabsolutehumidity*)ϑ0=-2π;(*initialanglebetweenperihelionandsolstice.TheinitialtimeisonJanuary1st,theEarthisatperiheliononJanuary4th*)system=(C0d+C0vU[t])TA'[t]αA(1+pτArG+qτArW)SolRad[t]+αAT*σB*(pϵG+qϵW)-σB(ϵAu+ϵAd)+phAG(TG[t]-TA[t])+qhAW(TW[t]-TA[t])-Λ*Min[μ(-U[t])-νU[t],0],ClTG'[t]τAαGSolRad[t]+σB(ϵAd-ϵG)-hAG(TG[t]-TA[t])+qG-ΛMax[μ(-U[t])-νU[t],0],CwTW'[t]τAαWSolRad[t]+σB(ϵAd-ϵW)-hAW(TW[t]-TA[t])+qW-Λ*Max[μ(-U[t])-νU[t],0],U'[t]μ(-U[t])-νU[t],ϑ'[t];solution=NDSolve[{system,TG[0]TG0,TW[0]TW0,TA[0]TA0,U[0]U0,ϑ[0]ϑ0},{TG,TW,TA,U,ϑ},{t,0,1.01Year}];While[Abs[TA0-First[TA/.solution][365*Day]]>||Abs[TG0-First[TG/.solution][365*Day]]>||Abs[TW0-First[TW/.solution][365*Day]]>||Abs[U0-First[U/.solution][365*Day]]>,TA1=First[TA/.solution][365*Day];TG1=First[TG/.solution][365*Day];TW1=First[TW/.solution][365*Day];U1=First[U/.solution][365*Day];solution=NDSolve[{system,TG[0]TG1,TW[0]TW1,TA[0]TA1,U[0]U1,ϑ[0]ϑ0},{TG,TW,TA,U,ϑ},{t,0,1.01Year}];TA0=TA1;TG0=TG1;TW0=TW1;U0=U1;];
3
365
4
(TG[t])
4
(TW[t])
4
(TA[t])
0.0666TA[t]-23.96
E
4
(TA[t])
4
(TG[t])
0.0666TA[t]-23.96
E
4
(TA[t])
4
(TW[t])
0.0666TA[t]-23.96
E
0.0666TA[t]-23.96
E
2π
Year
3
Sqrt[1-]
2
ϵ
2
(1+ϵCos[ϑ[t]])
-1
10
-1
10
-1
10
-3
10
Visualization
Visualization
Yearly temperature and humidity
Yearly temperature and humidity
The following lines plot the real averaged temperature and relative humidity together with the ones computed with the model
In[]:=
tottmp=Import[StringJoin["tmp_",code,".dat"]];tothmd=Import[StringJoin["hmd_",code,".dat"]];Tmp=Evaluate@Interpolation[tottmp,InterpolationOrder1];Hmdr=Evaluate@Interpolation[tothmd,InterpolationOrder1];US[T_]:=;Hmd[t_]:=Hmdr[t]US[Tmp[t]+273.15]realTmp=Plot[Tmp[x],{x,0,365},PlotPoints1000,PlotStyle{Black,Opacity[0.4]}];realHmdr=Plot[Hmdr[x],{x,0,365},PlotPoints1000,PlotStyle{Black,Opacity[0.5]}];solstices={{79,{Dashed,Green}},{172,{Dashed,Darker[Orange]}},{266,{Dashed,Brown}},{355,{Dashed,Cyan}}};Plot[{First[TA/.solution][tDay]-273.15},{t,0,365},PlotStyle{Red,Opacity[0.8]},PlotPoints5000,ImageSize400,PlotRangeAll];Show[{%,realTmp},GridLines{solstices,None},GridLinesStyleThickness[0.003],PlotRangeAll,AxesStyleArrowheads[0.03],AxesLabel{"days","C"},LabelStyleLarge,ImageSize700]Plot,{t,0,365},PlotStyle{Blue,Opacity[0.6]},PlotPoints3000,ImageSize400,PlotRange{0,1};Show[{%,realHmdr},PlotRange{0,1},GridLines{solstices,None},GridLinesStyleThickness[0.003],AxesStyleArrowheads[0.03],AxesLabel{"days",""},LabelStyleLarge,ImageSize700]
0.0666T-23.96
E
∘
First[U/.solution][tDay]
US[First[TA/.solution][tDay]]
Daily temperature and humidity
Daily temperature and humidity
The following lines, once given an integer from 0 to 365, plot the real averaged temperature and relative humidity together with the ones computed with the model in a period of three days starting form the chosen day. The spring equinox takes place at day 79, the summer solstice takes place at day 172, the autumn equinox takes place at day 266, the winter solstice takes place at day 355.
In[]:=
day=150;col[i_]:={Dashed,If[EvenQ[i],Blue,Darker[Orange]]};threedaysTmp=Plot[Tmp[x],{x,day,day+3},PlotPoints1000,PlotStyle{Thickness[0.003],Black}];threedaysHmdr=Plot[Hmdr[x],{x,day,day+3},PlotPoints1000,PlotStyle{Thickness[0.003],Black}];threedaysTmpsim=PlotFirst[TA/.solution][tDay]-273.15,{t,day,day+3},GridLinesTable-,col[i],{i,0,{t,day,day+3},GridLinesTable-,col[i],{i,0C"},LabelStyleLarge,ImageSize500],Show[{threedaysHmdrsim,threedaysHmdr},PlotRange{0,1},AxesStyleArrowheads[0.04],GridLinesStyleThickness[0.004],Ticks{Table[i,{i,0,400}],{0.2,0.4,0.6,0.8,1}},AxesLabel{"days",""},LabelStyleLarge,ImageSize500]}]
i
2
long
360
,
365},None,PlotStyle{Thickness[0.004],Red},PlotPoints3000,ImageSize400,PlotRange{-95,95};threedaysHmdrsim=PlotFirst[U/.solution][tDay]
US[First[TA/.solution][tDay]]
i
2
long
360
,
2*365},None,PlotStyle{Thickness[0.004],Blue},PlotPoints3000,ImageSize400,PlotRange{0,1};Row[{Show[{threedaysTmpsim,threedaysTmp},PlotRange{0,35},AxesStyleArrowheads[0.04],GridLinesStyleThickness[0.004],Ticks{Table[i,{i,0,400}],Automatic},AxesLabel{"days","∘
We create a spline of the average temperatures for every day of the year. We triplicate the list of temperatures and choose the central year to avoid aberrations at the beginning or at the end of the graph due to the spline algorithm.
Analysis
Analysis
L1 distance
L1 distance
To determine the accuracy of the simulated temperature, we define a function that computes a discretized L1-distance. On average the real temperature and the simulated one differ by the amount indicated by the output
In[]:=
1
Length[tottmp]
This line determines the accuracy of the simulated relative humidity. On average the real relative humidity and the simulated one differ by the amount indicated by the output
In[]:=
1
Length[tothmd]
First[U/.solution][Daytothmd[[ind,1]]]
US[First[TA/.solution][Daytothmd[[ind,1]]]]
Lag of season
Lag of season
In[]:=
Needs["Splines`"]tmpYmean=Transpose[Map[Mean,SplitBy[Map[{Floor[#[[1]]],#[[2]]}&,tottmp],#[[1]]&]]][[2]];tmpYmeanlong=Join[tmpYmean,tmpYmean,tmpYmean];tmpYmeanspline=BezierFunction[tmpYmeanlong];
We do the same for the amount of daily energy. This requires the computation of 365 numerical integration which is long.
In[]:=
solenerg=Flatten[Table[NIntegrate[SolRad[tDay]/.solution,{t,d,d+1}],{d,0,364}]];solenerglong=Join[solenerg,solenerg,solenerg];solenergspline=BezierFunction[solenerglong];
We hence compute the maximum of solar energy and the maximum of temperature during the year. The algorithm requires and initial guess that should possibly be changed.
In[]:=
guessmaxenerg=180;guessmaxtmp=220;maxenerg=FindArgMaxsolenergspline,{t,guessmaxenerg}[[1]];maxtmp=FindArgMaxtmpYmeanspline,{t,guessmaxtmp}[[1]];
t+365
3*365
t+365
3*365
We finally plot the graphs and the maxima determined above.
In[]:=
TwoAxisPlot[{f_,g_},{x_,x1_,x2_},{a_,b_,c_,d_}]:=Module[{fgraph,ggraph,frange,grange,fticks,gticks},{fgraph,ggraph}={Plot[f,{x,x1,x2},AxesTrue,PlotStyleDarker[Orange]],Plot[g,{x,x1,x2},AxesTrue,PlotStyleBlack]};{frange,grange}=(PlotRange/.AbsoluteOptions[#,PlotRange])[[2]]&/@{fgraph,ggraph};fticks=N@FindDivisions[frange,5];gticks=Quiet@Transpose@{fticks,ToString[NumberForm[#,3],StandardForm]&/@Rescale[fticks,frange,grange]};Show[fgraph,ggraph/.Graphics[graph_,s___]Graphics[GeometricTransformation[graph,RescalingTransform[{{0,1},grange},{{0,1},frange}]],s],AxesFalse,FrameTrue,LabelStyleLarge,FrameStyle{{Darker[Orange],Black},{Black,Black}},FrameTicks{{fticks,gticks},{Automatic,{{Round[maxenerg,1],Style[Round[maxenerg,1],Darker[Orange]]},Round[maxtmp,1]}}},RotateLabelFalse,FrameLabel{a,b,c,d},ImageSize750]]TwoAxisPlotsolenergspline,tmpYmeanspline,{t,0,365},"days","","","C"
t+365
3*365
t+365
3*365
KWh
2
m
∘
We do the same with the simulated temperatures
In[]:=
tmpYmeansim=TableNIntegrate[First[TA/.solution][t],{t,dayDay,(day+1)Day}]-273.15,{day,0,364};tmpYmeansimlong=Join[tmpYmeansim,tmpYmeansim,tmpYmeansim];tmpYmeansimspline=BezierFunction[tmpYmeansimlong];
1
Day
In[]:=
maxtmpsim=FindArgMaxtmpYmeansimspline,{t,220}[[1]];figlagb=PlottmpYmeanspline,tmpYmeansimspline,{t,0,365},AxesFalse,FrameTrue,PlotStyle{Black,Red},FrameStyle{{Black,Black},{Black,Black}},RotateLabelFalse,FrameLabel{"days","","","C"},FrameTicks{{None,All},{Table[50i,{i,0,8}],{{Round[maxtmpsim,1],Style[Round[maxtmpsim,1],Red]}}}},LabelStyleLarge,ImageSize650
t+365
3*365
t+365
3*365
t+365
3*365
∘
Lag of noon
Lag of noon
The following line computes the average time of the day in which the radiation from the sun is strongest. The evaluation of these cells requires a long time.
In[]:=
MaxSolRadDay=TableArgMax[{SolRad[t]/.solution[[1]],t>Day*date,t<Day*(date+1)},t],{date,0,364,1};MaxSolRadDayMOD=Mean[Mod[24*(MaxSolRadDay-MaxSolRadDay[[1]])+12,24]]-12+24MaxSolRadDay[[1]]
1
Day
The following line computes the average time of the day in which the real temperature is the highest.
In[]:=
splitted=SplitBy[tottmp,Floor[#[[1]]]&];maxDayReal=Table[splitted[[iter+1]][[Position[Transpose[splitted[[iter+1]]][[2]],Max[Transpose[splitted[[iter+1]]][[2]]]][[1,1]]]],{iter,0,364}];maxTempRealDay=Transpose[maxDayReal][[1]];maxTempRealDayMOD=Mean[Mod[24*(maxTempRealDay-maxTempRealDay[[1]])+12,24]]-12+24maxTempRealDay[[1]]
The following line computes the average time of the day in which the simulated temperature is highest. The evaluation of these cells requires a long time.
In[]:=
MaxTempSimDay=TableArgMax[{First[TA/.solution][t],t>Day*date,t<Day*(date+1)},t],{date,0,364,1};MaxTempSimDayMOD=Mean[Mod[24*(MaxTempSimDay-MaxTempSimDay[[1]])+12,24]]-12+24MaxTempSimDay[[1]]
1
Day
The real and the simulated lag of noons are hence
In[]:=
TimeObject[maxTempRealDayMOD]-TimeObject[MaxSolRadDayMOD]TimeObject[MaxTempSimDayMOD]-TimeObject[MaxSolRadDayMOD]
Cite this as: Gabriele Di Bona, Andrea Giacobbe, "A simple theoretical model for lags and asymmetries of surface temperature" from the Notebook Archive (2021), https://notebookarchive.org/2021-05-1vzle31
Download
